Niche overlap and Hopfield-like interactions in generalized random Lotka-Volterra systems

I. INTRODUCTION

The foundations of the theory of disordered systems date back to the early 1970s . Initially, the aim was to understand certain magnetic states in condensed-matter physics (“spin glasses”) . However, it became clear that applications of the tools developed for disordered systems had a reach far beyond the boundaries of physics. Methods such as replica theory or dynamic generating functionals were quickly adapted and used to answer questions in neural networks and to study the Minority Game (sometimes presented as a simple model of a financial market) or indeed evolutionary bimatrix games and so-called Nash equilibria .

The defining feature of disordered systems is the presence of quenched disorder. That is, the system is made up of many constituents, and the interactions between these are determined by coefficients that are drawn at random at the beginning but then remain fixed as the dynamics of the system unfolds. The disorder leads to complicated energy landscapes. The number of local minima can grow exponentially in the size of the system and is often organized in a hierarchical manner. Dynamic phenomena in disordered systems include ergodicity breaking and so-called ageing .

Ideas and methods from the physics of disordered systems have also been used to study complex ecosystems . The word “complex” in this context indicates that the ecosystem is composed of a large number of species and that these species are subject to randomly drawn interaction coefficients. In this paper we continue this line of work and focus on a Lotka-Volterra system with Hopfield-like interactions . More specifically, we are interested in a set of $𝑁$ species ( $𝑁 ≫ 1$ ), whose abundances develop in time following a generalized Lotka-Volterra equation (details will follow in Sec. ). This involves an $𝑁 \times 𝑁$ matrix $𝑎 𝑖 𝑗$ of interaction coefficients. Existing work on the statistical physics of complex ecosystems has mostly focused on the case in which the interaction matrix is drawn from distributions with either no correlations between different matrix elements or only correlations between diagonally opposed entries . There is also work on cases in which the matrix is composed of blocks and where the elements in different blocks have different statistics . One common element shared by many existing random Lotka-Volterra models is that the finest level of modeling is set by the interaction coefficients. No further assumptions are made about the properties of the species and how the species interactions come about from these properties.

The Hopfield model is inspired by structures first used in neural networks . Translated into the language of ecology, the starting point is now a set of species and a set of traits. Each species can either possess or not possess a given trait. This assignment of traits to species, in turn, determines how species will interact. Broadly speaking, the interaction between two species will be more competitive the more traits they share (i.e., the more similar the two species are). This type of interaction structure has also been studied in models combining resources and consumers, both in economics and in ecology . A particularly notable model is that by MacArthur and collaborators . Analyses of random replicator systems with “Hebbian” interactions have shown interesting statistical mechanics and in particular types of phase transition that are different from those seen in replicator systems with Gaussian couplings.

In this paper, we set out to characterize the behavior of a Lotka-Volterra system with Hopfield-like interactions, where we allow for a degree of “mild” dilution (the system is not fully connected, but each species still interacts with an extensive number of other species). A system of replicator equations with such interactions was studied in Ref. . Our aim is to calculate the statistics of fixed points in the phase where such fixed points are attained and identify the onset of instability. As in the system with Gaussian interactions, we find that only a proportion of the initial species survive at stable fixed points. Recent work on Gaussian systems has shown that the reduced interaction matrix (the matrix of interaction coefficients among the surviving species) has intricate statistics. Specifically, its bulk and outlier eigenvalues can be related to different types of dynamic phase transitions. As we will show, the types of phase transition seen in our model differ from those in the Gaussian model. One aim of the current paper is therefore to establish (in simulations) how these transitions relate to the spectra of the interaction matrix of the extant species.

The remainder of the paper is organized as follows. In Sec. we define the model and introduce the necessary notation. Section then contains the mathematical analysis. This is based on so-called generating functionals and dynamic mean-field theory. The phase diagram and further behavior of the model are then discussed in Sec. . In Sec. we finally turn to a study of the spectra of the reduced interaction matrix and their relation to the phase diagram. We conclude the paper with a discussion and an outlook in Sec. .

II. MODEL DEFINITIONS

We will study the following generalized Lotka-Volterra equation (gLVE)

̇ 𝑥 𝑖 (𝑡) = 𝑥 𝑖 (𝑡) [𝐾 𝑖 - 𝑢 𝑖 𝑥 𝑖 + \sum 𝑗 \neq 𝑖 𝑐 𝑖 𝑗 𝐽 𝑖 𝑗 𝑥 𝑗],

(1)

where the $𝑥 𝑖 \geq 0$ represent the abundances (or population densities) of different species, $𝑖 = 1, \dots 𝑁$ . We always assume initial conditions for which all $𝑥 𝑖$ are strictly positive.

The quantities $𝑢 𝑖 > 0$ denote the strength of intraspecific competition, and the $𝑎 𝑖 𝑗 = 𝑐 𝑖 𝑗 𝐽 𝑖 𝑗$ represent the interspecific interactions. The $𝐾 𝑖$ (together with the $𝑢 𝑖$ ) set the carrying capacities of the species in the absence of interactions between different species ( $𝑥 𝑖$ then tends to $𝐾 𝑖 / 𝑢$ in the long run). We focus on the case $𝑢 𝑖 \equiv 𝑢$ for all $𝑖$ , noting that $𝑢$ controls the time scale on which the noninteracting system approaches the fixed point $𝑥 𝑖 \equiv 𝐾 𝑖 / 𝑢$ . We allow for general positive values of $𝑢$ throughout our analysis, but in an effort to keep the number of parameters manageable we set $𝐾 𝑖 \equiv 1$ .

The dilution variables $𝑐 𝑖 𝑗 \in {0, 1}$ ( $𝑖 \neq 𝑗$ ) determine which species interact with one another, i.e., they set the topology of the interaction network. For each pair $𝑖 < 𝑗$ , the coefficients $𝑐 𝑖 𝑗$ and $𝑐 𝑗 𝑖$ are chosen from a Bernoulli distribution with

⟨ 𝑐 𝑖 𝑗 ⟩ = 𝑐, ⟨ 𝑐 𝑖 𝑗 𝑐 𝑗 𝑖 ⟩ - 𝑐 2 = Γ 𝑐 (1 - 𝑐) .

(2)

We thus have $𝑃 (𝑐 𝑖 𝑗 = 1) = 𝑐$ for all $𝑖 \neq 𝑗$ , i.e., $𝑐$ is the analog of what May called “connectance” . The parameter $Γ$ is restricted to the range from $- 1$ to 1 by construction, but we note that not all choices of pairs $(𝑐, Γ)$ are possible (see the Supplemental Material (SM) for details). We note that the choice of the diagonal coefficients $𝑐 𝑖 𝑖$ is irrelevant as we set $𝐽 𝑖 𝑖 = 0$ below

Throughout our paper, $𝑐$ is chosen not to scale with $𝑁$ $[𝑐 = 𝒪 (𝑁 0)]$ . This means that each species interacts with an $𝒪 (𝑁)$ number of other species. We are therefore not studying a “dilute” system in the sense of random matrix theory. The extensive connectivity allows us to use established methods from dynamic mean-field theory. Truly dilute systems with $𝑐 = 𝒪 (1 / 𝑁)$ (and where, consequently, each species only interacts with a finite number of other species) can be expected to behave very differently, see, e.g., Ref. , and a theoretical analysis would be much more intricate.

Interaction links in our system are directed, that is, an effect of the presence of species $𝑗$ on the dynamics of $𝑖$ does not necessarily imply the reverse. The parameter $Γ$ measures the correlations between $𝑐 𝑖 𝑗$ and $𝑐 𝑗 𝑖$ . A choice of $Γ = - 1$ implies $𝑐 𝑖 𝑗 = 1 - 𝑐 𝑗 𝑖$ with probability 1, and $Γ = 1$ means that $𝑐 𝑖 𝑗 = 𝑐 𝑗 𝑖$ with probability 1.

