Marginally stable equilibria in critical ecosystems

Ecosystems can be incredibly rich and their understanding remains an outstanding challenge. A number of aspects of the problem make the tools of statistical mechanics uniquely suitable for the challenge: highly-diverse ecological communities involve many degrees of freedom (e.g., the species), with intricate interactions. At the same time, there is special interest in system-wide, aggregate quantities, such as the total biomass, overall diversity, system stability and fluctuations, and the response to large-scale perturbations.

The use of statistical mechanics in a growing body of works has provided insight into the behavior of diverse ecosystems. In particular, they have highlighted role of collective phenomena, and the emergence of phases with distinct qualitative behaviors, with phase-transitions between them [1–10].

Here we study a canonical model for a community of interacting species. We identify a phase which is critical, in the language of condensed matter physics. This means that all stable equilibria are in fact very close to marginal stability. Beyond the implications to ecology, discussed throughout this paper, this observation makes contact with many complex systems in Nature that are poised just at the edge of stability [11–19]. One important common trait of all examples is that they are formed by strongly interacting units—species, neurons, agents and particles depending on the situation. The possible explanations of such phenomenon are varied. They include the need for flexibility and adaptiveness to time-varying conditions [12, 20], balance between functionality and stability [20], self-organized criticality [21], self-organized instability [22], and continuous constraints satisfaction [11]. In this work we provide exact results which uncover an underlying mechanism in a many-variable, interacting system.

We focus on the generalized Lotka–Volterra (LV) equations. They provide a simple and general setting to study assemblies of interacting degrees of freedom; as such they are used in several fields [23–26]. In particular, they provide a canonical model for ecosystems, with growing connections to systems across biology [23, 24, 27]. The study of stability of equilibria and their properties using LV equations and generalizations has become a very active research subject. Several important results were obtained recently; in particular general techniques to count the number of equilibria and their properties have been developed [8], and criticality and glassiness have been found to be emergent properties of ecosystems [2–7]. A phase-transition similar to the one we find was identified in the related replicator equations [2, 3], but the properties of the new phase were not analyzed. More recently, a phase-diagram for the LV equations was suggested in [7], see the Discussion section for the relation to the present work.

By mapping the model to condensed matter system, our approach reveal the generality of the results beyond LV models. Moreover, by transferring knowledge developed in glass physics, we propose ways to directly check the degree of criticality using experimental data on abundance distributions, see the Discussion section below.

In the model we consider, an ecological community is assembled from a pool of available species. We focus on the cases where the number of species is large. Since the detailed parameters of all interactions are not known in the majority of cases, and in any case not all details are expected to matter [28], we follow the long tradition pioneered by May in ecology [1] and Wigner in physics [29], and sample the interactions randomly. As recently shown in [28], this framework encompasses a wide array of ecological models, from resource competition to predation and mutualism. It is important to remark that we go beyond May’s classical work since randomness is here introduced at the level of interactions between all possible species, while the community self-organizes by choosing which species are present. In other words, the number and identity of the species that are present in the community is selected dynamically [2, 30]. Understanding the emergent stability of the equilibria reached dynamically and its dependence on the external parameters is the main purpose of this work.

We find, in agreement with [5, 6, 9], that when the interactions are weak or highly uniform, only one equilibrium is present and is determined mainly by self-regulation within each species. For stronger and more heterogeneous interactions, multiple equilibria emerge. Our main result is that when this happens, all possible states of the system are close to be marginally stable for large number of species and this determines the diversity of the ecosystem, see figure 1. Marginal stability has several important consequences, in particular it leads to extreme susceptibility to small perturbations. This situation is referred as ‘critical’ in the physics literature [31]. In a well-known work, May suggested that complexity and interactions limit the stability of ecosystems [1]. Our results provide a complementary perspective: complex ecological communities reduce dynamically their instability through a reduction of the possible number of coexisting species, i.e. diversity, and eventually reach a marginally stable state saturating May’s bound. Since this phenomenon stems from a dynamical process, it holds for a broad range of system parameters. It is robust against a range of variations in the model, including noise and different functional forms of responses and interactions (see figure 1 for two possible forms of the dynamics, and figure 6 for additional functional forms). Although in many physical cases criticality emerges only at phase transitions, i.e. for very special values of the parameters, there also exist critical phases of matter which instead cover a wide portion of their phase diagram. By relating the LV model to systems studied in condensed matter physics, the multiple equilibria regime is shown to be akin to a critical spin-glass phase. This connection to a phase that exists in a broad range of situations suggests that the applicability of our results goes well beyond the LV model and it offers a possible explanation of why so many different systems are found at the edge of stability: they are in a critical marginally stable phase. It also makes clear that this result, while general, is expected to have a well-defined regime of validity, as we shall explain at the end of this work. Finally it makes predictions on some distinctive features of the dynamical behavior of ecological systems at criticality, which will be interesting to test.

Zoom In Zoom Out Reset image size

The LV model we focus on is defined as follows. There are S species in the regional pool, whose abundance is N_i ≥ 0. The dynamical equations read

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{N}}_{i}}{{\rm{d}}{t}}=\displaystyle \frac{{r}_{i}}{{K}_{i}}{N}_{i}({K}_{i}-{N}_{i})-{N}_{i}\displaystyle \sum _{j,(j\ne i)}{\alpha }_{{ij}}{N}_{j}+\sqrt{{N}_{i}}{\eta }_{i}(t)+\lambda ,\end{eqnarray} \tag{ 1 }$$

