Complex interactions can create persistent fluctuations in high-diversity ecosystems.

When can ecological interactions drive an entire ecosystem into a persistent non-equilibrium state, where many species populations fluctuate without going to extinction? We show that high-diversity spatially heterogeneous systems 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 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.

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/ (M − 1). The A ij,u are independent for different (i, j) pairs (except in Appendix G).
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 mean (A ij,u ) = 0.3, std (A ij,u ) = 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 ρ = corr [A ij,u , A ij,v ] = 0.95 for u = 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 For each i, j, the M values A ij,u=1..M are drawn simultaneously from a multi-variate normal distribution with correlation matrix C uv = ρ + (1 − ρ) δ uv , using standard numerical methods (e.g., as implemented in the Matlab function mvnrnd). In Fig. 5(a), the size of the fluctuations are calculated from var (ξ u ) = ξ 2 u (t) = σ 2 C N,u (t, t), with C N,u (t, t) = N 2 u (t) − lim t−t →∞ N u (t) N u (t ) . For more details on the averaging, see Appendix B, Fig. A. Fig. 5(a) shows the strength of noise at different diversities. Extinctions beyond the time shown in simulations (t = 2 · 10 5 ) take extremely long times to happen, so reaching these lower diversities in simulations is unfeasible. Instead, we remove species that are most likely to go extinct. Recalling that the time to extinction is τ (1/N c ) 2M N * eff /W , we remove species with the lowest N eff . This is done by running the system for time ∆t = 15·10 3 , calculating N eff as for Fig. 6, and removing the 5 species with the lowest values of N eff . This process is repeated. Other protocols for species removal where attempted, such as increasing N c in time; they give similar results. The results are averaged over 3 runs, with independent sampling of interactions and initial conditions. Fig. 7 shows the linear stability λ stab of fixed points. To calculate this, the equations of motion, Eq. (1), are linearized around a fixed point N * i,u . This gives the matrix equation d − → n /dt = G − → n , where − → n is a vector with one entry for each patch of a surviving species, with G is a square matrix, which depends on the equilibrium abundances N * i,u . For a state with fluctuating abundances, we use the time average of N i,u (t) for N * i,u . Let {λ j } be all the eigenvalues of G. Then λ stab = max j (λ j ), where (..) denotes the real part. Note that a fixed point is stable if λ stab < 0. The results are averaged over the same 3 runs as in Fig. 5(a).

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 [1]. 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 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 • The interaction matrix can be decomposed as A ij,u = µ/S + σ/ √ S a ij,u . a ij,u are standard random variables with mean zero, variance one, and correlation: where we used the Kronecker symbol δ ik , and ρ uv = ρ + (1 − ρ)δ uv is a uniform correlation ρ between patches.
With these conventions, we rewrite Eq. (1) in the following way: is the mean abundance in patch u, and η i,u (t) = −σS −1/2 S j=1 a ij,u N j,u (t) is a species-and-patch-dependent noise.
The DMFT equation can be obtained by following Ref. [2]: 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.
is a deterministic function, and η 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 η u (t).
To make this point crystal clear, let us introduce two different averages: • Y averages over the stochastic process in Eq. (B1): over the stochastic noise η u and over the distribution P(X u ); • E S (Y ) denotes the statistical average over the deterministic multi-species system. E S (Y ) = S i=1 Y i , and therefore also includes sampling of X u i . 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 [2][3][4]: for a given observable Y : Y = lim S→∞ 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 selfconsistency. 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: where we used a last average · over the stochastic noise only, in order to define its covariance. Henceforth we use the notation C N uv (t, t ) = ρ uv N u (t) N v (t ). These relations exactly take into account the correlations that emerge between the abundances and the interactions.
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.
The correlation C N uv (t − t ) 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 η u (t) = z u + ξ u (t) such that z u and ξ u (t) are independent zero-mean Gaussian variables and processes verifying: where N * u = 1 − µ 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 In this case, we can write the self-consistent closure as follows: 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 µ, the standard deviation σ of interactions, and the distribution 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 µ u and σ u is Figure A. Covariance of the abundances in distinct patches. We use the general notation cov(Yu, . Left: In full lines we show the abundance covariance within a patch u = v, and across patches u = v in dotted lines. The correlation in abundances across patches is mainly static: dotted lines are reasonably flat. In other words, the correlation of ξu with ξv for u = v is very small. Right: The covariance in ξ is shown to reach a TTI state. It only depends on t − t after t = 10 5 : the colored curves collapse. In this data, 100 distinct simulations were averaged, with parameters (S, µ, σ | M, ρ, d, Nc) = (400, 10, 2 | 8, 0.95, 10 −10 , 10 −15 ). quite straightforward. So is the generalization to patchdependent correlation ρ uv . With the assumptions described below, the fraction of coexisting species S * /S is finite when S is large, so that resident diversity S * in each community is also large.
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 µ/C and variance σ 2 /C and set it to A ij . In the large connectivity limit, C → ∞, 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 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 → ∞ limit has been thoroughly studied in physics of disordered systems in the last twenty years [5].
In addition, the paper focuses on the case where migration connects all patches to one another. But the basic DMFT framework, Eq. (B1), is valid even if only certain patches are connected, and migration is zero otherwise. This can allow for analysis of different spatial connectivities, such as lattices reperesenting finite-dimensional space, and is an interesting direction for future work.

Appendix C: Extinction probability of a species
Here the probability of extinction of a species is presented, at the limit N c D 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 ξ u is uncorrelated between different patches, see Appendix B, Fig. A. 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 ≡ ln N u . The equations of the DMFT, Eq. (B1), become Here D = d/ (M − 1). We look for a rare realization of {ξ u } that makes all the x u reach x c = ln N c , in the case when the cut-off is low, x c → −∞. The calculation proceeds within the framework of large-deviation theory [6]. First, one defines the "action" with ξ u substituted with its value from Eq. (C1), 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 ∼ e Jmin with J min the action J minimized over all population trajectories {x u (t)} u=1..M that start at t → −∞ at the typical value of x u , obtained by the zero-noise fixed point of Eq. (C1), 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 . If N * 1 < 0 the species is extinct. On the other hand, if N * 1 > 0, then we obtain the known result [7,8] The result for all M is a generalization of this result, of the form To describe the calculation of N * eff , order the patches so that N * and where w satisfies: w ≤ N * u for all u ≤ m, and w > N * u for all u > m. Such a partition can be shown to always exist. Then The derivation will be given elsewhere.
Consider the case of low migration, D → 0 + . We now develop a theory assuming that the amplitude of the endogenous fluctuations, remains finite in the limit D → 0 + . Assume the species survives, i.e. there is at least one patch with N * Taking the time average of the above equation and therefore N u = N * u + O (D). The previous arguments lead to the conclusion that in the D → 0 + limit 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 v∼u Nv Nu − 1 will be finite, so that D v∼u Nv Nu − 1 will indeed be small. This might break if N u can be small while some other N v remains O (1). 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 − M D N * u , its carrying capacity alone. If patch u fluctuates alone, then For any given N * u this is finite. It diverges as N * u → 0. Therefore the migration term is negligible only if DW (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 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" otherwise 1 . In order to understand the correlation between the sources in the different communities, we unpack N * u using Appendix B. Taking the timeaverage is equivalent to averaging over the dynamical noise ξ. Therefore, in patch u for species i, z i,u = −σS −1/2 j a ij,u N j,u = −σS −1/2 j,+ a ij,u N * j,u . The 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: where we used the Kronecker symbol δ ik . Therefore: The 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 . But the contribution of this migration is small and not required for its persistence.
where we recall m u = N i,u = N * i,u + . We can now compute the different moments of the multivariate Gaussian random variable N * u , using equation (D1). We obtain the closure: When u = v, as ρ uu = 1, we find the expected single community result. In particular, mean [N * u ] and variance [N * u ] do not depend on the patch u. We numerically solve the closure in a self-consistent way: start with a guess for N * i,u N * i,v + , and then (1) Produce many samples of the vector N * u=1..M and (2) calculate the next estimate for N * i,u N * i,v + , 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 covariance [N * u , N * v ], 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 ρ N * of the N * u 's. We use here our simple case of a uniform correlation ρ a between patches ρ uv = ρ a + (1 − ρ a )δ uv . We introduce the notation ρ a instead of 'ρ' in this section in order to avoid confusion with ρ N * .
The results are surprising: even when ρ a → 1, the overlap between communities is not perfect (ρ N * < 1), so the total diversity is larger than the one in each patch. This happens exactly at the transition to chaos at σ c = √ 2, see Fig. B.
On Fig. D, 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 σ, and for ρ a closer to one. To look further into this difference, it is useful to study diversity as a function of the value of N * 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 * eff , which are precisely the species that are more likely to go extinct, with good agreement with theory at higher values of N * eff . This is demonstrated in Fig. C, which shows that the theoretical prediction for the number of species with N * 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,  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 corr [A ij , A ji ] = 0, results in σ c = √ 2 and with 1/2 of the species being sourced in each patch [9]. From top to bottom, we consider three different observables: the diversity, the covariance in the abundances across distinct patches, and this covariance rescaled by the one patch variance. By varying σ, we can control the state of the system: on the left (σ = 1), we show the results for fixed points; on the right (σ = 2 > √ 2), we show the results for persistent dynamical fluctuations. In dotted lines, we plot the theory predictions, as functions of the correlation between patches' interactions ρa. We compare them to simulations with parameters (S, M, µ) = (400, 8, 10), and obtained by simulations run until final time t f = 10 4 . We eventually vary the couple (D, Nc). We use 50 distinct samples of the simulations for each combination of parameters, in order to get error bars and relevant statistics. The cut-off is implemented via patch-wise extinctions when the abundance goes below the threshold in each particular patch, in which case migration out of the patch is turned off while still allowing inward migrations. On the left side, we can see that the theory is exact in the fixed point regime. In this regime, as ρa → 1, the predictions are equivalent to the one patch M = 1 theory, as all patches are the same. In the persistent fluctuation state, the theory is only a good approximation. More precisely, the predictions become more accurate as D and Nc go to zero, as expected. In addition, the agreement gets worse when ρa → 1, because synchronization can occur. In the top right figure, we show that the prediction for diversity is an upper bound. In the bottom right figure, we see that indeed the prediction for ρn is still far from 1 when ρa → 1, for the values of D, Nc used in the simulations. 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. E and Fig. F we show results of a numerical solution [2] 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 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   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 longlived 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 = 0 if and only if A ji,u . We define γ the correlation of the non-zero elements γ = corr [A ij,u , A ji,u ] Aij,u =0 .  shows two simulations, one with γ > 0 and the other with γ < 0. In both cases the system relaxes to a long-lived state with fluctuating abundances, without further loss of diversity up to time 2 · 10 5 . They are intended solely to demonstrate that conditions with γ = 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