The matrix $𝐽 𝑖 𝑗$ determines the strength of the effects of the presence of species $𝑗$ on the dynamics of the abundance of species $𝑖$ . Positive values of $𝐽 𝑖 𝑗$ imply that the population of species $𝑗$ is beneficial to the growth of species $𝑖$ , while negative values imply a detrimental interaction. In ecological terms the signs of the pair of interactions $(sgn 𝐽 𝑖 𝑗, sgn 𝐽 𝑗 𝑖)$ determine whether two species are in a mutualistic relation $(+, +)$ , whether they compete with one another $(-, -)$ , or whether there is an antagonistic predator-prey relation between $𝑖$ and $𝑗 (\pm, \mp)$ .

In this work, the interspecific interaction is chosen according to a niche overlap heuristic (see e.g., Ref. ). We assume each species is described by a set of binary traits, labeled $𝜇 = 1, \dots, 𝑃$ , and that a pair of species will compete in proportion to the similarity between the two species (i.e., the number of traits which both or neither species possess). We write $𝜉 𝜇 𝑖 = + 1$ , if species $𝑖$ has trait $𝜇$ , and $𝜉 𝜇 𝑖 = - 1$ if the species does not possess the trait. Interactions are then assumed to be of the form

𝐽 𝑖 𝑗 = {- 1 𝑐 𝑁 \sum 𝛼 𝑐 𝑁 𝜇 = 1 𝜉 𝜇 𝑖 𝜉 𝜇 𝑗 𝑖 \neq 𝑗 0 𝑖 = 𝑗 .

(3)

We have here set $𝑃 = 𝛼 𝑐 𝑁$ , with $𝛼 > 0$ a model parameter (in simulations $𝑃$ is restricted to integer values.) That is, we assume that the number of traits is proportional to the number of species in the system. The interaction in Eq. is reminiscent of the Hopfield model, used in the context of neural networks . This suggests interesting phase behavior when $𝛼 = 𝒪 (𝑁 0)$ , which is the regime we focus on. We have normalized the interaction strength by $𝑐 𝑁$ , the mean number of species that any one species will interact with. We will refer to the random variables $𝑐 𝑖 𝑗$ and $𝜉 𝜇 𝑖$ as the disorder of the system.

The traits $𝜉 𝜇 𝑖$ are chosen to be $𝜉 𝜇 𝑖 = \pm 1$ with equal probability, and there is no correlation between the different $𝜉 𝜇 𝑖$ . This implies that the distribution of the $𝐽 𝑖 𝑗$ approaches a Gaussian as $𝑁 \to \infty$ , reminiscent of the model studied for example in Refs. . We note, however, that the Hopfield structure introduces correlations between the different $𝐽 𝑖 𝑗$ , which are different from the correlations studied in the earlier literature. We highlight again the structural similarity to MacArthur's consumer-resource model (see also Refs. for statistical physics studies), noting though that the latter model is more sophisticated, with dynamical equations both for consumers and resources.

III. GENERATING FUNCTIONAL ANALYSIS AND STABILITY

A. Generating functional and effective process

We analyze the system in Eq. using dynamic generating functionals, an established method in the theory of disordered systems . This leads to an effective “dynamic mean-field theory.” Similar approaches have been used to study Lotka-Volterra systems with Gaussian random couplings . We note that an alternative approach is based on the so-called cavity method . We also add that the dynamics admit a Liapunov function when $𝑐 = 1$ (leading to a symmetric interaction matrix). Methods from the equilibrium statistical physics of disordered systems such as the replica approach can be used in this special case (for examples see Refs. ).

The outcome of the application of these techniques is an effective stochastic process for a “representative species.” The ensemble of realizations of stochastic processes is statistically equivalent to the set of single-species trajectories $𝑥 𝑖 (𝑡)$ of the disordered dynamical system in Eq. . The dynamic mean-field description becomes exact in the thermodynamic limit ( $𝑁 \to \infty$ ). Overall, in this limit, the infinite-dimensional deterministic dynamical system in Eq. is traded for an effective single-species process which is nonlocal in time (it involves retarded self-interaction) and contains colored noise.

The generating functional analysis begins from

̇ 𝑥 𝑖 (𝑡) = 𝑥 𝑖 (𝑡) [1 - 𝑢 𝑥 𝑖 + \sum 𝑗 \neq 𝑖 𝑐 𝑖 𝑗 𝐽 𝑖 𝑗 𝑥 𝑗 - ℎ 𝑖 (𝑡)],

(4)

where we have introduced the perturbation fields $ℎ 𝑖 (𝑡)$ in order to calculate linear response functions. These fields are not actually part of the model and are set to zero at the end of the calculation, as well as in all simulations shown in the paper. For more details see the SM. The generating functional of this dynamical system is given by

𝒵 [𝜓 𝑖 (𝑡), ℎ 𝑖 (𝑡)] = ⟨ e x p [𝑖 \sum 𝑖 \int 𝑑 𝑡 𝑥 𝑖 (𝑡) 𝜓 𝑖 (𝑡)] ⟩ paths,

(5)

where the average is over paths $[𝑥 1 (𝑡), \dots, 𝑥 𝑁 (𝑡)]$ of the dynamics in Eq. . The $𝜓 𝑖 (𝑡)$ constitute a source field. The generating functional in Eq. is the Fourier transform of the probability measure in the space of paths generated by Eq. .

The final outcome of the generating-functional analysis is a set of equations for the dynamic macroscopic order parameters of the problem. For the Lotka-Volterra model these are

𝑀 (𝑡) = l i m 𝑁 \to \infty 1 𝑁 𝑁 \sum 𝑖 = 1 ‾ ⟨ 𝑥 𝑖 (𝑡) ⟩ 0, 𝐶 (𝑡, 𝑡') = l i m 𝑁 \to \infty 1 𝑁 𝑁 \sum 𝑖 = 1 ‾ ⟨ 𝑥 𝑖 (𝑡) 𝑥 𝑖 (𝑡') ⟩ 0, 𝐺 (𝑡, 𝑡') = l i m 𝑁 \to \infty 1 𝑁 𝑁 \sum 𝑖 = 1 𝛿 ‾ ⟨ 𝑥 𝑖 (𝑡) ⟩ 0 𝛿 ℎ 𝑖 (𝑡'),

(6)

where $𝛿$ denotes a functional derivative and $⟨ \dots ⟩ 0$ stands for an average over random initial conditions. The overbar $‾ \dots$ represents the average over the disorder, i.e., over the $𝑐 𝑖 𝑗$ and $𝜉 𝜇 𝑖$ . The order parameters can be obtained from the disorder-averaged generating functional as derivatives with respect to the fields $𝜓 𝑖 (𝑡)$ and/or $ℎ 𝑖 (𝑡)$ , evaluated at $𝜓 𝑖 (𝑡) \equiv 0$ and $ℎ 𝑖 (𝑡) \equiv 0$ .

The order parameters in Eqs. are determined self-consistently from an effective process for a single representative (“mean-field”) species. The procedure to derive the effective equations is well documented ; therefore, we only report the final result (a more detailed derivation can be found in the SM). The effective single-species process for the model is given by

̇ 𝑥 (𝑡) = 𝑥 (𝑡) {1 - 𝑢 𝑥 (𝑡) - 𝛼 \int 𝑡 0 𝑑 𝑡' [𝑐 𝐺 (𝕀 - 𝐺) - 1 + Γ (1 - 𝑐) 𝐺] (𝑡, 𝑡') 𝑥 (𝑡') - 𝜂 (𝑡)},

(7)

where $𝕀$ is the identity operator and $𝜂 (𝑡)$ is colored Gaussian noise with zero mean and correlations in time, given by

⟨ 𝜂 (𝑡) 𝜂 (𝑡') ⟩ = 𝛼 [𝑐 (𝕀 - 𝐺) - 1 𝐶 (𝕀 - 𝐺 𝑇) - 1 + (1 - 𝑐) 𝐶] (𝑡, 𝑡') .

(8)

The order parameters in Eqs. are to be obtained self-consistently from the following expressions:

𝑀 (𝑡) = ⟨ 𝑥 (𝑡) ⟩ *, 𝐶 (𝑡, 𝑡') = ⟨ 𝑥 (𝑡) 𝑥 (𝑡') ⟩ *, 𝐺 (𝑡, 𝑡') = 𝛿 𝛿 𝜂 (𝑡') ⟨ 𝑥 (𝑡) ⟩ *,

(9)

where the average $⟨ \dots ⟩ *$ is performed over realizations of the process in Eq. . Equations form a closed system and have to be solved self-consistently.

B. Fixed-point analysis

There is no realistic prospect for a general analytical solution of the effective dynamics in Eq. . One alternative is to use Monte Carlo methods to construct sample paths for the effective process and solutions for the dynamic order parameters, for example via the Eissfeller-Opper procedure or using the more recent approach in Ref. . The latter reference explicitly discusses applications to random Lotka-Volterra systems with Gaussian disorder.

Here we will instead follow Refs. and focus on analytical solutions in the parameter regime in which the dynamics approach stable fixed points. This is motivated by observations from the numerical integration of Eq. . We find that, for certain parameters, the system tends to a unique fixed point, which is independent of initial conditions. Figure shows examples of parameter regions in which this is the case. Broadly speaking, we observe two different types of behavior: (i) the population densities converge to a fixed point or (ii) they diverge. These types of behavior occur in different regions of parameter space (Fig. ). There is a thin boundary between the two regions where other behavior (e.g., periodic behavior or persistent irregular motion) can appear, as evidenced by the occasional green or light blue pixel in Fig. . We attribute this to the fact that the system size $𝑁$ is necessarily finite in numerical experiments, and we expect that this behavior will become increasingly rare as $𝑁 \to \infty$ .

FIG. 1.

Fixed points and divergences of the dynamics in Eq. . Each panel illustrates the behavior of the model for different choices of $𝑐$ and $Γ$ . The heatmap indicates the fraction of samples that converge to a fixed point after numerical integration of the gLVE. The criteria for the identification of convergence or divergence are described in Appendix . The dashed lines are theoretical predictions for the onset of divergence (see Sec. )

We will thus assume that each path in the ensemble of trajectories of the effective process eventually arrives at a unique fixed point, $𝑥 = l i m 𝑡 \to \infty 𝑥 (𝑡)$ . Each realization of the noise variable $𝜂 (𝑡)$ in Eq. also approaches a stationary value $𝜂$ . We note that $𝑥$ and $𝜂$ will be random variables, differing across realizations of the effective dynamics. We can then write

𝑀 = ⟨ 𝑥 ⟩ *, 𝐺 (𝜏) = l i m 𝑡 \to \infty 𝐺 (𝑡 + 𝜏, 𝑡), 𝑞 = l i m 𝑡 \to \infty 𝐶 (𝑡 + 𝜏, 𝑡) = ⟨ 𝑥 2 ⟩ * .

(10)

These relations can be understood as follows: If all realizations of the effective dynamics approach stationary values, then $𝑀 (𝑡)$ will approach a constant, given by $⟨ 𝑥 ⟩ *$ . Furthermore, we assume that the response function $𝐺 (𝑡, 𝑡')$ becomes time-translation invariant for large $𝑡$ , i.e., $𝐺 (𝑡, 𝑡') = 𝐺 (𝑡 - 𝑡')$ . Causality implies that $𝐺 (𝑡 - 𝑡') = 0$ for $𝑡 < 𝑡'$ . Finally, given that all trajectories of the effective dynamics approach fixed points, the correlation function $𝐶$ loses all time dependence and so we have written $𝐶 (𝑡, 𝑡') = 𝑞$ . This is consistent with Eq. ; the noise variables $𝜂 (𝑡)$ also approach a random but time-independent value for all realizations. The mean of the random variable $𝜂$ is zero and using Eq. , its variance is given by

⟨ 𝜂 2 ⟩ = 𝛼 𝑞 [𝑐 (1 - 𝜒) 2 + 1 - 𝑐],

(11)

where

𝜒 = \int \infty 0 𝑑 𝜏 𝐺 (𝜏) .

(12)

From now on, we will write $𝜂 = \sqrt 𝑞 Σ 𝑧$ , where $𝑧$ is a standard Gaussian random variable, and

Σ 2 = 𝛼 [𝑐 (1 - 𝜒) 2 + 1 - 𝑐] .

(13)

Setting the time derivative on the left-hand side of Eq. to zero and using Eqs. we find

𝑥 {1 - 𝑢 𝑥 - \sqrt 𝑞 Σ 𝑧 - 𝛼 𝑥 [𝑐 𝜒 1 - 𝜒 + Γ (1 - 𝑐) 𝜒]} = 0 .

(14)

For a given value of $𝑧$ this is to be solved for $𝑥$ , subject to the constraint that abundances are non-negative, i.e., $𝑥 (𝑧) \geq 0$ . Irrespective of the value of $𝑧$ , Eq. always has the solution $𝑥 (𝑧) = 0$ . Additionally, a second non-negative solution is possible for some values of $𝑧$ . As we will confirm in simulations, the physically meaningful solution is given by

𝑥 = max {0, 1 - \sqrt 𝑞 Σ 𝑧 𝑢 + 𝛼 [𝑐 𝜒 1 - 𝜒 + Γ (1 - 𝑐) 𝜒]} .

(15)

For given order parameters $𝑞$ and $𝜒$ , the function $𝑥 (𝑧)$ in Eq. is therefore piecewise linear, with one piece equal to zero. The denominator in Eq. always comes out positive. Therefore, we have the solution $𝑥 (𝑧) > 0$ when $𝑧 < Δ$ , and $𝑥 (𝑧) = 0$ when $𝑧 \geq Δ$ , with

Δ = 1 Σ \sqrt 𝑞 .

(16)

Given that $𝑧$ is a Gaussian random variable, the abundances of extant species at the fixed point follow a clipped Gaussian distribution. This is similar to what was reported in other random Lotka-Volterra models; see e.g., Ref. . An explicit example of a species abundance distribution can be found for instance in Ref. .

Using this fixed-point ansatz, the relations for the order parameters in Eq. can be written in the following form :

𝑀 = \int Δ - \infty 𝐷 𝑧 𝑥 (𝑧), 𝜒 = 1 \sqrt 𝑞 Σ \int Δ - \infty 𝐷 𝑧 𝜕 𝑥 (𝑧) 𝜕 𝑧, 𝑞 = \int Δ - \infty 𝐷 𝑧 𝑥 (𝑧) 2,

(17)

where $𝐷 𝑧 = 𝑑 𝑧 \sqrt 2 𝜋 𝑒 - 𝑧 2 2$ . It is now convenient to introduce the following functions:

𝑓 𝑛 (Δ) : = \int Δ - \infty 𝐷 𝑧 (Δ - 𝑧) 𝑛,

(18)

for $𝑛 = 0, 1, 2$ . We then find from Eqs.

- 𝜒 {𝑢 + 𝛼 [𝑐 𝜒 1 - 𝜒 + Γ (1 - 𝑐) 𝜒]} = 𝑓 0 (Δ), 𝑀 𝑢 + 𝛼 [𝑐 𝜒 1 - 𝜒 + Γ (1 - 𝑐) 𝜒] \sqrt 𝛼 𝑞 [𝑐 (1 - 𝜒) 2 + (1 - 𝑐)] = 𝑓 1 (Δ), {𝑢 + 𝛼 [𝑐 𝜒 1 - 𝜒 + Γ (1 - 𝑐) 𝜒]} 2 𝛼 [𝑐 (1 - 𝜒) 2 + (1 - 𝑐)] = 𝑓 2 (Δ) .

(19)

Equations together with Eq. form a closed set for the set of unknowns $𝑞$ , $𝜒$ , $𝑀$ , and $Δ$ , which is to be solved as a function of the model parameters $𝑢$ , $𝑐$ , $Γ$ , and $𝛼$ .

Recalling that $𝑥 (𝑧) > 0$ if, and only if, $𝑧 < Δ$ we identify $𝑓 0 (Δ)$ as the fraction of surviving species,

𝜙 \equiv 𝑓 0 (Δ) = \int Δ - \infty 𝐷 𝑧 .

(20)

Equations can be solved parametrically. We fix $𝑢$ , $𝑐$ , $Γ$ , and $Δ$ and then solve for the set of $𝜒$ , $𝑞$ , $𝑀$ , and $𝛼$ .

In detail, we find the following cubic equation for $𝜒$ , valid for $𝑐 < 1$ ,

0 = 𝑓 0 (𝑐 (Γ - 1) 𝑓 0 - Γ 𝑓 0 - 𝑓 2) + 𝜒 (𝑓 20 [𝑐 + 2 Γ - 2 𝑐 Γ] + 2 𝑓 0 𝑓 2 [1 - 𝑐] - 𝑓 2 𝑢) + 𝜒 2 (𝑐 - 1) (Γ 𝑓 20 + 𝑓 0 𝑓 2 - 2 𝑓 2 𝑢) + 𝜒 3 𝑢 𝑓 2 (𝑐 - 1) .

(21)

Further, we have from Eqs. ,

𝑀 = 𝜒 𝑓 21 𝑓 0 (𝑓 0 - 𝑓 2), 𝑞 = 𝜒 2 (𝑓 1 𝑓 2 - 𝑓 0) 2 𝑓 2 𝑓 20, 𝛼 = 𝑓 20 𝑓 2 1 𝜒 2 [1 - 𝑐 + 𝑐 (1 - 𝜒) 2],

(22)

where the $𝑓 𝑛$ are to be evaluated at $Δ$ .

The relations in Eqs. and are also valid for $𝑐 = 1$ and can then be simplified as outlined in Appendix .

The validity of the predictions from Eqs. and is confirmed by direct numerical integration of the gLVE in Fig. .

FIG. 2.

Test of analytical predictions for the order parameters against numerical simulations. The figure shows the fraction of surviving species $𝜙$ , the mean abundance $𝑀$ , and the variance of abundances ( $𝑞 - 𝑀 2$ ) as a function of the model parameter $𝛼$ (where $𝑃 = 𝛼 𝑐 𝑁$ is the number of traits each species in the original system possesses). Lines are from the theory, derived in Eqs. and , markers from numerical integration of the gLVE ( $𝑁 = 1000$ , $𝑡 m a x = 30$ , averaged over 10 realizations of the disorder). The remaining model parameters are $𝑐 = 0.5$ and $Γ = 0.3$ . Vertical dashed lines indicate the onset of divergence as determined from the theory in Sec. .

C. Stability analysis

1. Diverging abundance

a. Model with $𝑐 < 1$

The first and second relations in Eqs. indicate that the order parameters $𝑀$ and $𝑞$ both diverge in the system with $𝑐 < 1$ when $𝑓 0 (Δ) = 𝑓 2 (Δ)$ . The latter implies $Δ = 0$ . The value of $𝛼$ for which this occurs can (for a given choice of $𝑐$ , $𝑢$ , and $Γ$ ) be obtained from the third relation in Eqs. , with $𝜒$ being the relevant root of Eq. . Using Eq. the susceptibility $𝜒$ is found to remain finite at the transition. We note that $𝑓 0 (Δ) > 0$ for all relevant values of $Δ$ .

b. Model with $𝑐 = 1$

The fully connected system also shows two types of divergences: (i) The quantities $𝑀$ and $𝑞$ both diverge when $𝑓 0 (Δ) = 𝑓 2 (Δ)$ ; see Eqs. . The susceptibility then remains finite; (ii) Eqs. further indicate that $𝑀$ and $𝑞$ also diverge in the model with $𝑐 = 1$ when $𝑓 20 (Δ) = 𝑢 𝑓 2 (Δ)$ . This latter condition results in $𝛼 = 𝑢$ . From Eqs. the susceptibility $𝜒$ is then seen to diverge as well (the divergences of $𝑀$ , $𝑞$ , and $𝜒$ take place simultaneously).

We note that the divergencies resulting from $𝑓 0 = 𝑓 2$ and $𝑓 20 = 𝑢 𝑓 2$ can take place at different locations in parameter space for the model with $𝑐 = 1$ . If this is the case, and starting in the stable phase, then the divergence that occurs first will determine the loss of stability in the fully connected system. For $𝑢 < 1 / 2$ the transition of type (ii) takes place first as $𝛼$ is increased ( $𝑀$ , $𝑞$ , and $𝜒$ diverge), and for $𝑢 > 1 / 2$ the transition of type (i) is instead observed ( $𝑞$ and $𝑀$ diverge and $𝜒$ remains finite). At present we do not have any further intuition regarding any significance or special role of the value $𝑢 = 1 / 2$ .

2. Linear instability

The system also shows a linear instability which can be identified using the procedure established in Refs. . We write $𝑥 (𝑡) = 𝑥 + 𝑦 (𝑡)$ and $𝜂 (𝑡) = \sqrt 𝑞 Σ 𝑧 + 𝜁 (𝑡)$ , where $𝑦 (𝑡)$ and $𝜁 (𝑡)$ are small perturbations about the fixed point of the trajectories of the effective process in Eq. . Expanding to first order in these perturbations we find that

̇ 𝑦 (𝑡) = 𝑥 [- 𝑢 𝑦 (𝑡) - 𝛼 \int 𝑡 - \infty 𝑑 𝑡' 𝐾 (𝑡, 𝑡') 𝑦 (𝑡') - 𝜁 (𝑡)] + 𝑦 (𝑡) [1 - 𝑢 𝑥 - 𝛼 𝑥 \int 𝑡 - \infty 𝑑 𝑡' 𝐾 (𝑡, 𝑡') - \sqrt 𝑞 Σ 𝑧],

(23)

with $𝐾 (𝑡, 𝑡') = [𝑐 𝐺 (𝕀 - 𝐺) - 1 + Γ (1 - 𝑐) 𝐺] (𝑡, 𝑡')$ . We also have the self-consistency relation

⟨ 𝜁 (𝑡) 𝜁 (𝑡') ⟩ = 𝛼 [𝑐 (𝕀 - 𝐺) - 1 𝐷 (𝕀 - 𝐺 𝑇) - 1 + (1 - 𝑐) 𝐷] (𝑡, 𝑡'),

(24)

where $𝐷 (𝑡, 𝑡') = ⟨ 𝑦 (𝑡) 𝑦 (𝑡') ⟩ *$ .

When $𝑥 = 0$ , Eq. becomes

̇ 𝑦 (𝑡) = 𝑦 (𝑡) (1 - \sqrt 𝑞 Σ 𝑧) .

(25)

Equation , together with the observation that the denominator in this equation is strictly positive, implies that $1 - \sqrt 𝑞 Σ 𝑧 < 0$ when $𝑥 = 0$ . This allows us to conclude that perturbations on extinct species decay and do not contribute to any linear instability.

For fixed points $𝑥 > 0$ we find from Eqs. and that

̇ 𝑦 (𝑡) = - 𝑥 [𝑢 𝑦 (𝑡) + 𝛼 \int 𝑡 - \infty 𝑑 𝑡' 𝐾 (𝑡 - 𝑡') 𝑦 (𝑡') + 𝜁 (𝑡)] .

(26)

To identify the onset of linear instability we follow Refs. . We move to Fourier space, writing $𝜔$ for the variable conjugate to time $𝑡$ and using tildes to indicate Fourier transforms.

Focusing on the mode with $𝜔 = 0$ and following steps similar to those in Refs. we then find from Eq.

⟨ | ̃ 𝑦 (0) | 2 ⟩ = {𝜙 - 1 [𝑢 + 𝛼 𝑐 𝜒 1 - 𝜒 + 𝛼 Γ (1 - 𝑐) 𝜒] 2 - 𝛼 [𝑐 (1 - 𝜒) 2 + 1 - 𝑐]} - 1 .

(27)

The left-hand side is manifestly non-negative, so a change of sign of the expression inside the square bracket on the right-hand side indicates an inconsistency (and divergence of $⟨ | ̃ 𝑦 (0) | 2 ⟩$ ). Using Eqs. this is shown to occur when

𝛼 [𝑐 (1 - 𝜒) 2 + 1 - 𝑐] (𝑓 0 - 𝑓 2) = 0 .

(28)

For $𝑐 < 1$ the expression in the square brackets is never zero. This leaves us with the condition $𝑓 0 = 𝑓 2$ , which is the same as we obtained for the divergence of $𝑀$ and $𝑞$ . If $𝑐 = 1$ , then the term in the square bracket is zero if $𝜒 \to \infty$ , which using Eq. we can write as $𝑓 20 - 𝑢 𝑓 2 = 0$ .

From this, we conclude that in our model the linear instability is always accompanied by the instability with diverging mean abundance. This is markedly different from the behavior of the gLVE model with Gaussian random interactions. In this Gaussian model there are instances where the linear instability sets in as the variance of interactions is increased, but where abundances remain finite and the divergence only occurs at a later point at even higher variance of the interactions. This leads to a phase with multiple attractors between the two transitions . Our analysis indicates that the model with Hopfield-like couplings does not have such a multiple-attractor regime.

IV. PHASE DIAGRAM AND FURTHER BEHAVIOR OF THE MODEL

A. Phase diagram for the fully connected system ( $𝑐 = 1$ )

The phase diagram of the fully connected model is shown in Fig. . We recall that, for $𝑐 = 1$ , the only model parameters are the self-interaction coefficient $𝑢$ and the ratio of the number of traits to the number of interspecies interactions in the original pool ( $𝛼 = 𝑃 / 𝑐 𝑁$ ). For a fixed value of $𝛼$ , the system shows a unique stable fixed point for $𝑢 > 𝑢 𝑐 (𝛼)$ , where $𝑢 𝑐 (𝛼)$ marks the onset of instability. The line in Fig. , obtained from Eq. , shows the phase boundary between the stable and unstable regions. At this boundary $𝑀$ and $𝑞$ diverge, and if $𝑢 𝑐 < 1 / 2$ we also observe a divergence of $𝜒$ .

FIG. 3.

Phase behavior of the fully connected model ( $𝑐 = 1$ ). (a) Phase diagram for the model with $𝑐 = 1$ , the only model parameters are then $𝑢$ and $𝛼$ . The system is stable to the right of the lines. At the dot-dashed line ( $𝑢 < 1 / 2$ ) $𝑞$ , $𝑀$ , and $𝜒$ all diverge, and at the dashed line $𝑀$ and $𝑞$ diverge, but $𝜒$ remains finite. (b) Illustration of the behavior of the abundances of individual species in the two different phases (convergence to a fixed point shown in green and red, diverging abundances in orange and blue).

The two types of trajectory in the stable and divergent phases are illustrated in Fig. . In the stable phase the system reaches a fixed point for any one realization of the interaction matrix (two examples are illustrated in green and red, respectively).

Figure also shows two examples in which the species abundances diverge (blue and orange). The divergence occurs at a finite time. We will discuss this further in Sec. .

B. Phase diagram for connectivity $𝑐 < 1$

Figure shows how the phase diagram for the system with $𝑐 < 1$ depends on the connectivity $𝑐$ and the symmetry parameter $Γ$ . In all cases there is a single phase boundary, where the divergence of $𝑀$ and $𝑞$ and the onset of linear instabilities coincide. This phase boundary separates a region where trajectories converge to a single globally stable fixed point (phase to the right of the line), from a region where trajectories are unbounded and diverge in finite time (phase to the left).

FIG. 4.

Phase diagram for different choices of the connectivity $𝑐$ , and the symmetry parameter $Γ$ . The colored lines in each panel indicate where the linear instability occurs. The instability coincides with the divergence of $𝑀$ and $𝑞$ . The system is stable to the right of the line, abundances diverge on the left.

The phase diagrams in Figs. and show that the system is in the stable phase for small values of $𝛼$ (i.e., a small number of traits relative to the number of species in the initial pool) or large values of $𝑢$ (i.e., large negative self-interaction). This is the consequence of two competing effects, the self-interaction (parametrized by $𝑢$ ) which stabilizes the system and the interaction between species (induced by competition of similar species) which promotes instability. When $𝑢$ is large and/or $𝛼$ is small, the stabilizing effect of the intra-species interaction dominates over the interactions across species. In the extreme limit $𝛼 = 0$ (no interaction between different species), each abundance follows a separate logistic equation, $̇ 𝑥 𝑖 = 𝑥 𝑖 (1 - 𝑢 𝑥 𝑖)$ , and converges to $𝑥 𝑖 = 1 / 𝑢$ . When $𝛼$ is small but nonzero, the system consists of weakly interacting species. The effect of the interactions between species is then a small perturbation to the logistic behavior of individual species and does not change the convergence to a fixed point. This can be confirmed from Eqs. and by taking the limit $𝛼 \to 0$ , which results in all species surviving with fixed-point abundance $𝑥 * 𝑖 = 1 / 𝑢$ ( $𝜙 \to 1$ , $𝑀 \to 1 𝑢$ , and $Var [𝑥] \to 0$ ). A similar result is obtained for $𝑢 ≫ 1$ at fixed value of $𝛼$ .

Conversely, for low values of $𝑢$ or large values of $𝛼$ the system is unstable. In this situation, the stabilizing self-interaction is not sufficient to overcome the destabilizing effect of the random interactions between species.

The most interesting behavior takes place at the phase boundary, where the effect of the intraspecific and interspecific nonlinearities are of comparable magnitude. From Eqs. and we can conclude that $𝜙 \to 1 / 2$ , and $𝑀 \sim (𝛼 𝑐 - 𝛼) - 1$ as the system approaches the instability (from the stable phase). Further details can be found in Appendix .

We further note that decreasing the value of the symmetry parameter $Γ$ increases the range of the stable region in the phase diagrams in Fig. . This is similar to the effect of increasing the fraction of predator-prey interactions in Lotka-Volterra models with Gaussian interactions . Indeed, the effect of a reduction of $Γ$ is to increase the fraction of species pairs $𝑖, 𝑗$ with $𝑐 𝑖 𝑗 = 1$ and $𝑐 𝑗 𝑖 = 0$ , that is the proportion of unidirectional interactions.

Interestingly, the effect of varying the “connectance” $𝑐$ is not straightforward. As can be seen in Fig. an increased connectivity can, depending on the other model parameters, turn a previously stable system into an unstable one or, vice versa, stabilize a previously unstable system.

C. Finite-time divergence of the mean abundance

As mentioned earlier, the divergence of the abundances in the divergent phase occurs at finite time. This has previously been reported in the model with Gaussian interactions and can be justified heuristically from the Lotka-Volterra equations. Indeed, Eq. has a second-order nonlinearity in the abundances $𝑥 𝑖$ . This can lead to dynamics of the form $̇ 𝑥 \sim 𝑥 2$ , which in turn implies a solution of the form $𝑥 (𝑡) = (𝑐 - 𝑡) - 1$ , where $𝑐$ is an integration constant. This results in a divergence at finite time.

Figure shows the time, $𝑡 d i v$ , at which the divergence occurs for different choices of the model parameters. This time grows as one approaches the stability line (from inside the unstable phase). When the stability line is crossed (into the stable phase), the time-to-divergence diverges itself ( $𝑡 d i v \to \infty$ ), i.e., the divergence no longer occurs. Results from the numerical integration of the gLVE suggest that the divergence of the abundances is of the form $𝑀 \sim (𝑡 d i v - 𝑡) 𝜈$ , where $𝜈 \approx 1$ , as shown in Fig. . This behavior appears to be independent of initial conditions, the values of the parameters $𝑢$ and $𝛼$ , and the initial number of species $𝑁$ .

FIG. 5.

Finite-time divergence of abundances. The heatmaps indicate the time, $𝑡 d i v$ , at which abundances diverge for initial conditions $𝑥 𝑖 (0) = 1$ . Data are obtained from numerical integration of the gLVE. The dotted line is the phase boundary predicted by the theory. To the right of the phase boundary the system is in the stable phase, so that no divergence occurs.

FIG. 6.

Divergence of $𝑀$ for initial conditions uniformly distributed in $(0, 10 - 3)$ and parameters $𝑁 = 1000$ , $𝑢 = 1$ , and different values of $𝛼$ . The dotted black line corresponds to $(𝑡 \infty - 𝑡) - 1$ . The deviation from $𝑀 \sim (𝑡 \infty - 𝑡) - 1$ close to the divergence is attributed to numerical error.

V. REDUCED INTERACTION MATRIX AND ITS EIGENVALUE SPECTRUM

Reference recently established a close connection between different instabilities in the Gaussian random Lotka-Volterra model and the eigenvalue spectrum of the interaction matrix of the surviving species. More specifically, the spectrum of this reduced interaction matrix is composed of a bulk region and a potential outlier eigenvalue. As parameters are changed (starting from within the stable phase), either the bulk spectrum or the outlier eigenvalue can cross into the right half of the complex plane. In the Gaussian model, the crossing of the outlier is associated with a transition marked by the divergence of abundances and a crossing of the bulk with a linear instability.

In this section, we explore in numerical simulations how the different transitions in the gLVE model with Hopfield-like interactions relate to the eigenvalue spectrum of the matrix of interactions between surviving species.

A. Spectrum of the original interaction matrix

Before we discuss the spectra of the reduced interaction matrix, we make a few remarks on the initial interaction matrix $𝛼 𝑖 𝑗 = 𝑐 𝑖 𝑗 𝐽 𝑖 𝑗$ among all species. Throughout this section we set the diagonal elements of this matrix to zero, the only effect of self-interaction (the term $- 𝑢 𝑥 𝑖$ ) is a simple shift of this spectrum. In the large- $𝑁$ limit the central limit theorem applies to $𝐽 𝑖 𝑗 = - 1 𝑐 𝑁 \sum 𝛼 𝑐 𝑁 𝜇 = 1 𝜉 𝜇 𝑖 𝜉 𝜇 𝑗$ , so each off-diagonal entry $𝛼 𝑖 𝑗$ of the interaction matrix is either a Gaussian random variable (if $𝑐 𝑖 𝑗 = 1$ ) or equals zero (if $𝑐 𝑖 𝑗 = 0$ ). The variance of $𝛼 𝑖 𝑗$ is

Var (𝛼 𝑖 𝑗) = 1 𝑐 𝑁 2 𝛼 𝑐 𝑁 \sum 𝜇, 𝜇' ⟨ 𝜉 𝜇 𝑖 𝜉 𝜇' 𝑖 𝜉 𝜇 𝑗 𝜉 𝜇' 𝑗 ⟩ = 𝛼 𝑁 .

(29)

Calculating the correlations between pairs of elements we obtain

Corr [𝛼 𝑖 𝑗, 𝛼 𝑛 𝑚] = ⟨ 𝛼 𝑖 𝑗 𝛼 𝑛 𝑚 ⟩ - ⟨ 𝛼 𝑖 𝑗 ⟩ ⟨ 𝛼 𝑛 𝑚 ⟩ Var (𝛼 𝑖 𝑗) = {Γ (1 - 𝑐) + 𝑐 (𝑖, 𝑗) = (𝑚, 𝑛) 1 (𝑖, 𝑗) = (𝑛, 𝑚) 0 else,

(30)

where we have used Eq. and the fact that $𝐽 𝑖 𝑗$ is symmetric. This means that only diagonally opposed pairs of elements are correlated and that their correlation is determined by both $Γ$ and $𝑐$ .

Based on a theory that only takes into account correlations between diagonally opposed matrix entries, one might then expect an elliptic spectrum , with support given by the ellipse

[𝑥 \sqrt 𝛼 (1 + 𝜏)] 2 + [𝑦 \sqrt 𝛼 (1 - 𝜏)] 2 = 1,

(31)

with $𝜏 = Γ (1 - 𝑐) + 𝑐$ . However, as illustrated in Fig. , this is an approximation to the true spectrum at best for large values of $𝛼$ . For intermediate values of $𝛼$ (an example is shown in orange in the figure), the eigenvalue spectrum appears to have a triangular shape, and for small values of $𝛼$ (shown in green), the spectrum becomes even more skewed and eventually appears to consist of two separate components (example shown in red). While we cannot fully exclude finite-size effects (the spectra in Fig. are for $𝑁 = 5000$ ), we believe that the deviations from an elliptical spectrum in Eq. are due to higher-order correlations between entries of the interaction matrix. For example, it has been shown in Ref. that cyclic correlations can result in eigenvalue spectra with shapes similar to the ones in Fig. . For $𝑐 = 1$ the interaction matrix is a (scaled and shifted) Wishart matrix, so its spectrum follows the Marcenko-Pastur law.

FIG. 7.

Examples of the eigenvalue spectrum of the original interaction matrix. The dashed lines are the naïve predictions of Eq. . Model parameters are $𝑐 = 0.4$ , $Γ = 𝑐 / (𝑐 - 1)$ , and $𝑁 = 5000$ .

B. Eigenvalues of the reduced interaction matrix

We now conclude the analysis of the model with a numerical study of the spectra of the reduced interaction matrix, that is, the interaction matrix between species that survive at the fixed points of the gLVE.

Figure shows the spectra of this matrix for the case $𝑐 < 1$ and for a choice of $Γ$ less than one. This means that the initial interaction matrix is not symmetric. The reduced matrix is not symmetric either, and as a consequence its eigenvalues will generally be complex. As seen in Fig. , the spectrum is not elliptic, and we have found no evidence of an outlier eigenvalue in this scenario ( $𝑐 < 1$ ). In the figure we have fixed $𝑢$ , and varied $𝛼$ . The data suggest that the phase transition at $𝛼 = 𝛼 𝑐 (𝑢)$ coincides with the point at which the rightmost bulk eigenvalue crosses the imaginary axis into the right half-plane.

FIG. 8.

Eigenvalue spectrum of the reduced interaction matrix for $𝑢 = 5$ , $𝑐 = 0.5$ , and $Γ = 𝑐 / (𝑐 - 1)$ , for different choices of the model parameter $𝛼$ . The vertical dashed lines indicate the real part of the rightmost eigenvalue.

In Fig. we study the fully connected system for two different values of the self-interaction strength $𝑢$ . The original interaction matrix in the fully connected model is symmetric by construction and so is the reduced interaction matrix. As a consequence, all eigenvalues are real.

FIG. 9.

Eigenvalue spectrum of the reduced interaction matrix in the fully connected system. The matrices are symmetric, and their eigenvalues are therefore real valued. Panel (a) is for $𝑢 = 0.7$ ; panel (b) is for $𝑢 = 0.3$ .

Figure focuses on the case $𝑢 > 1 / 2$ . We find no signs of outlier eigenvalues, and again the data indicate that the transition to instability occurs when the leading bulk eigenvalue crosses into the positive half of the real axis.

Figure shows a scenario in which $𝑢 < 1 / 2$ . In contrast with the situation in (a), an outlier eigenvalue now becomes apparent, and the transition to instability in the gLVE at $𝛼 = 𝛼 𝑐 (𝑢)$ now appears to coincide with the point at which the outlier becomes positive.

The connection between the different types of transition and the behavior of the spectrum of the reduced interaction matrix is shown in Table . We recall that the mean abundance and the second moment of the abundances diverge at all transitions and that the onset of the linear instability always coincides with the point of diverging abundances. There are thus only two types of transition, one in which the susceptibility remains finite ( $𝜒 < \infty$ ) and another for which it diverges ( $𝜒 \to \infty$ ). The table indicates that the former transition ( $𝜒$ finite) appears to coincide with the bulk spectrum of the reduced matrix crossing into the right half of the complex plane. The transition at which $𝜒 \to \infty$ (along with the divergences of $𝑀$ and $𝑞$ ), on the other hand, seems to be seen when the outlier eigenvalue of the reduced matrix in the fully connected system reaches the origin.

TABLE I.

Types of phase transition in the gLVE model with Hopfield-like interactions. The table summarizes the different transitions, giving details about the nature of the divergence at the transition, and the associated behavior of the spectrum of the reduced interaction matrix.
		$𝑞, 𝑀$ diverge	bulk spectrum
$𝑐 < 1$		$𝜒$ remains finite	crosses axis
$𝑐 = 1$	$𝑢 > 1 / 2$	$𝑞, 𝑀$ diverge	bulk spectrum
		$𝜒$ remains finite	crosses axis
	$𝑢 < 1 / 2$	$𝑞, 𝑀, 𝜒$	outlier eigenvalue
		all diverge	crosses axis
		(at $𝑢 = 𝛼$ )

We stress that these are numerical observations and that these findings should therefore be seen mostly as conjectures at this point. In principle, the spectrum of the reduced interaction matrix can likely be calculated in our model, adapting the method used in Ref. . However, this involves a substantial calculation and is beyond the scope of the current paper.

VI. DISCUSSION

To summarize, we have carried out a generating functional analysis of a random generalized Lotka-Volterra system with interactions determined by niche overlap. Species interactions in the model are governed by Hopfield-like couplings subject to mild dilution (the remaining connectivity is still extensive). We have computed the statistics of surviving species in the stable fixed-point phase, and we have analytically determined the onset of instability. Similarly to the gLVE with Gaussian interactions, asymmetry in the connectivity matrix promotes stability. That is, the system becomes more stable when there is a larger fraction of unidirectional interactions ( $𝑐 𝑖 𝑗 = 1$ , but $𝑐 𝑗 𝑖 = 0$ ). In contrast with the Gaussian model, the linear instability against small perturbations cannot be separated from an instability at which species abundances diverge. As a consequence, there is no phase with multiple stable fixed points our model. Despite some common features, the statistical mechanics of the Gaussian and Hopfield-like models are therefore rather distinct.

Our analysis shows two types of transitions to divergent abundances, one in which the integrated response $𝜒$ remains finite and another in which $𝜒$ diverges. This raises interesting questions about the exact nature of memory onset in the system (a diverging integrated response indicates persistent memory of perturbations). Future work could focus on the precise shape of the response function, where the numerical methods in Ref. might prove particularly useful. Given that the fully connected system has symmetric couplings it would also be interesting to see how crossing each of the different types of transition affects the energy landscape. A natural approach here might be the replica method and suitable levels of replica symmetry breaking .

Numerical simulations provide evidence that the transition at which the integrated response remains finite ( $𝜒 < \infty$ ) is associated with the bulk spectrum of the reduced interaction matrix (the matrix of interactions between extant species) crossing the axis. The transition at which $𝜒$ diverges on the other hand appears to be signalled by an outlier eigenvalue crossing the imaginary axis.

These findings in simulations reinforce the intriguing analytical result obtained recently in Ref. . Namely, the eigenvalues of the interaction matrix in the community of surviving species can be used to decide the stability of feasible equilibria, that is, fixed points with non-negative species abundances. In the traditional approach to ecosystem stability by Robert May , based on the eigenvalue spectra of random matrices, no actual dynamics are specified, and the feasibility of the assumed equilibria remains unclear. Any fixed point of the generalized Lotka-Volterra model on the contrary is feasible by construction. The study of the spectra of reduced interaction matrices resulting from Lotka-Volterra dynamics can therefore contribute to establishing how May's approach can be adapted to include feasible equilibria. Noting that previous work has shown that the statistics of the reduced interaction matrix in random Lotka-Volterra models can be quite different from those of the original interaction matrix, it would be interesting to study the statistics of the $𝐽 𝑖 𝑗$ among survivors in the present model in more detail. In particular, the niche overlap between surviving species.

On a broader level, our study highlights two common facets of work on the statistical physics of complex systems, which were also seen for example in the early 1980s to 1990s when physicists studied neural networks or the early 2000s to 2010 when a number of physicists worked on the Minority Game. On the one hand, tools from physics can make a difference in problems from other disciplines. In our system (and other models of complex ecosystems more generally) this is the study of feasible equilibria with methods from spin-glass physics. At the same time, studying problems arising in other areas can reveal new types of physics and complexity, which one would perhaps not find within the strict boundaries of traditional physics. In our case, these are the different types of phase transition in the generalized Lotka-Volterra model. We think that this mutually beneficial relation of physics and adjacent disciplines is what makes the field of complex systems particularly attractive.

ACKNOWLEDGMENTS

Partial financial support has been received from the Agencia Estatal de Investigación and Fondo Europeo de Desarrollo Regional (FEDER, UE) under Project APASOS (PID2021-122256NB-C21/PID2021-122256NB-C22) and the Maria de Maeztu program for Units of Excellence, CEX2021-001164-M, funded by MCIN/AEI/10.13039/501100011033.

Appendices

APPENDIX A: DETAILS OF NUMERICAL PROCEDURES

For the numerical integration of the gLVE we use scypi's solve_ivp function, which uses a RK45 integration scheme.

To determine the fraction of survivors we count the number of species above a threshold abundance of $10 - 4$ . There are two sources of systematic error associated with this method. The most relevant is the overestimation of the fraction of survivors if the system is not close enough to the equilibrium configuration. This can be addressed by extending the simulation time.

The second source of error comes from the fact there is no “gap” between zero and the lowest nonzero abundance [see the clipped Gaussian Eq. ]. This implies that for any value of the threshold, there is a nonzero probability of finding equilibrium abundances below it. A possible solution, making use of the facts that in simulations $𝑁$ is finite and that we know the abundance distribution analytically, is to choose the threshold value so that the expected number of surviving species with an abundance below the threshold is small (e.g., smaller than 1). The chosen value of $10 - 4$ provides good results in the parameter ranges we have explored.

As part of our measurements, it is necessary to detect divergences in the species' abundances. To detect this divergence we have used the failure of the integration method as an indicator. Indeed, as the abundances grow with each iteration, so does the estimated error used to adapt the step size. This causes the solver to lower the time step until it eventually drops below machine precision, at which point integration is stopped. The agreement of the theoretical and numerical phase boundaries in Fig. confirms the validity of this method.

APPENDIX B: ORDER PARAMETERS AT THE FIXED POINT OF THE FULLY CONNECTED SYSTEM

The parametric solution for the order parameters of the fully connected system ( $𝑐 = 1$ ) in the fixed-point phase can be obtained from the following relations:

𝛼 = (𝑢 + 𝑓 0) 2 𝑓 2 (𝑓 0 + 𝑓 2) 2, 𝑀 = 𝑓 21 (𝑓 0 + 𝑓 2) (𝑓 0 - 𝑓 2) (𝑓 20 - 𝑢 𝑓 2), 𝑞 = 𝑓 21 𝑓 2 (𝑓 0 + 𝑓 2) 2 (𝑓 0 - 𝑓 2) 2 (𝑓 20 - 𝑢 𝑓 2) 2, 𝜒 = 𝑓 0 (𝑓 0 + 𝑓 2) 𝑓 20 - 𝑢 𝑓 2 .

(B1)

The functions $𝑓 𝑛 (Δ)$ on the right provide $𝛼$ , $𝑀$ , $𝑞$ , $𝜒$ as implicit functions of $Δ$ .

Keeping Eq. in mind one sees that $𝑀$ and $𝑞$ can only diverge if $𝑓 0 = 𝑓 2$ or $𝑓 20 = 𝑢 𝑓 2$ , as indicated in the main text.

APPENDIX C: LIMITING BEHAVIOR OF THE ORDER PARAMETERS

1. Limit $𝛼 \to 0$

The weak interaction limit $𝛼 \to 0$ corresponds to $Δ \to \infty$ . [This can be seen from Eq. , keeping in mind that $𝑓 0 > 0$ .] From the definition in Eqs. we have, in this limit,

𝑓 0 (Δ) = 𝜙 = 1 + 𝑒 - Δ 2 2 \sqrt 2 𝜋 [- 1 Δ + 𝒪 (Δ - 3)],

(C1)

𝑓 1 (Δ) = Δ + 𝑒 - Δ 2 2 \sqrt 2 𝜋 [1 Δ 2 + 𝒪 (Δ - 4)],

(C2)

𝑓 2 (Δ) = 1 + Δ 2 + 𝑒 - Δ 2 2 \sqrt 2 𝜋 [- 2 Δ 3 + 𝒪 (Δ - 5)] .

(C3)

Next we compute the value of $𝜒$ . Since only $𝑓 0$ and $𝑓 2$ are present in Eq. , and only $𝑓 2$ is divergent, we group in terms proportional to $𝑓 2$ to obtain

0 = 1 - 𝜒 [2 (1 - 𝑐) - 𝑢] - 𝜒 2 (𝑐 - 1) (1 - 2 𝑢) - 𝜒 3 𝑢 (𝑐 - 1),

which we can check always has $𝜒 = - 1 / 𝑢$ as its negative solution. Finally, from Eq. we obtain $𝑀 = 1 / 𝑢$ and $Var [𝑥] = 0$ .

As expected, these values are independent of $𝑐$ and $Γ$ , since in the limit of absent interactions Eq. becomes a set of independent logistic maps. In this case, all species survive with abundance $1 / 𝑢$ , which is what we obtain.

2. Limit $𝛼 \to 𝛼 𝑐$

There are two different scenarios for the limit $𝛼 \to 𝛼 𝑐$ (where $𝛼 𝑐$ is the location of the phase transition).

(1) If $𝑐 = 1$ , $𝑢 < 1 / 2$ , then the divergence takes place as $𝑢 𝑓 2 \to 𝑓 20$ . Using Eq. we see that $𝛼 𝑐 = 𝑢$ , and

𝛼 𝑐 - 𝛼 = 𝑢 - 𝑓 2 (𝑓 0 + 𝑓 2) 2 (𝑓 20 - 𝑢 𝑓 2) .

(C4)

This implies that both $𝜒$ and $𝑀$ diverge as $(𝛼 𝑐 - 𝛼) - 1$ and $𝑞$ diverges as $(𝛼 𝑐 - 𝛼) - 2$ .

(2) The other type of transition occurs when $𝑓 0 \to 𝑓 2$ , which implies $Δ \to 0$ . In this case $𝜒$ remains finite and we have near $Δ = 0$ ,

𝑓 0 (Δ) = 12 + Δ \sqrt 2 𝜋 - Δ 3 6 \sqrt 2 𝜋 + 𝒪 (Δ 4) 𝑓 1 (Δ) = 1 \sqrt 2 𝜋 + Δ 2 + Δ 2 2 \sqrt 2 𝜋 + 𝒪 (Δ 4) 𝑓 2 (Δ) = 12 + \sqrt 2 𝜋 Δ + Δ 2 2 + Δ 3 3 \sqrt 2 𝜋 + 𝒪 (Δ 4) .

(C5)

We note from these expansions that $𝑓 0 - 𝑓 2 \propto Δ$ as $Δ \to 0$ . Using Eq. we then find $𝑀 \sim Δ - 1$ . Similarly, we have $𝛼 𝑐 - 𝛼 \propto Δ$ [this can be seen from expanding $𝑓 20 / 𝑓 2$ in the third relation in Eq. ], so we can conclude that $𝑀 \sim (𝛼 - 𝛼 𝑐) - 1$ . These results are consistent with simulations (see for example Fig. ).

Niche overlap and Hopfield-like interactions in generalized random Lotka-Volterra systems

Abstract

Physics Subject Headings (PhySH)

Article Text

I. INTRODUCTION

II. MODEL DEFINITIONS

III. GENERATING FUNCTIONAL ANALYSIS AND STABILITY

A. Generating functional and effective process

B. Fixed-point analysis

C. Stability analysis

1. Diverging abundance

a. Model with $𝑐 < 1$

b. Model with $𝑐 = 1$

2. Linear instability

IV. PHASE DIAGRAM AND FURTHER BEHAVIOR OF THE MODEL

A. Phase diagram for the fully connected system ( $𝑐 = 1$ )

B. Phase diagram for connectivity $𝑐 < 1$

C. Finite-time divergence of the mean abundance

V. REDUCED INTERACTION MATRIX AND ITS EIGENVALUE SPECTRUM

A. Spectrum of the original interaction matrix

B. Eigenvalues of the reduced interaction matrix

VI. DISCUSSION

ACKNOWLEDGMENTS

Appendices

APPENDIX A: DETAILS OF NUMERICAL PROCEDURES

APPENDIX B: ORDER PARAMETERS AT THE FIXED POINT OF THE FULLY CONNECTED SYSTEM

APPENDIX C: LIMITING BEHAVIOR OF THE ORDER PARAMETERS

1. Limit $𝛼 \to 0$

2. Limit $𝛼 \to 𝛼 𝑐$

Supplemental Material

References (40)

Niche overlap and Hopfield-like interactions in generalized random Lotka-Volterra systems

Abstract

Physics Subject Headings (PhySH)

Article Text

I. INTRODUCTION

II. MODEL DEFINITIONS

III. GENERATING FUNCTIONAL ANALYSIS AND STABILITY

A. Generating functional and effective process

B. Fixed-point analysis

C. Stability analysis

1. Diverging abundance

a. Model with 𝑐<1

b. Model with 𝑐=1

2. Linear instability

IV. PHASE DIAGRAM AND FURTHER BEHAVIOR OF THE MODEL

A. Phase diagram for the fully connected system (𝑐=1)

B. Phase diagram for connectivity 𝑐<1

C. Finite-time divergence of the mean abundance

V. REDUCED INTERACTION MATRIX AND ITS EIGENVALUE SPECTRUM

A. Spectrum of the original interaction matrix

B. Eigenvalues of the reduced interaction matrix

VI. DISCUSSION

ACKNOWLEDGMENTS

Appendices

APPENDIX A: DETAILS OF NUMERICAL PROCEDURES

APPENDIX B: ORDER PARAMETERS AT THE FIXED POINT OF THE FULLY CONNECTED SYSTEM

APPENDIX C: LIMITING BEHAVIOR OF THE ORDER PARAMETERS

1. Limit 𝛼→0

2. Limit 𝛼→𝛼𝑐

Supplemental Material

References (40)

a. Model with $𝑐 < 1$

b. Model with $𝑐 = 1$

A. Phase diagram for the fully connected system ( $𝑐 = 1$ )

B. Phase diagram for connectivity $𝑐 < 1$

1. Limit $𝛼 \to 0$

2. Limit $𝛼 \to 𝛼 𝑐$