where r_i is the intrinsic growth rate of species i, and K_i is the carrying capacity. It corresponds to the equilibrium abundance to which species i would self-regulate in absence of interaction. For sake of clarity, in the following we focus on the case where r_i and K_i are constants (that we set equal to one by rescaling $N\to {NK},t\to t/{rK}$ and absorbing those changes in redefinition of α_ij, η_i(t), λ). Later, we shall consider the effect of variability in r_i and K_i, and also different functional responses by replacing ${N}_{i}({K}_{i}-{N}_{i})$ with more general forms, such as ${N}_{i}f({N}_{i})$, and models with nonlinear interactions. In equation (1) the interaction between species is encoded in the matrix α_ij. We also add a small (infinitesimal) immigration rate λ to ensure that all invadable species exist, thus avoiding absorbing states that are unstable under small migration. Note that a species dependent λ_i would not change the results of our analysis. Finally, η_i(t) is a white noise with variance 2ω², and $\sqrt{{N}_{i}}$ captures the scaling of the noise due to the population size. We chose to use Ito’s convention for the multiplicative noise since it correctly captures the fact that a species with N_i = 0 remains at zero abundance also in presence of noise. We consider a symmetric interaction matrix α_ij = α_ji, corresponding to competitive (or weakly mutualistic) interactions; we will discuss in the conclusion the effect of asymmetry. Except for this constraint, no additional structure is incorporated, and the entries α_ij are taken to be independent identically distributed random variables. This provides a null-model, to which additional structure, such as trophic levels or space, can be subsequently added. Note that, as already anticipated above, the assumption on the randomness is done at the level of the pool and not of the community. The random variable α_ij can be drawn from any distribution without long tails, all that matters are its mean and variance. It turns out that the parameters that play a role in the final theory are the average number of links, C, per site and the first two moments of α_ij though the combination μ = C mean [α_ij ] and σ² = C var[α_ij]. For the sake of clarity, we now focus on the case C = S in which all species interact, extensions are discussed at the end of this paper⁵ .

The LV equations (1) can be rewritten in a way that makes their relationship with stochastic equations studied in physics more transparent:

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{N}}_{i}}{{\rm{d}}{t}}=-{N}_{i}\left[{{\rm{\nabla }}}_{{N}_{i}}{V}_{i}({N}_{i})+\displaystyle \sum _{j,(j\ne i)}{\alpha }_{{ij}}{N}_{j}\right]+\sqrt{{N}_{i}}{\eta }_{i}(t)+\lambda ,\end{eqnarray} \tag{ 2 }$$

where the ‘potential’ V_i(N_i) is equal to $\tfrac{{r}_{i}}{{K}_{i}}(-{K}_{i}{N}_{i}+{N}_{i}^{2}/2)$. Without noise these equations admit a Lyapounov function, i.e. a function whose value increases in the dynamics [32], dL/dt ≥ 0, with

$$\begin{eqnarray*}&&L=-\displaystyle \sum _{i}{V}_{i}({N}_{i})-\displaystyle \frac{1}{2}\displaystyle \sum _{i\ne j}{\alpha }_{{ij}}{N}_{i}{N}_{j}+\lambda \displaystyle \sum _{i}\mathrm{log}{N}_{i}.\end{eqnarray*}$$

In presence of noise, equations (2) are generalized Langevin equations. In the SI we show that they represent equilibrium dynamics of a thermal system with temperature T = ω² and characterized by the following effective Hamiltonian, or energy, $H=-L+{\sum }_{i}T\mathrm{log}{N}_{i}$. As a consequence, the long-time stationary probability distribution is the Boltzmann law: P = e^−H/T/Z, where the partition function Z guarantees the normalization. This result reveals that understanding the equilibria and the dynamics associated with the LV equations (1) can be exactly reformulated as a problem of statistical physics of thermal disordered systems, in which the N_i represent the degrees of freedom interacting via random couplings α_ij and subjected to individual potentials ${V}_{i}({N}_{i})+(T-\lambda )\mathrm{log}{N}_{i}$. Therefore, we can deduce properties of the equilibria reached dynamically by a thermodynamical analysis. In particular, without noise, i.e. at zero temperature, the equilibria correspond to the minima of the energy function, i.e. to the ground state and the metastable states In order to obtain the full solution of this model we shall use tools developed in statistical physics of disordered systems, including replica computations, discussed in detail in the SI. Analogous strategies have been already used in the past in similar contexts. Equilibrium thermodynamics was used to study the related replicator equations without noise in [2, 3]. More recently, in [5], the LV dynamics were studied in the context of the niche versus neutral debate. An approximate Hamiltonian was suggested, and the analytical solution of the random energy model was employed. Note, however, that in the present context such approximations would not allow to obtain the marginal stability we report here.

