Can endogenous fluctuations persist in high-diversity ecosystems?

TL;DR

This paper develops a theoretical framework to understand how high-diversity, spatially-extended ecosystems can sustain persistent endogenous chaotic fluctuations, which promote biodiversity and differ from well-mixed systems.

Contribution

It introduces a dynamical mean-field theory to predict conditions for endogenous fluctuations and explores their implications for biodiversity and ecosystem dynamics.

Findings

Endogenous fluctuations can persist for long times in high-diversity systems.

Fluctuation amplitude is closely linked to species diversity.

Spatial extension influences the nature and persistence of fluctuations.

Abstract

When can complex ecological interactions drive an entire ecosystem into a persistent non-equilibrium state, where species abundances keep fluctuating without going to extinction? We show that high-diversity spatially-extended systems, in which conditions vary somewhat between spatial locations, can exhibit chaotic dynamics which persist for extremely long times. We develop a theoretical framework, based on dynamical mean-field theory, to quantify the conditions under which these fluctuating states exist, and predict their properties. We uncover parallels with the persistence of externally-perturbed ecosystems, such as the role of perturbation strength, synchrony and correlation time. But uniquely to endogenous fluctuations, these properties arise from the species dynamics themselves, creating feedback loops between perturbation and response. A key result is that the fluctuation amplitude and species diversity are tightly linked, in particular fluctuations enable dramatically more species to coexist than at equilibrium in the very same system. Our findings highlight crucial differences between well-mixed and spatially-extended systems, with implications for experiments and their ability to reproduce natural dynamics. They shed light on the maintenance of biodiversity, and the strength and synchrony of fluctuations observed in natural systems.

Equations

$$\frac{d}{dt}{N}_{i,u}$$

$$+\sum _v{D}_{i,uv}\left(N_{i,v}-N_{i,u}\right).$$

$$\frac{d{N}_u}{dt}=N_u\left(N_u^{\ast }-N_u+{\xi }_u\right)+\sum _v{D}_{uv}\left(N_v-N_u\right).$$

$$⟨{\xi }_u\left(t\right){\xi }_u\left(t^{\prime }\right)⟩={\sigma }^2\left[C_{N,u}\left(t,t^{\prime }\right)-C_{N,u}^{\infty }\right]$$

$$C_{N,u}\left(t,t^{\prime }\right)=⟨N_u\left(t\right)N_u\left(t^{\prime }\right)⟩,$$

$$C_{N,u}^{\infty }=\underset{t-t^{\prime }\to \infty }{\lim }⟨N_u\left(t\right)N_u\left(t^{\prime }\right)⟩.$$

$$W\equiv \int dt{C}_{\xi }\left(t,t^{\prime }\right).$$

$$\mathbb{E}\left[a_{ij,u}{a}_{kl,v}\right]={\delta }_{ik}{\delta }_{jl}{\rho }_{uv}$$

$$\frac{d}{dt}{N}_{i,u}$$

$$\displaystyle\ \ \ +\sum_{v}D_{i,uv}\left(N_{i,v}-N_{i,u}\right)\$$

$$\frac{d}{dt}{N}_u$$

$$\displaystyle\ \ \ +\sum_{v}D_{uv}\left(N_{v}-N_{u}\right)\$$