The phase diagram without noise, i.e. at zero temperature, and for small migration, $\lambda \to {0}^{+}$, had been obtained beforehand [9]. We reproduce it in the SI for completeness (supplementary figure 1 available online at stacks.iop.org/NJP/20/083051/mmedia) and show that it coincides, as it should, with the one obtained by the replica method. One finds that when interaction strengths are all identical, all species coexist in the community when μ > 0. (The range −1 < μ < 0 is not of interest for this work because the multiple equilibria phase we want to study is not present.) By increasing the variability in the interactions, more and more species are driven out of the community. They are characterized by N_i = 0 for $\lambda \to {0}^{+}$ (we will call them ‘extinct’ henceforth). The equilibrium reached dynamically is stable under perturbations, that is by changing ${V}_{i}({N}_{i})\to {V}_{i}({N}_{i})-{\xi }_{i}{N}_{i}$, and there is a gap in the spectrum of the corresponding stability matrix [9], see figure 1(A). (This matrix, determining the stability under changes that affect the carrying capacities, is different from the one governing stability under the demographic noise η_i.) Note that in this regime the final community composition is unique, independent of the assembly history, e.g. the initial conditions for the dynamics. This picture persists up to ${\sigma }_{c}=1/\sqrt{2}$. By increasing randomness in the interactions above σ_c a transition to another regime, which is sharp for large systems, takes place. Our purpose is to study the phase reached when crossing the transition. Henceforth we continue to focus on the zero-noise case, the effect of the demographic noise is discussed at the end of the paper.

In physics terms, the single equilibrium regime corresponds to a ‘paramagnetic phase’ where the zero temperature values of the degrees of freedom N_i are mainly fixed by the external potential V_i(N_i) and a mean-field anti-ferromagnetic (competitive) interaction. By increasing the randomness in the α_ij variables the system undergoes a zero temperature phase transition toward a spin-glass phase, characterized by many local minima of the energy (or maxima of the Lyapounov) function and, hence, multiple equilibria. We have used the replica method to study it (see SI) and found that the regime with multiple equilibria corresponds, technically, to a full replica symmetry breaking (RSB) solution. On the basis of all previous analysis of mean-field spin-glasses [33, 34], we can then make general statements about the regime with multiple equilibria. Note that here the term mean-field refers to the fact that the underlying interaction network is fully connected, and not (as often used in ecology) that all interactions are identical in strength.

First, it is characterized by a large number of equilibria. These equilibria are minima of the energy, separated by regions with higher energies that form what are called barriers. The lowest equilibria are typically separated by barriers that diverge in the large S limit, while the higher ones by barriers of order of one, i.e. that do not scale with S [35]. Second, and central to our discussion, all minima display a stability matrix characterized by arbitrary small eigenvalues for large S, i.e. minima are marginally stable and characterized by flat directions in the energy landscape at quadratic order. The ground state has this property, and also the higher energy local minima which are usually less stable.

We now explain the main findings of our thermodynamic analysis and relate them to random matrix theory results, see SI for details. The two main observables we focus on are: (i) $\tfrac{\partial {N}_{i}^{* }}{\partial {\xi }_{i}^{* }}$, which is the response of a single species to a perturbation $V({N}_{i})\to V({N}_{i})-{\xi }_{i}{N}_{i}$, where the star indicates that only non-extinct species are considered and (ii) $\phi \equiv {S}^{* }/S$, which is the fraction of species present in the community, called diversity in what follows. For identical interaction strengths, i.e. at σ = 0, all species coexist (ϕ = 1) and $\tfrac{\partial {N}_{i}^{* }}{\partial {\xi }_{i}^{* }}=1$. Increasing σ we find that $\tfrac{\partial {N}_{i}^{* }}{\partial {\xi }_{i}^{* }}$ is constant across species and increasing. Concomitantly, the diversity decreases. As found beforehand [9], at ${\sigma }_{c}=1/\sqrt{2}$ the system undergoes a sharp transition from the single to the multiple equilibria regime. This corresponds physically to a phase transition to the spin-glass phase. In this phase we find that for all equilibria and all species i,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {N}_{i}^{* }}{\partial {\xi }_{i}^{* }}=\sqrt{\displaystyle \frac{1}{{\phi }_{c}{\sigma }_{c}^{2}}}=2\,\,\,\mathrm{and}\,\,\,\phi {\sigma }^{2}={\phi }_{c}{\sigma }_{c}^{2}=\displaystyle \frac{1}{4},\end{eqnarray} \tag{ 3 }$$

where ϕ_c is the value of the diversity at the transition. These identities hold throughout the entire multiple equilibria phase and they are consequences of the criticality of the spin-glass phase. We also find that for σ > σ_c the number of equilibria, i.e. of energy landscape minima, is exponential in S. In fact, the number of equilibria scales as e^hS where h goes to zero at the transition from one to multiple equilibria. One cannot make general statements on its order of magnitude, since it depends on the external parameters and on the model, in particular the functional response f(N). When comparing to numerical and experimental results, it is important to keep in mind that it can be small, as we found for instance for the standard LV model. In consequence even for rather large S ∼ 100 the number of equilibria may be modest, see SI. Figures 2 and 3 confirm our analytical predictions by showing, respectively, numerical results for $\tfrac{\partial {N}_{i}^{* }}{\partial {\xi }_{i}^{* }}$ and ϕ corresponding to a given equilibrium reached dynamically.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The second identity in equation (3) corresponds to a saturated form of May’s original bound, $\mathrm{var}({\alpha }_{{ij}})=1/(4{S}^{* })=1/(4{C}^{* })$, where C^* is the average number of interactions per surviving species [1] (the prefactor 4 comes from the symmetry of α_ij in the present model). In order to reveal this connection with random matrix theory we focus on the ${S}^{* }\times {S}^{* }$ stability matrix M^* associated to a given equilibrium, defined by the relation $({M}^{* }{)}_{{ij}}^{-1}=\tfrac{\partial {N}_{i}^{* }}{\partial {\xi }_{j}^{* }}$. Using the fixed point equation corresponding to equation (1) it is easy to check that

$$\begin{eqnarray}&&{M}_{{ij}}^{* }={\delta }_{{ij}}+{\alpha }_{{ij}},\end{eqnarray} \tag{ 4 }$$

where δ_ij is the Kronecker delta. In this equation the indices i, j have to be reduced to the surviving species since extinct species remain so if one adds an infinitesimal perturbation ξ_j, and do not contribute to the stability of the equilibrium reached dynamically In fact, in the limit of small migration, $\lambda \to {0}^{+}$, extinct species are those that cannot invade: ${{\rm{\nabla }}}_{{N}_{i}}{V}_{i}(0)+{\sum }_{j(\ne i)}{\alpha }_{{ij}}{N}_{j}^{* }\gt 0$. Adding an infinitesimal ξ_i does not change this property and thus the species remains extinct. Following procedures developed for mean-field spin-glasses [36], one can show (see SI) that the spectrum of ${M}_{{ij}}^{* }$ is identical for large S to the one of a ${S}^{* }\times {S}^{* }$ matrix with independent identically distributed Gaussian off-diagonal entries having the same first and second moment of α_ij. This is by no means trivial since the equilibrium reached dynamically [37], and hence the identity of the surviving species, depend on α_ij and induce correlations in the off-diagonal elements of ${M}_{{ij}}^{* }$. The relation with random matrices implies that the eigenvalue density of the stability matrix is a Wigner semi-circle with support $[-2\sigma \sqrt{\phi }+1,2\sigma \sqrt{\phi }+1]$, as we indeed find numerically, see figure 1. Moreover, this directly connects (3), which holds in the entire spin-glass phase, to marginal stability. We have therefore recovered May’s original stability bound but in a saturated form: the number of surviving species, Sϕ, is exactly the one guaranteeing that the system is poised at the edge of stability, similarly to what was proposed in the self-organized instability scenario [22].

Let us now discuss extensions and the range of validity of our results. We have verified that our conclusions on the multiple equilibria regime continue to hold for several different convex functional responses f(N), the standard $f(N)=N(1-N)$ being only an example among others, and with variability in the values of r_i and K_i. This is a direct consequence of the properties of the spin-glass phase to which the multiple equilibria regime is related to. In these more general cases, the critical character of this phase is encoded in the following identity valid for the average of the square of the single species response in the whole multiple equilibria regime:

$$\begin{eqnarray}&&\phi {\sigma }^{2}\left(\displaystyle \frac{1}{{S}^{* }}\displaystyle \sum _{i=1}^{{S}^{* }}{\left(\displaystyle \frac{\partial {N}_{i}^{* }}{\partial {\xi }_{i}^{* }}\right)}^{2}\right)=1.\end{eqnarray} \tag{ 5 }$$

Note that equation (5) reduces to the first identity in equation (3) when single species response are identical, as previously discussed. We show in figure 4 the numerical results obtained for $f(N)=1-N-{(N-1)}^{2}/4$ confirming this prediction: the rhs of equation (5) is less than one in the single equilibrium phase, it reaches one at the transition and then remains stuck to this value in the whole multiple equilibria phase. As before, the criticality of the spin-glass phase implies marginal stability. Indeed, similarly to the standard LV case, the spectrum of the stability matrix ${M}_{{ij}}^{* }=V^{\prime\prime} ({N}_{i}^{* }){\delta }_{{ij}}+{\alpha }_{{ij}}$ can be shown to be identical for large S to that of a ${S}^{* }\times {S}^{* }$ matrix with independent identically distributed Gaussian off-diagonal entries having the same first and second moment of α_ij, and independent identically distributed diagonal entries with the same statistics of $V^{\prime\prime} ({N}_{i}^{* })$. As we show in the SI, the condition that the left edge of the eigenvalues density touches zero for this class of random matrices is $\phi {\sigma }^{2}\left(\tfrac{1}{{S}^{* }}{\sum }_{i=1}^{{S}^{* }}{(({M}^{* }{)}_{{ii}}^{-1})}^{2}\right)=1$, which using $\tfrac{\partial {N}_{i}^{* }}{\partial {\xi }_{i}^{* }}=({M}^{* }{)}_{{ii}}^{-1}$ turns out to be identical to equation (5). Therefore we obtain that the multiple equilibria regime is indeed generically characterized by marginal stability and, by doing so, we derive a new generalized version of May’s bound (equation (5)). Remarkably, these properties hold despite the fact that for this general class of random matrices the density of eigenvalues is no longer a shifted semi-circle and the singularity at the left edge is not necessarily a square root, as shown numerically in figure 1 for the $f{(N)=1-N-3/4(N-1)}^{2}$ case (which has $f^{\prime} (N)\gt 0$ for some values of N, corresponding to an Allee effect). In particular, the singularity at the edge depends on the statistics of the diagonal components ${M}_{{ii}}^{* }$, i.e. of $V^{\prime\prime} ({N}_{i}^{* })$. It is a square root if the distribution of the random variable $V({N}_{i}^{* })$ approaches the left edge of the support slower than linearly, in the other cases it may be a square root or instead inherit the singularity of $V^{\prime\prime} ({N}_{i}^{* })$ at the left edge [38]. Note that for general f(N) the phase diagram is modified compared to the standard LV model. For instance for $f{(N)=1-N-3/4(N-1)}^{2}$ the single equilibrium phase is absent even for infinitesimally small interactions. In conclusion, our investigations show that the class of f(N) leading to the marginally stable multiple equilibria phase, i.e. the phase boundaries of the critical spin-glass phase in the physics terminology, is quite large. Determining its boundaries is an important and interesting task that we leave for future studies. Based on previous results on mean-field glassy systems [39], it is possible that the property that the multiple equilibria probed by the system are marginally stable is robust, even though the detailed properties of the landscape might be different depending on the shape of f(N). This is due to the fact that in many different situations, the most numerous minima, which are the ones with the biggest basins of attraction, are marginally stable [39]. Therefore it is not necessary that all minima are marginally stable to find that the ones reached dynamically are like this. Another extension of our work worthy of future analysis concerns the role of the interactions network. As long as the connectivity C per species is large and the underlying structure rather homogeneous, e.g. no fat tails in the distribution of the local connectivity, the mean-field approach we developed is a very good approximation. In consequence, our results are expected to hold also in this more general setting where C ≫ 1 but C is not necessarily equal to S. As a first interesting extension one could consider LV models on random regular graphs [40].