$$\{\begin{array}{rl}& m_u(t)=\overline{N_u(t)}\\ & ⟨{\eta }_u(t){\eta }_v(t^{\prime })⟩={\sigma }^2{\rho }_{uv}\overline{N_u(t)N_v(t^{\prime })}\end{array}$$

$$\{\begin{array}{rl}& m_u=\overline{N_u(t)}\\ & C_{uv}^N(t,t^{\prime })=C_{uv}^N(t-t^{\prime })\end{array}$$

$$⟨z_u{z}_v⟩={\sigma }^2\underset{t-t^{\prime }\to \infty }{\lim }{C}_{uv}^N(t-t^{\prime })$$

$$\frac{d}{dt}{N}_u$$

$$\displaystyle\ \ \ +\sum_{v}D_{uv}\left(N_{v}-N_{u}\right)\$$

$$\{\begin{array}{rl}& m_u=\underset{t^{\prime }\gg 1}{\lim }\overline{N_u(t^{\prime })}\\ & ⟨z_u{z}_v⟩={\sigma }^2\underset{t\gg t^{\prime }\gg 1}{\lim }{C}_{uv}^N(t,t^{\prime })\\ & ⟨{\xi }_u(t){\xi }_v(t^{\prime })⟩={\delta }_{uv}{\sigma }^2\left[C_{uv}^N(t,t^{\prime })-\underset{t\gg t^{\prime }\gg 1}{\lim }{C}_{uv}^N(t,t^{\prime })\right]\end{array}$$

$${\partial }_t{x}_u=N_u^{\ast }-e^{x_u}-\sigma {\xi }_u+D\sum _v\left(e^{x_v-x_u}-1\right).$$

$$J=\frac{1}{2W}{\int }^{t_f}dt\sum _{u=1}^M{\xi }_u^2,$$

$$J_{\min }=\frac{2x_c}{W}{N}_1^{\ast }.$$

$$J_{\min }=\frac{2x_c}{W}M{N}_{\mathrm{eff}}^{\ast }.$$

$$w\equiv -\sqrt{\frac{1}{m}\sum _{u=1}^m{\left(N_u^{\ast }\right)}^2},$$

$$N_{\mathrm{eff}}^{\ast }\equiv \frac{1}{4M}\sum _{{\left\{u\right\}}_{+}}\frac{{\left(N_u^{\ast }-w\right)}^2}{w}.$$

$$\frac{1}{N_u}\frac{d{N}_u}{dt}$$

$$W\equiv \int dt{C}_{\xi }\left(t,t^{\prime }\right),$$

$$0=\overline{\frac{1}{N_u}\frac{d{N}_u}{dt}}=N_u^{\ast }-\overline{N_u}+O\left(D\right)$$

$$\frac{d{x}_u}{dt}=N_u^{\ast }+\xi \left(t\right)\Rightarrow P\left(x\right)\sim e^{\frac{2N^{\ast }x}{{\sigma }^{2}W}}$$

$$\overline{e^{-x_u}}$$

$$\mathbb{E}\left[a_{ij,u}{a}_{kl,v}\right]={\delta }_{ik}{\delta }_{jl}{\rho }_{uv}$$

$$N_{i,u}^{\ast }=1-\mu m_u-\sigma S^{-1/2}\sum _{j,+}{a}_{ij,u}{N}_{j,u}^{\ast }$$

$$\{\begin{array}{l}\mathrm{mean}\left[N_u^{\ast }\right]=1-\mu {⟨N_{i,u}^{\ast }⟩}_{+}\\ \mathrm{covariance}\left[N_u^{\ast },N_v^{\ast }\right]={\sigma }^2{\rho }_{uv}{⟨N_{i,u}^{\ast }{N}_{i,v}^{\ast }⟩}_{+}\end{array}$$

Full text

While large temporal variations are widespread in natural populations lundberg_population_2000 ; inchausti_relation_2003 , it is difficult to ascertain how much they are caused by external perturbations, or by the ecosystem’s internal dynamics, see e.g. (ellner_chaos_1995, ; scheffer_why_2003, ). In particular, both theoretical tools and empirical results come short of addressing a fundamental question: can we identify when fluctuations in species abundances arise from complex ecological interactions?

Our focus here is on high-diversity communities. Historically, studies of endogenous fluctuations have focused on single populations or few species may_biological_1974 ; allen_chaos_1993 . On the other hand, theories of many-species interaction networks often center on ecosystems that return to equilibrium in the absence of perturbations (may_will_1972, ). Some authors have even proposed that endogenous fluctuations are generally too rare or short-lived to matter, since they can be self-defeating: dynamics that lead to large erratic variations cause extinctions, leaving only species whose interactions are less destabilizing, with extinctions continuing until an equilibrium is reached (berryman_are_1989, ; nisbet_avoiding_1989, ).

Many-species endogenous fluctuations can only persist if they do not induce too many extinctions. Extinction rates are related to the amplitude of fluctuations giles_leigh_average_1981 ; lande_risks_1993 , their synchrony loreau_biodiversity_2003 and their correlation time ripa_noise_1996 . The peculiarity of endogenous fluctuations is that these properties arise from the species dynamics, and therefore feed back on themselves. A theory of these feedbacks is however lacking.

Here we propose a novel quantitative approach, and show that many-species endogenous fluctuations can persist for extremely long times. Furthermore, they can be realized in experimental conditions, and identified in these experiments by multiple characteristic features. Crucially, we show that states with higher species diversity have stronger fluctuations, and vice versa. We also offer reasons why they may not have been observed in previous studies, and directions in which to search. An important factor in maintaining a dynamically fluctuating state is the spatial extension of the ecosystem, here modeled as a metacommunity: multiple patches (locations in space) that are coupled by migration.

Our strategy is the following. We first propose and simulate experiments to show that persistent fluctuations can be very elusive in a single well-mixed community, yet attainable in a metacommunity via three main ingredients: the existence of multiple patches, moderate migration fluxes coupling them, and differences in conditions between patches. These three ingredients can mitigate the likelihood that large fluctuations within a patch will lead to overall extinctions (see Fig. 1), and make it possible for species to persist in highly fluctuating states.

We then offer a quantitative understanding of this phenomenon. We build on the analytical framework developed in roy_numerical_2019 (dynamical mean-field theory) that allows us to investigate, in a quantitative and predictive way, the conditions under which robust fluctuations can arise from complex interactions. This theory exactly maps a deterministic metacommunity (many-species dynamics over multiple spatial locations) to a stochastic representative metapopulation (single-species dynamics over multiple spatial locations). It predicts the distribution of abundance, survival and variability for a species subjected to “noise” that results from other species in the same community, rather than external perturbations. Dynamical mean-field theory allows us to analyze these fluctuations, and show that the effective stochasticity of species dynamics is a manifestation of high-dimensional chaos.

The intuitive picture that emerges from our analysis is the following: the persistence of endogenous fluctuations, which can be found in a wide range of realistic conditions, requires a balancing act between forces that stabilize and destabilize the dynamics, see Fig. 1. On the one hand, the system needs to preserve a high diversity (both in terms of species number and interaction heterogeneity), as it is known may_will_1972 ; gravel_stability_2016 that lower diversity leads to a stable equilibrium. On the other hand, the system also has to limit excursions towards very low abundances. This requires weeding out species that induce unsustainable fluctuations, and rescuing the others from sudden drops.

To accomplish that, the system relies on asynchronous dynamics between different spatial locations, and finite strength and correlation time of the abundance fluctuations. In addition, despite large fluctuations in the abundances of all species, that strongly affect the species growth in any given patch, long-lasting “sources” emerge for some of the species, i.e. patches where these species are more likely to remain away from extinction. Rare dynamical fluctuations leading to extinction in a given patch are hampered by migration from the other patches, which keeps the system in a non-equilibrium state. We show that, with moderate migration and some spatial heterogeneity, high-diversity dynamical states can be reached where species populations fluctuate over orders of magnitude, yet remain bounded for very long times above their extinction threshold.

Our findings allow us to paint a more precise picture of when persistent endogenous fluctuations can arise. We conclude with a discussion of the implications for biodiversity and ecosystem stability, and predictions for future experiments on community dynamics.

Proposed experiments

In the following, we introduce our results via a set of proposed experiments, realized in simulations, see Fig. 2. These results are later explained in the theoretical analysis. All parameters for simulations are detailed in Appendix A.

We focus on a meta-community which consists of $M$ patches (well-mixed systems) connected by migration, and isolated from the external world. We consider generalized Lotka-Volterra equations for the dynamics of the abundance $N_{i,u}$ of species $i$ in patch $u$:

[TABLE]

where $A_{ij,u}$ are the interactions coupling the species, $B_{i,u}$ represents the equilibrium abundance in absence of interactions and migration (known as the carrying capacity), and $D_{i,uv}$ are the migration rates between patches $u$ and $v$. In addition, an extinction threshold is implemented as follows: when a species’ abundance goes below a cutoff $N_{c}$ in all patches, the species is removed from the metacommunity and cannot return111We are interested in the regime where recolonization by migration between patches is fast compared to the rate of extinction events. In this regime, we expect (and checked in a few cases) that other implementations of the cut-off $N_{c}$ will lead to the same qualitative phenomena. For instance, we implemented patch-wise extinctions when the abundance goes below the threshold in one particular patch, while still allowing migrations in.. This threshold corresponds to the minimum sustainable number of individuals, hence $1/N_{c}$ sets the scale for the absolute population size ($P$) of the species. For simplicity, we take $D_{i,uv}=d/\left(M-1\right)$ and $N_{c}$ identical for all species and patches.

The species are assumed to have unstructured interactions (e.g. they belong to the same trophic level), meaning that $A_{ij,u}$ are sampled independently222In the main text we focus on the asymmetric case in which $A_{ij,u}$ and $A_{ji,u}$ are uncorrelated. We show in Appendix F that our results also hold when correlations are present. and identically for different $\left(i,j\right)$. For a given species pair, its interactions $A_{ij,u}$ vary somewhat with $u$; this variability corresponds to small differences in the conditions between the patches (graham_towards_2018, ). In the simulation examples we set all carrying capacities $B_{i,u}=1$; the phenomena described below are also found if $B_{i,u}$ vary between patches in addition to, or instead of the interaction coefficients.

Our proposed experiments, illustrated by dynamical simulations, are the following:

(A) First, we model a single patch, $M=1$ initially containing $S=250$ species333Each interaction coefficient is non-zero with probability $c=1/8$, and the non-zero interactions are Gaussian with $\operatorname{mean}\left(A_{ij,u}\right)=0.3,\,\operatorname{std}\left(A_{ij,u}\right)=0.45$. . Species go extinct until the system relaxes to a fixed point (stable equilibrium), see left panel of Fig. 2.

(B) We now take $M=8$ patches with the same initial diversity $S=250$ and interaction statistics as in (A). For each pair of interacting species, $A_{ij,u}$ varies slightly with location $u$, with a correlation coefficient $\rho=0.95$ between patches. The abundances now fluctuate without reaching a fixed point, see right panel of Fig. 2. At first the diversity decreases as species go extinct, but this process dramatically slows down, and the diversity is unchanged at times on the order of $10^{5}$, see middle panel of Fig. 2.

Three essential observations emerge from simulating these experiments, and repeating them for different parameters. First, species diversity and the strength of endogenous fluctuations are tightly bound, each contributing to the other’s maintenance. Second, as shown in Fig. 2, species trajectories first go through a transient phase where they fluctuate over many orders of magnitude, causing numerous extinctions which lead to a reduction of variability, until a fixed point (for $M=1$) or non-equilibrium state (for $M=8$) with weaker fluctuations is reached. Third, the qualitative difference between experiments A and B is robust to changes in parameter values. Changes in $N_{c}$ and $d$ affect only quantitatively the states that are reached in experiment B, see Fig. 3(top). For instance, by increasing the population size $P=1/N_{c}$, we can reach dynamically fluctuating states with higher long-time diversities, as shown in Fig. 3(bottom). When the population size is reduced by increasing $N_{c}$, the long-time diversity decreases, but remains high until $N_{c}\sim 10^{-2}-10^{-1}$, where it decreases dramatically. For example, the diversity shown in Fig. 2(right) is $80\%_{\pm 13\%}$ higher than that reached for fixed-points at higher $N_{c}$. Similarly, as long as the migration coefficient is in the range $d\lesssim 0.1$ the main qualitative results remain unaltered.

Theory

We now aim to understand which conditions allow a fluctuating state to be reached and maintained without loss of species.

Dynamical Mean Field Theory

We build on a powerful theory, known as Dynamical Mean Field Theory (DMFT) that exactly maps the deterministic meta-community problem (many species in multiple patches) to a stochastic meta-population problem (single species in multiple patches). When species traits and interactions are disordered, e.g. drawn at random from some probability distributions, all species can be treated as statistically equivalent (barbier_generic_2018, ). We can then describe the whole system by following the trajectory of a single species, randomly sampled from the community, and studying its statistics. In the DMFT framework, the effect of all other species on that single species is encapsulated by an “ecological noise” term generated by their fluctuations. This is analogous to the use, in physics, of thermal noise to represent interactions between an open system and its environment. Since species are statistically equivalent, the properties of this ecological noise can be self-consistently obtained from the dynamics of the single species.

While the theory applies to all times (roy_numerical_2019, ), as discussed in Appendix B, we only consider here the stationary state444We will see in the next section that, when endogenous fluctuations are present, this state is actually metastable, i.e. it is almost stationary on large but finite time-scales. reached after a long time, in which extinctions are already rare. In that state, observables such as the mean abundance are stable over time, and two-time measures, such as correlation functions between times $t$ and $t^{\prime}$, depend only on the difference $t-t^{\prime}$. This entails that each species abundance fluctuates with a finite correlation time, i.e. it tends to return to some constant characteristic value after a finite time.

The result of this mapping is that the abundance $N_{u}$ of a given species in patch $u$ undergoes stochastic dynamics,

[TABLE]

This equation models the dynamics of the target species, including its interactions with other species whose abundances are fluctuating. The contribution of interactions can be separated into a time-independent and a time-dependent parts. The time-independent part goes into $N_{u}^{\ast}$, the characteristic value around which the species abundance will fluctuate in the patch. It differs between species and between patches, due to interactions and to environmental preferences modeled by $B_{i,u}$ in Eq. (1), and follows a multivariate Gaussian distribution. The time-dependent part is encapsulated in $\xi_{u}\left(t\right)$, a Gaussian noise with a finite correlation time.

As the quantities $N_{u}^{\ast}$ and $\xi_{u}\left(t\right)$ result from interaction with other species, which are statistically equivalent to the target species, one can express their properties from the statistics of $N_{u}(t)$ itself, as shown in Appendix B. The most important features are that time-dependent Gaussian noise $\xi_{u}\left(t\right)$ has zero mean and a covariance $C_{\xi}\left(t,t^{\prime}\right)$, which is directly related555It verifies the self consistent equation:

$\left\langle\xi_{u}\left(t\right)\xi_{u}\left(t^{\prime}\right)\right\rangle=\sigma^{2}\left[C_{N,u}\left(t,t^{\prime}\right)-C_{N,u}^{\infty}\right]$

where

$C_{N,u}\left(t,t^{\prime}\right)=\left\langle N_{u}\left(t\right)N_{u}\left(t^{\prime}\right)\right\rangle\,\,,$

$C_{N,u}^{\infty}=\lim_{t-t^{\prime}\rightarrow\infty}\left\langle N_{u}\left(t\right)N_{u}\left(t^{\prime}\right)\right\rangle\ .$

The average $\langle..\rangle$ over the stochastic process corresponds to the average over species in the original Lotka-Volterra equations. to the time auto-correlation $C_{N}\left(t,t^{\prime}\right)$ of $N_{u}\left(t\right)$ within a patch, and to $\sigma^{2}=cS\thinspace\operatorname{var}\left(A_{ij}\right)$ the variance of interactions rescaled as in (may_will_1972, ). Moreover, the covariance of the $N_{u}^{\ast}$ is fixed by the time auto-correlation of $N_{u}\left(t\right)$, both within and in-between patches. In principle, the noise is also correlated between patches, but this is a small effect in the dynamical regime of interest to us, see the next section. Note that $C_{\xi}\left(t,t^{\prime}\right)$ vanishes when a stable equilibrium is reached.

The analysis of the DMFT equations clarifies the main effect of coupling patches by migration: patches with higher $N_{u}^{\ast}$ tend to act as sources, i.e. the species most often grows there, and migrates out to sites where it cannot grow (sinks). We show directly from simulations of the Lotka-Volterra equations in Fig. 4 that species have particular patches which tend to act as sources consistently over long times. This fact is counter-intuitive, as the abundances of all species may be fluctuating over orders of magnitude in any given patch, yet this patch will retain its identity as a source (or sink) when averaging over long time periods. The variability of the $N_{u}^{\ast}$s between patches thus leads to an insurance effect, since it is enough to have one patch acting as a source to avoid extinction of the species in the others.

Reaching and maintaining a dynamical state

Let us first consider a single community ($M=1$). For a species to survive for long periods of time, it follows from DMFT that it must have positive $N^{\ast}$, or else $N\left(t\right)$ decays exponentially until the species goes extinct. Even if $N^{\ast}>0$, there is still a probability (per unit time) of extinction, which depends on $N^{*},\,N_{c}$ and on the strength of the noise $\xi\left(t\right)$. Following extinctions, a remaining species interacts with fewer fluctuating other species, causing the strength of the noise to decrease, and with it the probability for extinction, see feedback loop in Fig. 5.

We can develop an analytical treatment for very small cut-off $N_{c}$ (large population size). In this case there is a large difference in time-scales between the short-term dynamics induced by endogenous fluctuations, and the long-term noise-diversity feedback cycle discussed above. In fact, the extinctions driving this feedback are due to rare events in which the abundance of species with a positive $N^{*}$ decreases below the (very small) cut-off $N_{c}$. For a species in an isolated patch ($M=1$), the time-scale for such an event is known giles_leigh_average_1981 ; lande_risks_1993 to be of order of $\tau\left(1/N_{c}\right)^{a}$ where $\tau$ is a characteristic time of the endogenous fluctuations, and $a=2N^{*}/W$ is independent666The expresson of the time scale is analogous to the Arrhenius law for activated processes in physics and chemistry: in this case, the counterpart of the energy barrier is $-[N^{*}\ln N_{c}]$ and fluctuation amplitude $W$ plays the role of the temperature. of $N_{c}$, with $W$ the amplitude of the endogenous fluctuations,

[TABLE]

The important point here is that, although endogenous fluctuations disappear eventually, there is a clear separation of time-scales between typical endogenous fluctuations, that are fast and lead to a quasi-stationary dynamical state, and rare extreme fluctuations that cause extinctions and push the ecosystem into a different state.

While a single community might in principle achieve long-lasting endogenous fluctuations, this however requires unrealistically large population sizes and species number, see Appendix E. Migration between multiple patches substantially enhances persistence due to the spatial insurance effect (loreau_biodiversity_2003, ): species are more unlikely to go extinct because they need to disappear everywhere at once. The time scale for such an event is $\tau\left(1/N_{c}\right)^{Ma_{\mathrm{eff}}}$ with $a_{\mathrm{eff}}=2N_{\mathrm{eff}}^{\ast}/W$ where $N_{\mathrm{eff}}^{\ast}$ is an effective value for $N^{*}$ of a species across patches, an expression for which is given in Appendix C. This result is thus similar to the one identified above for one patch, raised to the power $M$. These results assume that $W$ is finite, and that the noise acting on a species is independent between patches (asynchrony). Indeed, for moderate values of $D$ and $\rho$ not too close to one, simulations show that $C_{\xi}\left(t,t^{\prime}\right)$ is a well-behaved function of $t$ so that $W$ is finite, and the correlation between patches is found to be very small, see Appendices B,C.

These expressions provide a quantitative description of the feedback cycle in Fig. 5. Endogenous fluctuations disappear on the time scale at which species with characteristic abundance $N_{\mathrm{eff}}^{\ast}$ of order one would go extinct. We must further account for the vanishing strength of the noise $W$ as species disappear. This is shown in Fig. 5(a), where the strength of the fluctuations is tightly linked to species diversity, and is zero at the diversity of fixed points. Hence, extinctions significantly increase $a_{\mathrm{eff}}$, reducing the chance for further extinctions.

This picture agrees well with the analytical predictions, which can be obtained for very small but positive $N_{c}$ and $D$ by DMFT. These give the distribution of $N_{\mathrm{eff}}^{\ast}$, whose integral over positive $N_{\mathrm{eff}}^{\ast}$ amounts to the maximal total diversity at long times. We show in Fig. 6 that extinct species are generally those that have lower values of $N_{\mathrm{eff}}^{\ast}$: we compare the analytically predicted distribution of $N_{\mathrm{eff}}^{\ast}$ to the one observed in simulations, and see that most missing species had low values of $N_{\mathrm{eff}}^{\ast}$. The difference becomes smaller for lower $N_{c}$. Due to these differences, the obtained diversities are lower than the theoretical maximum, see Fig. 3(top).

As stressed above, the asynchrony of fluctuations in different patches is crucial: it allows some species to survive with positive characteristic abundance $N_{i}^{*}$ in at least one of the patches. This leads to a higher total number of long-term persisting species, decreases the likelihood of fluctuations to small abundances, and hence increases the stability of a dynamically fluctuating state. If the migration rate $D$ is too strong or $\rho$ very close to one, dynamics in the different patches synchronize, quickly annulling the insurance effect. However, minor (few percent) changes in interaction coefficients or carrying capacities between patches are enough to maintain this effect, see Fig. 5(bottom); we don’t need to impose coexistence artificially, e.g. by requiring that every species has at least one refuge (a patch so favorable to it that it always dominates there). These little variations in the interaction coefficients are highly plausible, as interaction strength can vary with many factors, including resource availability (carstensen_beta_2014, ), or temperature and its influence on metabolism (kordas_community_2011, ). The heterogeneity $\rho$ required to reach a fluctuating state decreases with $M$, see Fig. 5(bottom).

In practice, maintaining a dynamical state seems unfeasible for only one patch, at least for reasonable values of population size $P=1/N_{c}$ and species number $S$, see discussion in Appendix E. Yet the combined effect of the two phenomena described above allows for very long-lived endogenous fluctuations in metacommunities, already for $M=2$ patches.

Conclusions

Species interactions can give rise to long-lasting fluctuating states, which both require and allow the maintenance of high species diversity. This can happen under a wide range of conditions, which we have illustrated in simulated experiments, and identified through an analytical treatment based on Dynamical Mean-Field Theory.

While we have drawn parallels with the theory of stability and coexistence in externally-perturbed ecosystems loreau_biodiversity_2003 ; giles_leigh_average_1981 ; lande_risks_1993 ; allen_chaos_1993 , our approach also highlights essential differences between environmentally-driven and endogenous fluctuations. We show that many-species dynamics induce feedback loops between perturbation and response, and in particular a tight relationship between fluctuation strength and species diversity, which are absent from externally-perturbed ecosystems. Moreover, while similar species can display correlated responses to environmental stochasticity loreau_biodiversity_2013 , we expect here that their trajectories will be starkly different and unpredictable, due to high-dimensional interactions which lead to complex dynamics. The resulting picture from DMFT is that the abundance of any given species undergoes stochastic dynamics with a finite correlation time. This means that the trajectory of the species abundance cannot be predicted after a time that is large compared to the correlation time–a hallmark of chaos, also found in other models of high-dimensional systems sompolinsky_chaos_1988 . Our theory paves the way for quantitative testing of these fingerprints of diversity-driven fluctuations in data.

In a counterpoint to classic results (may_will_1972, ), we have shown that, while highly diverse ecosystems are unstable, they might still persist: extinctions can be avoided and biodiversity maintained, despite species abundances fluctuating over multiple orders of magnitude. We do observe a negative feedback loop, in which endogenous fluctuations cause extinctions, and eventually lead to their own disappearance as the ecosystem reaches a lower-diversity stable equilibrium. But this self-suppression of fluctuations can be mitigated by a number of factors, among which space is particularly important.

In a single well-mixed community, we expect that persistent fluctuations might not be observed in practice: while theoretically possible, they may require unrealistic population sizes and species numbers. But spatial extension and heterogeneity can dramatically reduce these requirements, in a way that parallels the insurance effect against exogenous perturbations. When fluctuations are not synchronized across space, some patches can act as sources, from which failing populations will be rescued through migration (allen_chaos_1993, ; loreau_biodiversity_2003, ). Here, we find that the existence of sources is surprisingly robust: even if there is no location where the environment is favorable to a given species, source patches can arise from interactions, and endure for long times despite the large fluctuations in species abundances. By allowing fluctuations without extinctions, spatial heterogeneity helps maintain species diversity, and thus the fluctuations themselves. This result is robust over a wide range of parameters, as it only calls for moderate values of inter-patch migration: the rate $D$ must be such that, over the typical time scale of abundance fluctuations, many individuals can migrate out of a patch (allowing recolonization in the absence of global extinction), while representing only a small fraction of the population in that patch.

A crucial result is that this condition suffices to ensure that synchronization between patches is absent, and that the total strength and correlation time of the noise within patches ($W$ above) remain bounded for finite populations and finite migration rates between patches. This is in contrast to alternative scenarios where noise correlations decay slowly with time halley_ecology_1996 . This result is non-trivial for endogenous fluctuations, as the existence of feedbacks (encoded in the self-consistent equations of the DMFT framework) can potentially lead to synchronization and long-time correlations in the noise. Yet we demonstrate that synchrony is avoided, both through direct simulations, and by building an analytical theory based on these assumptions, whose predictions match simulations quantitatively.

In conclusion, non-equilibrium fluctuating states might be much more common than suggested by experiments and theory for well-mixed communities. And since these fluctuations permit the persistence of more species than could coexist at equilibrium, we might also expect significantly higher biodiversity in natural environments.

Acknowledgments - It is a pleasure to thank J.-F. Arnoldi, J.-P. Bouchaud, C. Cammarota, and M. Loreau for helpful discussions. We thank D. S. Fisher for valuable discussions and sharing the results of his parallel work dsf , and in particular for key inputs concerning the instability of chaos without migration and dynamical fluctuations in the chaotic phase in the limit of small migration from the mainland. G. Bunin acknowledges support by the Israel Science Foundation (ISF) Grant no. 773/18. M. Barbier was supported by the TULIP Laboratory of Excellence (ANR-10-LABX-41) and by the BIOSTASES Advanced Grant, funded by the European Research Council under the European Union’s Horizon 2020 research and innovation programme (666971). F. Roy acknowledges support by Capital Fund Management - Fondation pour la Recherche. G. Biroli was partially supported by the Simons Foundation collaboration Cracking the Glass Problem (No. 454935).

Appendix A Model parameters used in simulations

For convenient reference, this Appendix includes the parameters for all simulations. The model is given in Eq. (1). All $B_{i,u}=1$ and all $D_{i,uv}=d/\left(M-1\right)$. The $A_{ij,u}$ are independent for different $\left(i,j\right)$ pairs (except in Appendix 13).

In Fig. 2, the probability of $A_{ij,u}$ to be non-zero is $c=1/8$, and the non-zero elements are sampled from a normal distribution with $\operatorname{mean}\left(A_{ij,u}\right)=0.3,\,\operatorname{std}\left(A_{ij,u}\right)=0.45$. The same elements $A_{ij,u}$ are non-zero across all patches $u$. The correlation coefficient between non-zero $A_{ij,u}$ in different patches is $\rho=\operatorname{corr}\left[A_{ij,u},\,A_{ij,v}\right]=0.95$ for $u\neq v$. (The correlation is $0.964$ when interactions with $A_{ij,u}=0$ are also counted.) The initial (pool) diversity is $S=250$. In Fig. 2(A), $M=1$. In Fig. 2(B), $M=8$ patches and $d=10^{-3}$. The cutoff is $N_{c}=10^{-15}$.

Fig. 3(bottom), uses the same parameters as Fig. 2, but with a range of values for $d,S$ and $N_{c}$.

Fig. 4 uses the runs shown in Fig. 2(B). Standard deviation and mean are estimated from 1601 time points during the time period $t=[10^{4},2\cdot 10^{5}]$. Fig. 6 uses multiple runs, with the same parameters as 2, except for $d=10^{-4}$ and the values of $N_{c}$ that are detailed in the figure legend.

Fig. 5(top) uses the same parameters as 6. Fig. 5(bottom), shows the line where half of the runs are fixed points, and half continue to fluctuate until $t=2\cdot 10^{5}$. It uses same parameters as 6, except with $D=d/\left(M-1\right)=10^{-4}$.

In Fig. 5(top), the size of the fluctuations are calculated from $\operatorname{var}\left(\xi_{u}\right)=\left\langle\xi_{u}^{2}\left(t\right)\right\rangle=\sigma^{2}C_{N,u}\left(t,t\right)$, with $C_{N,u}\left(t,t\right)=\left\langle N_{u}^{2}\left(t\right)\right\rangle-\lim_{t-t^{\prime}\rightarrow\infty}\left\langle N_{u}\left(t\right)N_{u}\left(t^{\prime}\right)\right\rangle$. For more details on the averaging, see Appendix B, Fig. 7.

Appendix B DMFT equations

In this section, we present the full DMFT equations, and explain how they can be reduced to the steady-state equations quoted in the main text.

We consider as a starting point equation Eq. (1). For the sake of clarity, we derive DMFT under simplifying assumptions, but the result is much more robust and could be applied to different ecology models as well as real data [17]. DMFT for ecological models has a double valency analogous to the one of mean-field theories in physics: it is at the same time an exact theory for some simple models, and a powerful approximation largely applicable to a broad range of systems. For the sake of clarity, the derivation assumes a fully connected model (all interactions are non-zero), but the results hold for any connectivity $C$ as long as $C\gg 1$, see remark at the end of this Appendix.

The assumptions which make DMFT exact are the following: all constants $N_{i,u}(0)$, $B_{i,u}$, $D_{i,uv}$ and $A_{ij,u}$ are random variables, sampled from known distributions. More precisely:

• •

in each patch $u$ and for all species $i$, the parameters $X^{u}=\{N_{i,u}(0),\,B_{i,u},\,D_{i,uv}\}_{i=1}^{S}$ are drawn from a probability distribution $\mathbb{P}$ which is a product measure $\mathbb{P}_{u}(X^{u})=\prod_{i=1}^{S}\mathbb{P}(X_{i}^{u})$;

• •

the interaction matrix can be decomposed as $A_{ij,u}=\mu/S+\sigma/\sqrt{S}\;a_{ij,u}$. $a_{ij,u}$ are standard random variables with mean zero, variance one, and correlation:

[TABLE]

where we used the Kronecker symbol $\delta_{ik}$, and $\rho_{uv}=\rho+(1-\rho)\delta_{uv}$ is a uniform correlation $\rho$ between patches.

With these conventions, we rewrite Eq. (1) in the following way:

[TABLE]

where $m_{u}(t)=S^{-1}\,\sum_{i=1}^{S}N_{i,u}(t)$ is the mean abundance in patch $u$, and $\eta_{i,u}(t)=-\sigma S^{-1/2}\,\sum_{j=1}^{S}a_{ij,u}\,N_{j,u}(t)$ is a species-and-patch-dependent noise.

The DMFT equation can be obtained by following Ref. [14]: in the large-$S$ limit, it can be shown that the statistics of this multi-species deterministic process corresponds to the following one-species stochastic process, for each patch.

[TABLE]

where $\{N_{u}(0),\,B_{u},\,D_{uv}\}$ are sampled from the distribution $\mathbb{P}(X^{u})$, $m_{u}(t)$ is a deterministic function, and $\eta_{u}(t)$ is a zero-mean Gaussian noise. The variability from one species to another becomes in the DMFT setting the randomness contained in $\{N_{u}(0),\,B_{u},\,D_{uv}\}$ and $\eta_{u}(t)$.

To make this point crystal clear, let us introduce two different averages:

• •

$\overline{Y}$ averages over the stochastic process in Eq. (2): over the stochastic noise $\eta_{u}$ and over the distribution $\mathbb{P}(X^{u})$;

• •

$\mathbb{E}_{S}(Y)$ denotes the statistical average over the deterministic multi-species system. $\mathbb{E}_{S}(Y)=\sum_{i=1}^{S}Y_{i}$, and therefore also includes sampling of $X_{i}^{u}$.

DMFT represents in terms of a stochastic process the deterministic dynamical system governing the dynamics of the $S$ species in the ecosystem. In consequence, averages over the stochastic process coincide with average over species: for a given observable $Y$: $\overline{Y}=\lim_{S\to\infty}\mathbb{E}_{S}(Y)$. This is analogous to the representation of the environment of an open physical system in terms of thermal noise, as it is done e.g. in the case of the Langevin equation.

The second important aspect of DMFT is self-consistency. This is related to the fact that the noise is induced by the dynamics of the species themselves, so its properties can be obtained from dynamical averages:

[TABLE]

where we used a last average $\langle\cdot\rangle$ over the stochastic noise only, in order to define its covariance. Henceforth we use the notation $C_{uv}^{N}(t,t^{\prime})=\rho_{uv}\,\overline{N_{u}(t)\,N_{v}(t^{\prime})}$.

We now show how DMFT equations simplify for a time-translationally-invariant state of the system, which is in general reached after some transient time. In this state, all one-time observables become constant in time, and two-time observables become functions of the time difference only.

[TABLE]

The correlation $C_{uv}^{N}(t-t^{\prime})$ decays at large time differences to a non-zero constant, leading to a static contribution to the noise term. In order to disentangle the static part and the time-fluctuating part of the noise, we perform the decomposition $\eta_{u}(t)=z_{u}+\xi_{u}(t)$ such that $z_{u}$ and $\xi_{u}(t)$ are independent zero-mean Gaussian variables and processes verifying:

[TABLE]

and subsequently $\langle\xi_{u}(t)\xi_{v}(t^{\prime})\rangle=\sigma^{2}\,C_{uv}^{N}(t-t^{\prime})-\sigma^{2}\,\lim_{t-t^{\prime}\rightarrow\infty}C_{uv}^{N}(t-t^{\prime})$ which vanishes for $t-t^{\prime}\rightarrow\infty$.

Substituting this decomposition into Eq. (2), we obtain:

[TABLE]

where $N_{u}^{*}=1-\mu\,m_{u}+z_{u}$ is a Gaussian variable, whose statistics is described in Appendix D. We checked numerically that for small migration $D$, the noise is only correlated between patches through its static part: for $u\neq v$, $\xi_{u}(t)\,\xi_{v}(t^{\prime})\ll z_{u}\,z_{v}$, as presented in Fig. 7. In this case, we can write the self-consistent closure as follows:

[TABLE]

As explained above, DMFT can be implemented as an approximation for a large variety of systems. In this case one has to infer the average $\mu$, the standard deviation $\sigma$ of interactions, and the distribution $\mathbb{P}(X^{u})$ from the data (we remind that $X^{u}=\{N_{u}(0),\,B_{u},\,D_{uv}\}$) and use them as an input to define an effective model. The generalization to patch-dependent cumulants $\mu_{u}$ and $\sigma_{u}$ is quite straightforward. So is the generalization to patch-dependent correlation $\rho_{uv}$.

We have derived DMFT for a completely connected set of interactions $A_{ij}$. A different way to obtain DMFT is considering a finite connectivity network of interactions $A_{ij}$, e.g. the one produced by a Erdos-Renyi random graph with average connectivity per species $C$ or a regular random graph with connectivity $C$. In these cases, for each link $ij$ one generates a random variable with average $\mu/C$ and variance $\sigma^{2}/C$ and set it to $A_{ij}$. In the large connectivity limit, $C\rightarrow\infty$, each species interacts with a very large number of species and one can replace the deterministic interaction with an effective stochastic noise, as done for a completely connected lattice. Although the resulting DMFT equations are the same, the two cases are quite different: in the former a species interact with $C\ll S$ species whereas in the latter a species interacts with $C=S$ species. The equivalence of DMFT for completely connected lattices and finite connectivity ones in the $C\rightarrow\infty$ limit has been thoroughly studied in physics of disordered systems in the last twenty years [23].

Appendix C Extinction probability of a species

Here the probability of extinction of a species is presented, at the limit $N_{c}\ll D\ll 1$. More specifically, we assume that $N_{c}$ is small compared to the typical fluctuations of the abundances. In addition, in simulations we see that it is reasonable to assume complete lack of synchrony, namely that the noise $\xi_{u}$ is uncorrelated between different patches, see Appendix B, Fig. 7. We will therefore assume that in the following calculation. Finally, we assume that for at least one patch, $N_{u}^{*}>0$, otherwise the species goes quickly extinct.

Within DMFT, the problem thus becomes ones of calculating the extinction probability of a meta-population (single species), under environmental fluctuations, that are uncorrelated between the different patches. We only present the result here; a full account will be given elsewhere.

Let $x_{u}\equiv\ln N_{u}$. The equations of the DMFT, Eq. (2), become

[TABLE]

Here $D=d/\left(M-1\right)$. We look for a rare realization of $\left\{\xi_{u}\right\}$ that makes all the $x_{u}$ reach $x_{c}=\ln N_{c}$, in the case when the cut-off is low, $x_{c}\rightarrow-\infty$. The calculation proceeds within the framework of large-deviation theory [24]. First, one defines the “action”

[TABLE]

with $\xi_{u}$ substituted with its value from Eq. (4), and $W$ defined as in the main text. Here we approximated the noise correlations by white noise, which is justified here as the extinction event takes a time which is long compared to the correlation time. We assume that $D$ is small.

Then the mean time to the occurrence of such an event scales as $P\sim e^{J_{min}}$ with $J_{min}$ the action $J$ minimized over all population trajectories $\left\{x_{u}\left(t\right)\right\}_{u=1..M}$ that start at $t\rightarrow-\infty$ at the typical value of $x_{u}$, obtained by the zero-noise fixed point of Eq. (4), and terminate at $t_{f}$ at $x_{c}=\ln N_{c}$ .

We first describe the result for $M=1$. In this case there is only one patch, $u=1$, with $N_{1}^{\ast}$. If $N_{1}^{\ast}<0$ the species is extinct. On the other hand, if $N_{1}^{\ast}>0$, then we obtain the known result [11, 25]

[TABLE]

The result for all $M$ is a generalization of this result, of the form

[TABLE]

To describe the calculation of $N_{\mathrm{eff}}^{\ast}$, order the patches so that $N_{1}^{\ast}\geq N_{2}^{\ast}\geq..\geq N_{M}^{*}$. Then there exists $1\leq m\leq M$ such that

[TABLE]

and where $w$ satisfies: $w\leq N_{u}^{\ast}$ for all $u\leq m$, and $w>N_{u}^{\ast}$ for all $u>m$. Such a partition can be shown to always exist. Then

[TABLE]

The derivation will be given elsewhere. We illustrate the result by considering two cases. First, in the $M=1$ example, since $N_{1}^{\ast}>0$, the partition is $\left\{u\right\}_{+}=\left\{1\right\}$ and $\left\{u\right\}_{-}$ the empty set. Indeed, this gives $w=-N_{1}^{\ast}$, so $w\leq N_{1}^{*}$. Then $N_{\mathrm{eff}}^{\ast}=N_{1}^{\ast}$ and $J_{\min}=\frac{2x_{c}}{W}N_{1}^{\ast}$, so the result for $M=1$ is reproduced. Another simple case is when there are $M$ patches with identical carrying capacities $N_{u=1..M}^{\ast}=N^{\ast}$. Here $\left\{u\right\}_{+}=\left\{1,..,M\right\}$ and $w=-N^{\ast}$. Then $N_{\mathrm{eff}}^{\ast}=N^{\ast}$, and $J_{\min}=\frac{2x_{c}}{W}MN_{\mathrm{eff}}^{\ast}=\frac{2x_{c}}{W}MN^{\ast}$. This result is intuitively clear: to go extinct, the species must go extinct in all patches at once, so the probability is $P\sim\exp\left(\frac{x_{c}}{W}MN^{\ast}\right)\sim\left(P_{1}\right)^{M}$, where $P_{1}$ is the result for $M=1$.

Appendix D Diversity and stability at low migration rates

We use notations from Appendix B. Within the time-translational-invariant state:

[TABLE]

Consider the case of low migration, $D\rightarrow 0^{+}$. We now develop a theory assuming that the amplitude of the endogenous fluctuations,

[TABLE]

remains finite in the limit $D\rightarrow 0^{+}$. Assume the species survives, i.e. there is at least one patch with $N_{u}^{*}>0$. If $N_{u}^{*}<0$ then $N_{u}=O\left(D\right)$. If $N_{u}^{*}>0$ then $N_{u}=O\left(1\right)$ and therefore $\sum_{v\sim u}D_{uv}\left(\frac{N_{v}}{N_{u}}-1\right)=O\left(D\right)$. Taking the time average of the above equation

[TABLE]

and therefore $\overline{N_{u}}=N_{u}^{*}+O\left(D\right)$.

The previous arguments lead to the conclusion that in the $D\rightarrow 0^{+}$ limit $\overline{N_{u}}=N_{u}^{*}$ if $N_{u}^{*}>0$ and is equal to zero otherwise. In the following we provide more detail more this argument and its possible limitations. For this last equality to be valid, we need that $\sum_{v\sim u}\left(\frac{N_{v}}{N_{u}}-1\right)$ will be finite, so that $D\sum_{v\sim u}\left(\frac{N_{v}}{N_{u}}-1\right)$ will indeed be small. This might break if $N_{u}$ can be small while some other $N_{v}$ remains $O\left(1\right)$. An estimate for that proceeds by noting that the carrying capacity of patch $u$ in the presence of other patches is larger or equal to $N_{u}^{*}-MD\simeq N_{u}^{*}$, its carrying capacity alone. If patch $u$ fluctuates alone, then

[TABLE]

This gives for $\overline{1/N_{u}}$

[TABLE]

For any given $N_{u}^{*}$ this is finite. It diverges as $N_{u}^{*}\rightarrow 0$. Therefore the migration term is negligible only if $\frac{DW}{1-D}\simeq DW\ll N_{u}^{*}$. (Note that migration itself would limit $N_{u}$ going below much below $DN_{v}$, which would make this term smaller.) The main approximation (or limitation) of our approach is the assumption that $W$ remains finite in the small $D$ limit. This is shown to hold in simulations presented in Appendix B. It breaks down if the noise develops long-lasting correlations in time. Our approximation will be nevertheless good for large $|N_{u}^{*}|$ and for weak endogenous fluctuations.

We now used the relationship discussed above between $\overline{N_{u}}$ and $N_{u}^{*}$ to determine the statistics of $N_{u}^{*}$. We shall use the term “source” for patches where $N_{u}^{*}>0$, and “sink” otherwise777The term “source” is used here so as to include patches (sometimes referred to as pseudo-sinks) where a species might still receive migration from patches with even larger $N_{u}^{\ast}$. But the contribution of this migration is small and not required for its persistence.. In order to understand the correlation between the sources in the different communities, we unpack $N_{u}^{*}$ using Appendix B. Taking the time-average is equivalent to averaging over the dynamical noise $\xi$. Therefore, in patch $u$ for species $i$, $z_{i,u}=-\sigma S^{-1/2}\sum_{j}a_{ij,u}\,\overline{N_{j,u}}=-\sigma S^{-1/2}\sum_{j,+}a_{ij,u}\,N_{j,u}^{*}$. The sum $\sum_{j,+}$ means that we only sum over $N_{j,u}^{*}>0$. Here, we recall that $a_{ij,u}$ are standard random variables with mean zero, variance one, and correlation between patches:

[TABLE]

where we used the Kronecker symbol $\delta_{ik}$.

Therefore:

[TABLE]

where we recall $m_{u}=\left\langle\overline{N_{i,u}}\right\rangle=\left\langle N_{i,u}^{*}\right\rangle_{+}$. We can now compute the different moments of the multivariate Gaussian random variable $N_{u}^{*}$, using equation (6). We obtain the closure:

[TABLE]

When $u=v$, as $\rho_{uu}=1$, we find the expected single community result. In particular, $\operatorname{mean}\left[N_{u}^{*}\right]$ and $\mathrm{variance}\left[N_{u}^{*}\right]$ do not depend on the patch $u$.

We numerically solve the closure in a self-consistent way: start with a guess for $\left\langle N_{i,u}^{*}\,N_{i,v}^{*}\right\rangle_{+}$, and then (1) Produce many samples of the vector $N_{u=1..M}^{*}$ and (2) calculate the next estimate for $\left\langle N_{i,u}^{*}\,N_{i,v}^{*}\right\rangle_{+}$, by averaging only over $N_{i,u}^{*}$ and $N_{i,v}^{*}$ that are both positive. For stability of this numerical scheme, we only replace half the samples at each iteration. We use $10^{5}$ samples and 1000 iterations. The algorithm is always found to converge to the same solution.

Given $\mathrm{covariance}\left[N_{u}^{*},N_{v}^{*}\right]$, the distribution of $N_{i,u}^{*}$ is completely specified: it is a multivariate Gaussian in $u$, has the single-patch statistics of a single community, and a known covariance between patches. The solution can then also give the distribution of the number of sourcing patches.

In addition, we can compute the correlation coefficient $\rho_{N^{*}}$ of the $N_{u}^{*}$’s. We use here our simple case of a uniform correlation $\rho_{a}$ between patches $\rho_{uv}=\rho_{a}+(1-\rho_{a})\delta_{uv}$. We introduce the notation $\rho_{a}$ instead of ‘$\rho$’ in this section in order to avoid confusion with $\rho_{N^{*}}$.

[TABLE]

The results are surprising: even when $\rho_{a}\rightarrow 1$, the overlap between communities is not perfect ($\rho_{N^{*}}<1$), so the total diversity is larger than the one in each patch. This happens exactly at the transition to chaos at $\sigma_{c}=\sqrt{2}$, see Fig. 8.

On Fig. 10, we compare the theory predictions to simulations. In terms of diversity, the theory appears to give an upper bound to the simulations. The difference becomes larger at higher values of $\sigma$, and for $\rho_{a}$ closer to one. To look further into this difference, it is useful to study diversity as a function of the value of $N_{\mathrm{eff}}^{*}$ of each species. As shown in Fig. 6 in the main text, most of the difference in diversity is due to low values of $N_{\mathrm{eff}}^{*}$, which are precisely the species that are more likely to go extinct, with good agreement with theory at higher values of $N_{\mathrm{eff}}^{*}$. This is demonstrated in Fig. 9, which shows that the theoretical prediction for the number of species with $N_{\mathrm{eff}}^{*}>0.2$ is closer to simulation results than the predictions for total diversity. At the moment we do not know if remaining differences are because the theoretical argument is only approximate, or whether in principle, with exceedingly low values of $N_{c}$ and $D$, it could be approached by simulations for any $\sigma$.

To find the boundary of parameter space where fixed points loose their stability and the system becomes chaotic, we look at the linear stability of persistent species. When $D$ is small, the species that are not sourced in each patch do not affect the stability, and so the question simplifies to single patch stability, which when $\operatorname{corr}\left[A_{ij},A_{ji}\right]=0$, results in $\sigma_{c}=\sqrt{2}$ and with $1/2$ of the species being sourced in each patch [26].

Appendix E Single patch ($M=1$)

Here we show that in principle a single patch can reach and maintain a dynamically fluctuating state. However, this requires prohibitively large $S$, not attainable in practice. In Fig. 11 and Fig. 12 we show results of a numerical solution [14] to the DMFT equations detailed in Appendix B. At extremely low values of $N_{c}$ the system appears to reach a final diversity above the May bound and, hence, to be chaotic. DMFT however describes the behavior in the $S\gg 1$ limit. When full simulations of the model in Eq. (1) are carried out at finite $S$, they diversity falls somewhat below the DMFT final diversity, leading to a fixed point, rather than a chaotic state, see Fig. 12. This finite-size correction to the DMFT result are important since they show that maintaining a dynamically fluctuating state for realistic values of $S$ is not possible for $M=1$.

Appendix F Correlations of interactions in a pair of species

In the main text we assumed that $A_{ij,u}$ is sampled independently from $A_{ji,u}$. Here we show that the long-lived endogenous fluctuations can be found even if this assumption is relaxed. For this purpose, we consider a symmetric network of non-zero $A_{ij,u}$, namely $A_{ij,u}\neq 0$ if and only if $A_{ji,u}$. We define $\gamma$ the correlation of the non-zero elements $\gamma=\operatorname{corr}\left[A_{ij,u},\,A_{ji,u}\right]_{A_{ij,u}\neq 0}$. Fig. 13 shows two simulations, one with $\gamma>0$ and the other with $\gamma<0$. In both cases the system relaxes to a long-lived state with fluctuating abundances, without further loss of diversity up to time $2\cdot 10^{5}$. They are intended solely to demonstrate that conditions with $\gamma\neq 0$ exist, rather than a systematic exploration of such cases.

The parameters for the simulations (using the notation of Appendix A) are the following:

Run with positive $\gamma$: $\gamma=1/4$, $S=350$, $\operatorname{mean}\left(A_{ij,u}\right)=0.075$, $\operatorname{std}\left(A_{ij,u}\right)=0.175$, $c=0.357$, $M=8$, $d=10^{-3}$, $\rho=0$, $N_{c}=10^{-15}$.

Run with negative $\gamma$: $\gamma=-1/2$, $S=250$, $\operatorname{mean}\left(A_{ij,u}\right)=0.075$, $\operatorname{std}\left(A_{ij,u}\right)=0.358$, $c=0.5$, $M=8$, $d=10^{-3}$, $\rho=0$, $N_{c}=10^{-15}$.