Zoom In Zoom Out Reset image size

The properties (and the existence) of the multi-equilibria phase continue to hold also in presence of small noise. On the basis of previous studies on mean-field spin-glasses [34, 35], we can state several general facts. In presence of small noise the system moves around between multiple dynamical states. These correspond to low temperature spin-glass states associated with the local minima discussed above. Assuming that the free-energy landscape is akin to the one found for mean-field spin-glasses, e.g. the Sherrington–Kirkpatrick (SK) model [41] (more complicated free-energy landscapes have been found in generalized glass models [42]), we expect that only the low energy minima are able to trap the dynamics of the ecosystem over long periods of time, while the ones higher in energy are instead separated by small energy barriers; the transitions between them are so frequent that their identities as separate states disappear even for small noise. The stability of these dynamical states can described by the matrix ${{ \mathcal M }}_{{ij}}$, which is a generalization of the stability matrix and is defined as the (matrix) inverse of $\tfrac{\partial \langle {N}_{i}\rangle }{\partial {\xi }_{j}}$ ($\langle \cdot \rangle$ denotes the average over the noise). As the stability matrix in absence of noise, ${{ \mathcal M }}_{{ij}}$ is positive definite and has arbitrary small eigenvalues for large S, thus leading to marginal stability. Physically, this is a consequence of the fact that the spin-glass phase is not destroyed by small thermal fluctuations and is critical over a finite range of temperatures [33].

In summary, by mapping the LV models to thermal disordered systems and studying their thermodynamics, we find that marginal stability is a property of all communities that are reached dynamically by an ecosystem in the multiple equilibria phase. Equation (5), combined with our random matrix analysis, relates this property to the single species response. This is the main result of our work: it represents an exact statement of May’s stability bound [1], with three notable differences: 1. it follows from an exact analysis of the communities reached dynamically rather than from a priori assumptions on the stability equations, it is thus a property of the emergent community; 2. it is saturated, with an equality rather than an inequality; 3. it is more general, allowing to incorporate nonlinear $f(N)$.

Having established that marginal stability is a generic property of the multiple equilibria phase, we now discuss some of its consequences and propose measurable tests. The most striking effects are expected to appear in dynamical phenomena. Again, previous results on dynamics of mean-field spin glasses provide useful guidelines [34]. In particular, starting the LV dynamics from random initial conditions one expects slow relaxations toward the minima and history dependence for large S. Both phenomena are tightly linked to marginal stability which results in flat directions in configurations space. The response to perturbations is also expected to be very unusual: marginal stability should lead to strong and wildly fluctuating nonlinear responses [43] and avalanches of extinctions and invasions [11]. Working out the relevance of this phenomena in various ecological contexts certainly warrants future research. Note that, these avalanches of extinctions are different from cascades in trophic systems [44], as here no trophic structure is included. The results we found for symmetric interactions have also important consequences for cases where α_ij are asymmetric. Indeed, given that the multiple equilibria are marginally stable, we expect that adding asymmetry leads immediately to a chaotic behavior in which the system moves among the different regions of configurations space corresponding to the vestige of those equilibria [45–47]. Chaotic dynamics of LV equations have been observed in simulations in some region in parameter space in [7], and may be related to the phase we describe here.

Throughout this paper, we stressed that the unusual properties of the multiple equilibria phase are related to the criticality of the corresponding spin-glass phase. In the following, we show that this relationship also suggests new ways to test for marginal stability. Criticality corresponds to a state in which the microscopic degrees of freedom are all strongly correlated, which naturally leads to singular responses. The properties of the stability matrix in the multiple equilibria phase are one facet of this phenomenon; diverging fluctuations are another. Note, however, that simple fluctuations such as $\overline{\langle {N}_{i}^{2}\rangle -\langle {N}_{i}{\rangle }^{2}}$ or its time-dependent counterpart $C(t,t^{\prime} )=\overline{\langle {N}_{i}(t){N}_{i}(t^{\prime} )\rangle -\langle {N}_{i}(t)\rangle \langle {N}_{i}(t^{\prime} )\rangle }$ do not capture criticality [33] (the overbar denotes the average across the species). One needs to probe large scale fluctuations, hence inter-species correlations (this is the counterpart of correlations between different points in space in standard critical phenomena). As it was understood in the context of disordered systems, a good probe for this are four-points correlations [33] (see also SI for detailed computations on this point and the following). As a matter of fact, in the LV model we considered, which can be mapped onto thermal dynamics, diverging responses and fluctuations are exactly related by the fluctuation–dissipation relation $\tfrac{\partial \langle {N}_{i}\rangle }{\partial {\xi }_{j}}={\omega }^{-2}(\langle {N}_{i}{N}_{j}\rangle -\langle {N}_{i}\rangle \langle {N}_{j}\rangle )$. On this basis, and following previous work on glassy systems [48], we propose to probe criticality, or the collective nature of the equilibria and more generally of the dynamical states by looking at a fourth-order correlation function, called ${\chi }_{4}(t,t^{\prime} )$, which reads:

$$\begin{eqnarray*}&&{\chi }_{4}(t,t^{\prime} )=\displaystyle \frac{1}{S}\displaystyle \sum _{{ij}}[\langle \delta {N}_{i}(t)\delta {N}_{i}(t^{\prime} )\delta {N}_{j}(t)\delta {N}_{j}(t^{\prime} )\rangle -\langle \delta {N}_{i}(t)\delta {N}_{i}(t^{\prime} )\rangle \langle \delta {N}_{j}(t)\delta {N}_{j}(t^{\prime} )\rangle ],\end{eqnarray*}$$

where $\delta {N}_{i}(t)=({N}_{i}(t)-\langle {N}_{i}(t)\rangle )/\sqrt{C(t,t)}$. This function allows one to measure to what extent species are dynamically correlated and is therefore a way to quantitatively test criticality and marginal stability. In the LV model and for σ < σ_c the dynamics becomes stationary after a short transient. In this case ${\chi }_{4}(t,t^{\prime} )$ depends on $t-t^{\prime}$ only. Its long time limit ($t-t^{\prime} \gg 1$) is equal to:

$$\begin{eqnarray*}&&{\chi }_{{\rm{SG}}}=\displaystyle \frac{1}{S{\overline{\langle {N}_{i}^{2}{\rangle }_{c}}}^{2}}\displaystyle \sum _{i,j}{[\langle {N}_{i}{N}_{j}\rangle -\langle {N}_{i}\rangle \langle {N}_{j}\rangle ]}^{2},\end{eqnarray*}$$

whereas for $t=t^{\prime}$ and small temperatures χ₄(t, t) is equal to twice χ_SG. As shown in figure 5, where we check these analytical predictions by numerical simulations for S = 400, ${\chi }_{4}(t-t^{\prime} )$ is a decreasing function of $t-t^{\prime}$. The shape of ${\chi }_{4}(t,t^{\prime} )$ is different from the one found in glassy systems. The main difference which lies in the different nature of the degrees of freedom is that χ₄(t, t) becomes large approaching the transition, whereas in glassy systems it is featureless. In the condensed matter theory context ${\chi }_{{\rm{SG}}}$ is known as the spin-glass susceptibility and is known to diverge in the entire spin-glass phase. We indeed recover this result and connect it to marginal stability since

$$\begin{eqnarray}&&{\chi }_{{\rm{SG}}}={\left(\displaystyle \frac{\partial \langle {N}_{i}\rangle }{\partial {\xi }_{i}}\right)}^{-2}\int {\rm{d}}\lambda \displaystyle \frac{\rho (\lambda )}{{\lambda }^{2}}\sim {S}^{1/3},\end{eqnarray} \tag{ 6 }$$

where ρ(λ) is the density of eigenvalues of the matrix ${{ \mathcal M }}_{{ij}}$. The divergence of χ_SG as S^1/3 comes from the fact that the minimum eigenvalue of ${ \mathcal M }$ scales [49] as S^−2/3. In the inset of figure 5 we show the behavior of χ_SG obtained analytically in the large S limit. Note that if the singularity of ρ(λ) is milder than a square root, as it is the case for example for $f{(N)=1-N-3/4(N-1)}^{2}$, then one needs to consider high-order moments. The bottom line of this discussion is that if data on the time-dependence of abundances is available, the function ${\chi }_{4}(t,t^{\prime} )$ allows one to measure to what extent species dynamics are correlated and test directly for criticality and marginal stability. Even though in the simple LV case, measuring χ_SG would be sufficient for that purpose, measuring the time dependent four-point function ${\chi }_{4}(t,t^{\prime} )$ is the way to go in order to obtain information in more general cases, which may be neither stationary, nor related to thermal equilibrium dynamics. In particular, ${\chi }_{4}(t,t^{\prime} )$ can show history dependence of dynamic correlations [34, 48], by being a function of t and $t^{\prime}$ separately, and not only of $t-t^{\prime}$.

Zoom In Zoom Out Reset image size

These results provide measurable predictions useful to probe whether the system is at (or close to) marginality: First, the single-species response is predicted to follow equation (3) or more generally equation (5). Experimentally, such perturbations may be measured using press perturbation experiments [50]. Secondly, the fourth moment ${\chi }_{4}$ defined above is expected to be very large, diverging as S^1/3, see equation (6), whereas in the Unique Equilibrium phase it does not grow with the diversity S. Note that χ₄ can be computed directly from species abundance data.

In conclusion, our analysis of LV equations in the limit of large species shows a mechanism for which many systems in Nature are poised at the edge of stability: we have shown that when the parameters of an ecosystem cross the limit of stability, the system dynamically self-adapts to remain exactly marginally stable. It does so reducing the number of species in such a way to saturate May’s bound, which therefore emerges as a result of a dynamical process. This leads to a whole critical phase with multiple marginally stable equilibria, which is expected to be present for several different models and to display highly non-trivial dynamical behaviors that may be measured experimentally. Its consequences can be relevant and important in many fields [13, 16, 20, 51, 52].

Numerics

Models with nonlinear interactions

In this section we show that a marginal phase can also be found in cases where the interactions assume a nonlinear form. As in figure 1(B), and in contrast to figure 1(A), the distribution of eigenvalues of the equilibria reached dynamically by the system, is seen to touch the y-axis, indicating marginal stability. The two studied models are defined by

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{N}}_{i}}{{\rm{d}}{t}}={r}_{i}{N}_{i}\left\{1-{N}_{i}-g\left[\displaystyle \sum _{j}{\alpha }_{{ij}}f({N}_{j})\right]\right\},\end{eqnarray} \tag{ 7 }$$

where two different forms of $g(x)$ are taken into account. In the model in panel (A) of figure 6, $g(x)$ is taken to have a Holling type II form $g(x)=5x/(1+| x| )$ and $f(x)=x$. The model in panel (B) is defined via $f(x)=x/(1+x)$ and $g(x)=x$. The matrix α is symmetric with $\langle {\alpha }_{{ij}}\rangle =0.25,\sigma =0.7$ for panel (B) and $\langle {\alpha }_{{ij}}\rangle =0.125,\sigma =1.3$ for panel (A). All r_i = 1, and S = 400. In both cases ρ(λ) has a support whose left edge touches zero, corresponding to marginal stability, even when its shape is different compared to the shifted Wigner semi-circle of a standard LV model (figure 1). The parameters have not been fine-tuned, and can be varied within some limits—which define the phase boundaries—while maintaining the marginality property.

Zoom In Zoom Out Reset image size

Spin glasses

Many of the theoretical methods used in this work were originally developed to study disordered systems in physics, and in particular spin-glasses. Here we make some brief comments on such systems and how they relate to the present problem.

The behavior described in this work requires variability in the interaction strengths α_ij (as measured by their standard deviation σ). In physics, systems where interactions between the constituents exhibit analogous variations are known as ‘disordered systems’. In particular, a spin glasses is systems where magnetic interactions vary between the different pairs of atoms. Models of such systems, starting with [53], traditionally use binary variables to model the state of each magnetic spin σ ∈ {−1, 1}, while the complex interactions were modeled as i.i.d. Gaussian random variables. For example, the first spin glass model [53] is the Edwards and Anderson model where the Hamiltonian (the energy of a state) is given by

$$\begin{eqnarray*}&&{H}_{\mathrm{SK}}=-\displaystyle \sum _{(i,j)}{J}_{i,j}{\sigma }_{i}{\sigma }_{j},\end{eqnarray*}$$

where only terms involving pairs of nearest neighbors (i, j) are included in the sum, and J_i,j are randomly sampled with zero mean. The first model to be solved is its fully connected analog, the SK model where all pairs of spins are taken to interact [41, 54]. A number of other spin glass models were considered along the years [33, 55–57]. Among the most interesting ones, we mention those where interactions involve p > 2 variables at a time, called p-spin models [55]

$$\begin{eqnarray*}&&{H}_{p}=-\displaystyle \sum _{({i}_{1},{i}_{2},\ldots ,{i}_{p})}{J}_{{i}_{1},{i}_{2},\ldots ,{i}_{p}}{\sigma }_{{i}_{1}}{\sigma }_{{i}_{2}}...{\sigma }_{{i}_{2}}.\end{eqnarray*}$$

It turned out that spin glass models show interesting new phases and phase transitions and can be classified in different universality classes which correspond to different macroscopic behaviors.

The techniques developed for the solution of spin glasses are in general useful to describe systems with high level of frustration [58] (i.e. absence of optimal solutions for certain instances of the couplings). For this reason they are widely applied nowadays in different fields including condensed matter (magnetic systems, supercooled liquids), biology (protein folding, neural networks), social sciences, economics, computer sciences (optimization theory, machine learning).

The major obstacle that had to be tackled while dealing with the solution of spin glass models is represented by the task of performing disorder averages of quantities of physical interest to extract information on the macroscopic behavior of the system. The information about equilibrium is contained in the so-called free energy of the system

$$\begin{eqnarray*}&&F=-\displaystyle \frac{1}{\beta }{\overline{\mathrm{log}\displaystyle \sum _{\{{\sigma }_{i}\}}\exp [-\beta H]}}^{J},\end{eqnarray*}$$

where the overline represents the average over difference instances of the disordered couplings, the inner sum runs over all the possible configurations of spins, and the argument of the logarithm is commonly called partition function

$$\begin{eqnarray*}&&Z=\displaystyle \sum _{\{{\sigma }_{i}\}}\exp [-\beta H].\end{eqnarray*}$$

The operation of taking the average usually is reduced to perform a Gaussian integration. This would have required little computation effort indeed had the logarithm not been on the way. Yet its presence cannot be neglected nor the operation of taking the average simply performed on the logarithm’s argument (the last procedures is called annealed calculation but does not lead to the correct solution in the interesting regimes). To keep the order of the operations and yet end up with an analytically tractable problem the so-called replica trick was introduced [33]. It amounts to use the fact that

$$\begin{eqnarray*}&&\overline{\mathrm{log}x}=\mathop{\mathrm{lim}}\limits_{n\to 0}\displaystyle \frac{\mathrm{log}\overline{{x}^{n}}}{n},\end{eqnarray*}$$

which can be easily verified. Hence the free energy can be written in terms of the so-called replicated partition function Zⁿ as

$$\begin{eqnarray*}&&-\beta F=\mathop{\mathrm{lim}}\limits_{n\to 0}\displaystyle \frac{\mathrm{log}\overline{{Z}^{n}}}{n}=\mathop{\mathrm{lim}}\limits_{n\to 0}\displaystyle \frac{\mathrm{log}\overline{\{{\displaystyle \sum }_{\{{\sigma }_{i}\}}\exp [-\beta H]\}{}^{n}}}{n},\end{eqnarray*}$$

where the power n can be interpreted, before the limit $n\to 0$ is taken, as we were focusing on n independent copies of the same system in presence of a unique sample of random couplings. Averages over different realizations of the disorder are in this form straightforward. The copies in the spin glass jargon are usually called replicas.

To give an intuition of the physical meaning of the results that can be obtained within replica computations we must remember that frustrated systems are usually characterized by a multi minima structure of the energy, or any equivalent cost function that might be of interest. This arrangement of minima is uniquely associated to any instance of the random couplings. The role of replicas is the one of revealing the main features of this multi minima structure by independently probing different minima. In fact one of the most important piece of information that comes out from a replica computation is the average width of equilibrium minima and the average distance between pairs of them, or more in general the hierarchical arrangements of minima in the space of configurations. All this is contained in the structure of overlap Q_a,b (or similarity, which accounts for the inverse of distance in the space of configurations) between pairs of replicas a, b:

$$\begin{eqnarray*}&&{Q}_{a,b}=\displaystyle \frac{1}{N}\displaystyle \sum _{i}^{N}{\overline{\displaystyle \frac{{\displaystyle \sum }_{\{{\sigma }_{i}^{c}\}}\exp \left[-\beta {\displaystyle \sum }_{c}^{n}{H}_{c}\right]{\sigma }_{i}^{a}{\sigma }_{i}^{b}}{{\displaystyle \sum }_{\{{\sigma }_{i}^{c}\}}\exp \left[-\beta {\displaystyle \sum }_{c}^{n}{H}_{c}\right]}}}^{J}=\displaystyle \frac{1}{N}\displaystyle \sum _{i}^{N}{\overline{\langle {\sigma }_{i}^{a}{\sigma }_{i}^{b}{\rangle }_{\mathrm{AR}}}}^{J},\end{eqnarray*}$$

where $\langle {\rangle }_{{\rm{AR}}}$ denotes the measure over all replicas with the Boltzmann weight and after averaging over disorder (rhs of the first equation above). Along a typical replica computation a new conceptual obstacle arises when the free energy is rewritten in terms of an integral over all the possible choices of the overlap matrix. In the thermodynamic limit ($N\to \infty$) the Laplace method (or saddle point method) can be applied to evaluate the integral but it requires to find the overlap matrix that maximizes the argument of the integral. This operation requires the introduction of a good ansatz for the overlap matrix. The currently used scheme was proposed by Parisi [59, 60] and subsequently proved to be the one providing the correct saddle point result [61]. It is called RSB scheme and will be discussed in more details in the following sections.

Depending on the number of steps of breaking of the replica symmetry required to get a meaningful solution [62] (i.e. stable in the replica space), we could end up working with replica symmetric scheme, one step replica symmetry breaking, $\infty$ steps replica symmetry breaking (FRSB), just to mention the most relevant ones. This differentiation allows to classify spin glass models and characterize the features of their relevant phases and phase transition.

It turns out that the ecological model we consider in this work, at large values of σ, is characterized by a FRSB solution. The FRSB solution represents, as stressed in the main text, a critical phase. From the technical point of view, it is marginally stable, meaning that within the Laplace method it corresponds to an extremum with vanishing small eigenvalues.

We thank J-P Bouchaud, D S Fisher, J Kurchan, T Mora, V Ros, E Vanden-Eijnden, P Vivo, A Walczak for useful discussions. This work was partially supported by the grant from the Simons Foundation (♯454935, Giulio Biroli).

The authors declare that they have no competing financial interests.