Interactions and Migration Rescuing Ecological Diversity

I. INTRODUCTION

Community ecology explores how the interactions between different species shape the diversity-rich ecosystems that characterize the natural world. Understanding the main mechanisms at play is a challenge that spans different scientific fields and it is relevant for human health .

There are three salient facts that one has to take into account in this endeavor. Many ecosystems of interest are species-rich. The interactions between these large sets of species, and the induced ecological dynamics, can lead to complex dynamical behaviors such as chaos and a very large number of possible equilibria . Many ecosystems are spatially extended: the ecological dynamics takes place at some local scale, but individuals can then explore different spatial locations through migration . This can lead to the appearance of complex ecological phenomena, such as traveling activity fronts, pattern formation, and persistent chaotic dynamics . Ecosystems are subject to noise, particularly environmental and demographic noise (due to stochasticity in births and deaths). Both noises induce fluctuations that are a key factor in determining abundance distributions, and their time-dependence . Understanding the interplay between these three properties of ecosystems is essential for answering many central questions in community ecology.

In this work, we consider spatially extended species-rich ecosystems subject to demographic noise. We will consider populations that are large but spatially structured, so that demographic fluctuations globally average out, but they have an important effect on the local dynamics. This is the case, for example, in semiarid ecosystems: the total number of plants is such that global fluctuations are negligible, but at the local level stochasticity can play a fundamental role . Our aim is to understand how in these cases interactions and spatial migration can allow for large diversity and finite abundances despite the adversarial role of demographic noise. In fact, in an isolated community, demographic noise leads to extinctions, irreversibly reducing the ecosystem's diversity until there are no species left .

Previous works, following the classical theory of Island Biogeography by MacArthur and Wilson , proposed as a rescuing mechanism the immigration from a static reservoir (or “mainland,” when thinking of an island-mainland system) . Nevertheless, this approach simply shifts the question from how diversity is maintained on the island to its maintenance on the mainland. Here we use a different approach. We consider ecosystems as a network of ecological communities (a metacommunity) coupled by passive dispersal. In this case, the immigration rates are not externally imposed, but they are the result of the internal dynamics. If a species goes locally extinct in one of the communities, immigrants from the neighboring ones can reinvade, providing an “insurance” (or “storage”) effect . This makes the possibility of a global extinction much more unlikely, and it can allow the ecosystem to self-sustain at finite abundances and diversity. The stabilization of high-diversity states by spatial structure is a very general phenomenon: it can arise in the presence of spatial heterogeneity of environmental conditions or when abundances in different spatial locations exhibit unsynchronized fluctuations . Providing a theory for this mechanism for species-rich ecosystems subject to demographic noise, and assessing the role of interactions, is the main contribution of this work.

The situation is well understood in the case of a few species, in which, depending on the competition between migration and death-birth rates, the system is found to be either in a survival or in an extinct state. A transition separates the two regimes . This phase transition falls in the universality class of Directed Percolation, a second-order out-of-equilibrium transition studied in statistical physics and widely used to describe spreading phenomena, from forest fires to epidemics .

In a many-species metacommunity with constant competitive interactions, it was recently shown that a similar second-order phase transition takes place and that it also belongs to the Directed Percolation universality class . Because the transition is continuous with vanishing abundances, interactions, which are quadratic in the abundances, are subleading at the critical point. As a consequence, the main mechanism at play in this case is still the competition between migration and death-birth rates. We shall show that the scenario for heterogeneous interactions is different and goes beyond the directed percolation paradigm. The transition can become discontinuous. The ecosystem can exhibit global bistability and tipping points between drastically different alternative states. Upon small changes in the environmental condition, the system can therefore undergo catastrophic shifts from a state with large diversity and finite abundances to one in which all species are extinct. As in many other dynamical systems, from coral reefs to arid ecosystems and from Earth's climate to financial markets , it is important to find early warning signals of these kinds of transitions in order to prevent them. We have identified a specific probe, which is based on the response of the ecosystem to perturbations, and that can be monitored in experiments. Our analytical framework shows that interactions play a key role both in the overall scenario and in promoting a self-sustained survival state, in agreement with results obtained for constant mutualistic interactions . Remarkably, in our case, heterogeneous interactions of the pool of species are not necessarily mutualistic on average. It is the ecological dynamics that shapes the ecosystem in a self-sustained phase characterized by emergent mutualistic behavior among the nonextinct species.

In our work, we make use of several methods developed in statistical physics that are particularly well suited for species-rich ecosystems, which are complex systems formed by many interacting degrees of freedom undergoing stochastic dynamics. To model the heterogeneity in the interactions, we sample the coupling coefficients from a random ensemble. We have thus to deal with “disordered” ecosystems, which can be analyzed by transferring methods from spin-glass theory . This disorder approach, which dates back to May's seminal paper , has recently inspired a growing body of work and also received positive experimental confirmations . Previous works have explored within this framework the effect of heterogeneous interactions , demographic fluctuations , and spatial structure , but the analysis we present here is, to our knowledge, the first analytical study in which the three ingredients are combined.

The model we focus on is a disordered Generalized Lotka Volterra (GLV) system of metacommunity subject to demographic noise. For one community, the disordered GLV has been shown to have a rich phase diagram, and to display several dynamical regimes: single equilibrium, multistability, and chaos . We expect this complex behavior also in the case of spatially structured ecosystems . In this work, we focus on the moderate-heterogeneity regime in which there is a single stable equilibrium. This allows us to disentangle the multistability due to the fragmentation of the basins of attraction of the ecological dynamics at strong heterogeneity from the bistability of the feedback mechanism between abundance and immigration. Our analysis is performed using a mean-field approximation on the spatial fluctuations, which is equivalent to considering that the community network is a fully connected graph.

Note that because of their generality, Lotka-Volterra equations have been applied to a variety of fields besides their original ecological interpretation, from immunology to economics and game theory . Our results could therefore find applications beyond ecology, notably for the study of global bistability and economic crashes.

II. THE MODEL

We consider a metacommunity of $𝑆$ species living on a network of $𝐿$ discrete spatial locations, or patches. A graphical representation of the system is given in Fig. in the case of a fully connected network of three patches. Each species is characterized by its abundance in each patch, which is modeled by a continuous variable, $𝑁 𝑖, 𝑢$ , representing the total number of individuals divided by the typical size of the local population, $˜ 𝑁 typ$ . The abundance of species $𝑖$ in patch $𝑢$ evolves according to the stochastic differential equation:

̇ 𝑁 𝑖, 𝑢 = 𝑟 𝑘 𝑁 𝑖, 𝑢 (𝑘 - 𝑁 𝑖, 𝑢 - \sum 𝑗 𝛼 𝑢 𝑖 𝑗 𝑁 𝑗, 𝑢) + 𝐷 𝑐 \sum 𝑣 \in 𝜕 𝑢 (𝑁 𝑖, 𝑣 - 𝑁 𝑖, 𝑢) + 𝜂 𝑢 𝑖 (𝑡) \sqrt 𝑁 𝑖, 𝑢,

(1)

which corresponds to Lotka-Volterra dynamics, with constant growth rate $𝑟$ and carrying capacity $𝑘$ that are set to 1 throughout. The notation $𝜕 𝑢$ indicates the set of patch neighbors of $𝑢$ (from and to which species in the patch $𝑢$ can migrate). The growth of each species is influenced by the abundance of all the others through the interaction coefficients $𝛼 𝑢 𝑖 𝑗$ : if $𝛼 𝑢 𝑖 𝑗$ is positive, species $𝑗$ inhibits the growth of species $𝑖$ in patch $𝑢$ and vice versa. Positive $𝛼 𝑢 𝑖 𝑗$ and $𝛼 𝑢 𝑗 𝑖$ correspond to two species competing for resources, whereas $𝛼 𝑢 𝑖 𝑗$ and $𝛼 𝑢 𝑗 𝑖$ both negative correspond to mutualistic behavior. Predation leads to opposite signs.

FIG. 1.

A metacommunity of seven species living on three patches. Each individual interacts with the local community to which it belongs possibly migrating to neighboring patches with diffusion coefficient $𝐷$ .

To model the heterogeneity in the interactions of species-rich ecosystems, we follow and consider the disordered LV model. As already discussed in the Introduction, the disorder approach has attracted recently a lot of attention and also received positive experimental confirmations . In this framework, the interaction coefficients are random variables, with mean $𝜇 / 𝑆$ and variance $𝜎 2 / 𝑆$ . They are independent in each patch except for $𝛼 𝑢 𝑖 𝑗$ and $𝛼 𝑢 𝑗 𝑖$ , which have a correlation coefficient $𝛾$ . In the following, we will first focus on the symmetric interactions case ( $𝛾 = 1$ ), and then show that a small asymmetry does not qualitatively change the results. As the interactions between species can depend on the environmental conditions (temperature, humidity, resources availability, etc.) which differ in space, we consider interaction matrices fluctuating from one patch to another, i.e., they are not identical in different patches, but corresponding elements $𝛼 𝑢 𝑖 𝑗$ and $𝛼 𝑣 𝑖 𝑗$ have a correlation coefficient $𝜌$ .

We will restrict the choice of $𝜇$ and $𝜎$ to values for which an isolated Lotka-Volterra community only displays a single uninvadable equilibrium (the single equilibrium phase studied in Ref. ). Without spatial heterogeneity, the transition point is not modified by the introduction of a spatial structure , and spatial heterogeneity decreases the effective complexity of the interaction network , favoring the single equilibrium phase. Therefore, we also expect the metacommunity to be in the single equilibrium phase for all the allowed values of $𝜇$ and $𝜎$ . The effect of migration between patches in the strong heterogeneity regime with nonsymmetric interactions, in which a single community with fixed immigration exhibits chaotic dynamics , was studied in in the absence of demographic noise. It leads to complex dynamical behavior with long-lived persistent fluctuations. Combining strong heterogeneity, demographic noise, and spatial migration is a challenge left for future studies.

In the model defined by Eq. , individuals can migrate on the patches network through diffusion, with a constant diffusion coefficient $𝐷 / 𝑐$ , where $𝑐$ is the connectivity (or number of connections per site) of the network. We assume the network to be translationally invariant, therefore each site has the same connectivity. Migration is possible and equiprobable from patch $𝑢$ to any of its $𝑐$ nearest neighbors $𝑣 \in 𝜕 𝑢$ .

Each species is subject to a white demographic noise $𝜂 𝑢 𝑖$ , accounting for the stochasticity in birth and death events in a continuum setting . We follow Ito's convention, according to which fluctuations in birth and death at time $𝑡 + 𝑑 𝑡$ depend on the abundance at the previous time step. The noise is uncorrelated and of constant amplitude across species and patches:

⟨ 𝜂 𝑢 𝑖 (𝑡) 𝜂 𝑣 𝑗 (𝑡') ⟩ = 2 𝑇 𝛿 𝑖 𝑗 𝛿 𝑢 𝑣 𝛿 (𝑡 - 𝑡') .

(2)

The autocorrelation of the demographic noise defines the noise strength $𝑇$ , which depends on the birth and death rates and on the typical size of the local population; $𝑇$ scales as $𝑇 \propto 1 / ˜ 𝑁 typ$ : the larger the local populations, the more negligible are demographic fluctuations. In the $𝛾 = 1$ case, $𝑇$ can be interpreted as an effective temperature, as we shall show later.

Some further insights into the effect of the demographic noise can be obtained considering it in the absence of all the other terms. In this case, an exact solution to the associated Fokker-Planck equation is available, showing that starting from any initial condition the population goes to 0 abundance with some finite rate . Therefore, also in the continuous model, extinction is possible over finite times, and not only asymptotically, as would be the case, for example, with environmental noise.

The dynamics of species in the presence of birth and death has important connections with the celebrated directed percolation problem studied in out-of-equilibrium physics and statistical field theory . Directed percolation is a model of particles that hop on a network and are subjected to births and deaths; a graphical illustration of the process can be found in Fig. for a one-dimensional network. Directed percolation was originally introduced to model spreading phenomena, from forest fires to epidemics . In our case, the sites of the network represent spatial locations, or patches, on which (or from which) species can migrate; the particles indicate which sites are colonized by species. At each time step the particles can produce an offspring in a neighboring site, die, or just survive. In our case, this corresponds to colonization or extinction. Depending on the competition between death and birth rates, the activity can spread to the entire system and lead to a finite density of particles (active, self-sustaining state) or die out (absorbing, inactive state). Between these two phases, there is a continuous phase transition, characterized by universal critical behavior . We show in Fig. the phase diagram in the mean-field approximation (discussed in the next section). A direct link between DP and GLV is obtained by coarse-graining . In this way, the discrete DP occupation variable becomes a continuous quantity that represents the mean occupation, the competition between birth and death rates gives rise to a logistic growth, hopping is replaced by diffusion, and the stochastic fluctuations generate the demographic noise. This leads to a set of independent GLV Eqs. in the absence of interactions, one for each species. Each equation corresponds to an independent directed percolation process.

FIG. 2.

Directed percolation on an array of seven sites. Each row represents a different time step, green arrows indicate birth, gray arrows death, and orange arrows survival.

FIG. 3.

Phase diagram for Directed Percolation in the mean-field approximation: in green the active phase, in which at long times there is a finite density of particles, in white the inactive phase, in which all particles eventually die. $𝐷 0 (𝑇)$ indicates the transition line (see Sec. and Appendix for details).

The directed percolation transition can therefore be interpreted as a transition between a self-sustained phase where migration enables a finite abundance of species to persist, to a regime, characteristic of small (or zero) dispersal, where species go extinct due to demographic noise. The aim of this work is to develop a theory for these phenomena for species-rich ecosystems in the presence of heterogeneous interactions. Upon increasing the number of species in the pool and considering heterogeneous interactions, the set of directed percolation processes is no longer independent, and the complexity of the model increases considerably. In fact, the system becomes equivalent to the collection of an infinite number of directed percolation processes, coupled by random interactions—an interesting and open statistical physics problem.

III. METHODS

A. DMFT and coupled Directed Percolation processes

In this work, we aim to study systems in which both the number of species and the number of patches are very large. To obtain analytical results, we follow the statistical physics “way” and take the limit of an infinite number of species and an infinite number of patches. In this double limit (whose order is irrelevant) the macroscopic properties of the system do not depend on the particular realization of the demographic noise and of the interactions: the macroscopic properties are self-averaging in the jargon of disordered systems .

The large- $𝑆$ limit allows for an analytical treatment, as the dynamics of the $𝑆$ interacting degrees of freedom can be replaced by the effective dynamics for a single representative species, through Dynamical Mean Field Theory (DMFT) . The interaction effect with other species is captured by a noise term, which can be seen as an environmental noise (or a thermal bath) statistically defined in a self-consistent way. The DMFT procedure is analogous to the one used to derive Langevin's equation from Newtonian dynamics , with the difference that here the degrees of freedom that are integrated out, giving rise to the noise, are equivalent to the degree of freedom under consideration, thus allowing a self-consistent closure of the equations of motion. DMFT is a very powerful technique that has been employed in several different contexts from quantum many-body systems to glassy dynamics . Thanks to DMFT, we can map an infinite number of randomly coupled DP processes—a formidable problem—into a single DP process with additional terms to be determined self-consistently (a colored noise and a memory term).

Our derivation follows the one developed in Ref. for LV models, and it can be found in Appendix for generic values $𝜌$ of the spatial heterogeneity of the interactions. Here we outline the main steps in the special case of patch-independent interactions, $𝜌 = 1$ . In the following, we are interested in the steady states of the dynamics. In fact, we expect that after a transient, the system will settle in a time translationally invariant regime. For $𝑆 \to \infty$ , DMFT allows one to replace the interaction term $- \sum 𝑗 𝛼 𝑖 𝑗 𝑁 𝑗, 𝑢$ by a stochastic expression that has the same statistical properties:

- 𝜇 ℎ - 𝜎 ̃ 𝜉 𝑖 𝑢 (𝑡) + 𝜎 2 𝛾 \int 𝑡 0 \sum 𝑣 𝑅 𝑢 𝑣 (𝑡, 𝑠) 𝑁 𝑖 𝑣 (𝑠) 𝑑 𝑠 .

(3)

Since this allows us to decouple different species, for simplicity we will omit the species index $𝑖$ in the following. We now discuss the different contributions. Note that in the following, empirical averages over species will be denoted as $𝔼 [\cdot]$ .

The first term represents the average interaction with all other species. It is given by the product of the mean of the interaction strength and the mean abundance, $ℎ = 𝔼 [𝑁 𝑢]$ , that in the steady state does not depend on the patch $𝑢$ thanks to translational invariance.

The second term represents the fluctuation of the interaction with all other species. It is given by the product of the standard deviation of the interaction coefficients and Gaussian noise with zero mean and correlation matching the time autocorrelation of the single species abundances:

⟨ ̃ 𝜉 𝑢 (𝑡) ̃ 𝜉 𝑣 (𝑠) ⟩ = 𝔼 [𝑁 𝑢 (𝑡) 𝑁 𝑣 (𝑠)] : = 𝐶 𝑢, 𝑣 (𝑡 - 𝑠) .

(4)

The noise $̃ 𝜉 𝑢 (𝑡)$ is multiplied by the abundance in the LV equations. Henceforth we will call it environmental since its effect is to add fluctuations to the carrying capacity. Since the autocorrelation of the abundances generically decays to a positive plateau at large time separations , one can decompose the environmental noise into a fluctuating component and a static one. The former corresponds to the fluctuations due to ecological dynamics for a given species. The latter is characteristic of a given species and fluctuates from species to species . We decompose the noise by rewriting $̃ 𝜉 𝑢 (𝑡) = 𝑧 \sqrt 𝐶 \infty 𝑑 + 𝜉 𝑢 (𝑡)$ , where $𝐶 \infty 𝑑 = l i m 𝜏 \to \infty 𝐶 𝑢, 𝑢 (𝑡, 𝑡 + 𝜏)$ is the value of the correlation function within the same patch at infinite times, $𝑧$ is a static Gaussian variable with zero mean and unit variance, which now plays the role of quenched disorder, and $𝜉 𝑢 (𝑡)$ is a fluctuating noise whose covariance vanishes at long times. Again $𝑧$ and $𝐶 \infty 𝑑$ do not depend on the patch $𝑢$ thanks to translational invariance.

To distinguish the roles of fluctuating and static noises in the GLV equation, we introduce two kinds of averages: $⟨ \cdot ⟩$ indicates the average over the fluctuating noises $𝜉$ and $𝜂$ . It is an average over the ecological dynamics, or by ergodicity, over patches for a fixed species. In analogy with physical systems, we call it thermal average. The overline $‾ \cdot$ instead stands for the average over the static field $𝑧$ ; it corresponds to averaging over species or over different instances of the interaction matrix. Again in analogy with the physical system, we call it quenched disorder average.

The last term in the dynamical mean-field treatment of the interactions is due to a feedback mechanism: a fluctuation of the abundance of species $𝑖$ influences species $𝑗$ , which in turn influences species $𝑖$ . These contributions sum up because of the correlation between $𝛼 𝑖 𝑗$ and its reciprocal $𝛼 𝑗 𝑖$ , leading to the factor $𝛾$ . This feedback mechanism (the famous Onsager reaction in the spin-glass literature) generates a memory term, containing the response function of the abundance on patch $𝑢$ to a perturbation in the carrying capacity in patch $𝑣$ :

𝑅 𝑢, 𝑣 (𝑡, 𝑠) = 𝔼 [𝛿 𝑁 𝑢 (𝑡) 𝛿 𝜁 𝑣 (𝑠) | 𝜁 = 0] .

(5)

In the $𝑆 \to \infty$ limit, there is convergence in law between the statistics of the infinite number of randomly coupled DP processes and the effective one , i.e., the dynamics of a species satisfying the GLV equation is equivalent to the effective one of a single species living on the original spatial network:

̇ 𝑁 𝑢 = 𝑁 𝑢 (𝑘 - 𝑁 𝑢 - 𝜇 ℎ - 𝜎 (𝑧 \sqrt 𝐶 \infty 𝑑 + 𝜉 𝑢) + 𝜎 2 𝛾 \int 𝑡 0 \sum 𝑣 𝑅 𝑢 𝑣 (𝑡, 𝑠) 𝑁 𝑣 (𝑠) 𝑑 𝑠) + 𝐷 𝑐 \sum 𝑣 \in 𝜕 𝑢 (𝑁 𝑣 - 𝑁 𝑢) + 𝜂 𝑢 \sqrt 𝑁 𝑢 .

(6)

The DMFT closure consists then in replacing the empirical averages over species $𝔼 [\cdot]$ with the one with respect to the effective single-species one. Because the effective process itself depends on some averaged quantities, one ends up with a self-consistent stochastic equation.

Equation can also be interpreted as the Langevin equation associated with a Directed Percolation (DP) process, with the addition of a memory term (that is absent in the special case $𝛾 = 0$ ) and environmental noise. The effect of the static part of the environmental noise $𝑧$ is to change the control parameter of the DP process, determining whether this is subcritical or supercritical.

Interestingly, whereas a system of a few species interacting and diffusing on a network was established to boil down to a standard DP problem , the case of many species is fundamentally different and belongs to a different class. Indeed, a system of many species is equivalent to a family of many DP processes, characterized by different values of static and fluctuating noises and coupled through the common self-consistently determined mean, correlation, and response functions. Understanding the behavior of this self-consistent DP problem is an open challenge. In this work, we study whether the DP transition can fundamentally change nature due to this self-consistent coupling. Even if the transition remained qualitatively DP-like (continuous and from an absorbing state to a fluctuating one), critical properties could change. In fact, although an environmental noise can be shown to be an irrelevant perturbation of the associated field theory , within DMFT the environmental noise inherits the time dependence of the correlation function through the self-consistency. It can therefore develop long-range correlations in time at the critical point, possibly altering the critical behavior and leading to a new universality class.

B. Symmetric interactions, mean-field approximation, and mapping to a system in thermal equilibrium

Studying the coupled field theories is a formidable task. In the following, we simplify the problem by doing a mean-field approximation, which allows us to obtain a general theory independent of the underlying network of patches.

We replace the term $𝐷 𝑐 \sum 𝑣 \in 𝜕 𝑢 𝑁 𝑣$ by its thermal average. This amounts to $𝐷 𝑐 \sum 𝑣 \in 𝜕 𝑢 𝑁 𝑣 \to 𝐷 𝑁 *$ , where $𝑁 * = 1 𝑐 \sum 𝑣 \in 𝜕 𝑢 ⟨ 𝑁 𝑣 ⟩$ , and, using translation invariance, it simplifies to $⟨ 𝑁 𝑢 ⟩$ (which is time-independent since we are considering steady states). This procedure corresponds to a mean-field approximation of the spatially dependent DMFT Eqs. . Such a $D M F T 2$ approximation becomes exact for a fully connected network. In fact, in this case, taking the $𝐿 \to \infty$ limit, the empirical average of the abundances over the patches concentrates around the thermal average $𝑁 * = ⟨ 𝑁 𝑢 ⟩$ . From now on, we shall focus on this case.

By substituting $𝐷 𝑐 \sum 𝑣 \in 𝜕 𝑢 𝑁 𝑣$ with $𝑁 *$ in Eq. , one obtains an equation on $𝑁 𝑢$ only, with an additional parameter to be determined self-consistently. Note that $𝑁 *$ is obtained by averaging only over thermal fluctuations, and not over disorder: therefore, it will have to be determined as a function of $𝑧$ . This means that different species will have different immigration rates (here, for simplicity, we are still focusing on the $𝜌 = 1$ case; generalizations will be discussed later).

This substitution allows us to decouple stochastic processes for the abundance in different patches. Omitting for simplicity the index $𝑢$ , we now obtain (for large times, i.e., in the steady state)

̇ 𝑁 = 𝑁 (𝑘 - 𝑁 - 𝜇 ℎ - 𝜎 𝑧 \sqrt 𝐶 \infty 𝑑 - 𝜎 𝜉 (𝑡) + 𝜎 2 𝛾 (\int 𝑡 0 𝑅 𝑑 (𝑡 - 𝑠) 𝑁 (𝑠) 𝑑 𝑠 + 𝑁 * \int 𝑡 0 𝑅 0 (𝑡 - 𝑠) 𝑑 𝑠)) + 𝐷 (𝑁 * - 𝑁) + 𝜂 (𝑡) \sqrt 𝑁 .

(7)

Since all patches are equivalent on a fully connected lattice, the $𝑅 𝑢 𝑣$ and $𝐶 𝑢 𝑣$ matrices (of functions) only have two independent elements: the diagonal ones, $𝑅 𝑑$ and $𝐶 𝑑$ , and the off-diagonal ones, $𝑅 0 / 𝐿$ and $𝐶 0$ (see Appendix for the justification of the scaling with $𝐿$ of $𝑅 0 / 𝐿$ and $𝐶 0$ ).

In the case of symmetric interactions, $𝛾 = 1$ , one can show (see Appendix ) that the self-consistent solution maps to a thermal equilibrium process. In fact, one finds that the diagonal elements of the response and correlation functions obey a fluctuation-dissipation relation:

𝑅 𝑑 (𝜏) = - 1 𝑇 𝜕 𝜕 𝜏 𝐶 𝑑 (𝜏) .

(8)

The memory term and $𝜉$ therefore play the role of a friction term and the noise associated with a colored thermal bath at temperature $𝑇$ . The stochastic process maps then to a generalized Langevin equation whose stationary probability distribution is given by the Boltzmann distribution at temperature $𝑇$ and with the effective Hamiltonian:

𝐻 eff = (1 - 𝜎 2 𝑇 (𝐶 0 𝑑 - 𝐶 \infty 𝑑)) 𝑁 2 2 - (𝑘 - 𝐷 - 𝜇 ℎ - 𝑧 \sqrt 𝐶 \infty 𝑑 𝜎 + 𝜎 2 𝑁 * 𝑅 int 0) 𝑁 + (𝑇 - 𝐷 𝑁 *) l n 𝑁,

(9)

where $𝐶 0 𝑑$ is the equal-time correlation function, namely the second moment of the abundances over disorder and noise, $‾ ⟨ 𝑁 2 ⟩$ . The long-time limit of the correlation function, $𝐶 \infty 𝑑$ , represents instead the second moment of the thermal-averaged abundances, $‾ ⟨ 𝑁 ⟩ 2$ . $𝑅 int 0$ is the integrated off-diagonal response, which is the solution of the self-consistent equation (see Appendix ):

𝑅 int 0 = ‾ 𝑟 𝑑 (𝑧) 𝐷 𝜒 (𝑧) + 𝜎 2 𝑅 int 0 𝑟 𝑑 (𝑧) 1 - [𝐷 𝜒 (𝑧) + 𝜎 2 𝑅 int 0 𝑟 𝑑 (𝑧)] .

(10)

$𝜒 (𝑧)$ and $𝑟 𝑑 (𝑧)$ are the species-dependent response to a perturbation in the immigration rate or the carrying capacity, respectively:

𝜒 (𝑧) = ⟨ 𝑁 l o g 𝑁 ⟩ - ⟨ 𝑁 ⟩ ⟨ l o g 𝑁 ⟩,

(11)

𝑟 𝑑 (𝑧) = ⟨ 𝑁 2 ⟩ - ⟨ 𝑁 ⟩ 2 .

(12)

The self-consistent equations can be expressed as averages with respect to the Boltzmann distribution:

𝑁 * (𝑧) = ⟨ 𝑁 ⟩ = \int \infty 0 𝑑 𝑁 𝑁 𝑒 - 𝛽 𝐻 eff \int \infty 0 𝑑 𝑁 𝑒 - 𝛽 𝐻 eff,

(13)

ℎ = ‾ ⟨ 𝑁 ⟩ = \int 𝒟 𝑧 \int \infty 0 𝑑 𝑁 𝑁 𝑒 - 𝛽 𝐻 eff \int \infty 0 𝑑 𝑁 𝑒 - 𝛽 𝐻 eff,

(14)

𝐶 0 𝑑 = ‾ ⟨ 𝑁 2 ⟩ = \int 𝒟 𝑧 \int \infty 0 𝑑 𝑁 𝑁 2 𝑒 - 𝛽 𝐻 eff \int \infty 0 𝑑 𝑁 𝑒 - 𝛽 𝐻 eff,

(15)

𝐶 \infty 𝑑 = ‾ ⟨ 𝑁 ⟩ 2 = \int 𝒟 𝑧 (\int \infty 0 𝑑 𝑁 𝑁 𝑒 - 𝛽 𝐻 eff \int \infty 0 𝑑 𝑁 𝑒 - 𝛽 𝐻 eff) 2,

(16)

and analogously for $𝑅 int 0$ . $\int 𝒟 𝑧 = \int 𝑑 𝑧 \sqrt 2 𝜋 𝑒 - 𝑧 2 / 2$ indicates the average over the Gaussian field.

These equations can be solved iteratively: starting from a suitable initial condition for $𝑁 * (𝑧), ℎ, 𝐶 0 𝑑, 𝐶 \infty 𝑑$ , and $𝑅 int 0$ , one updates their values according to Eqs. until reaching a fixed point. Because very large values of $𝑧$ are exponentially suppressed by the Gaussian distribution, it is sufficient to determine $𝑁 * (𝑧)$ for $𝑧$ of $𝑂 (1)$ .

In conclusion, within the $D M F T 2$ approximation and for the symmetric case, the formidable self-consistent stochastic equations can be analyzed by studying a set of static self-consistent equations on four parameters $ℎ, 𝐶 0 𝑑, 𝐶 \infty 𝑑, 𝑅 int 0$ and one function $𝑁 * (𝑧)$ . Solving these equations (see the next section) allows us to obtain a general picture of the interplay between migration and demographic noise for spatially extended metacommunities. To show that such a picture is valid beyond the simplified case we focus on, we have also considered several extensions that we shall present below.

C. Extensions

1. Spatial heterogeneity

In the case of a generic value of the spatial heterogeneity of the interactions $𝜌$ , an analogous procedure can be implemented, with some important differences. The static disorder is now a patch-dependent and correlated variable, which we can decompose as $𝜌 \sqrt 𝐶 \infty 0 𝑧 + \sqrt 𝐶 \infty 𝑑 - 𝜌 2 𝐶 \infty 0 𝑤 𝑢$ , where $𝑧$ is constant and $𝑤 𝑢$ is independent across locations, and $𝐶 \infty 𝑑$ and $𝐶 \infty 0$ are the infinite time correlation function of the abundance on the same patch and on different patches, which coincide for $𝜌 = 1$ . Averaging the abundance across patches to obtain the immigration rate requires an additional step, i.e., averaging also over $𝑤 𝑢$ . The solution of the self-consistent equations, albeit conceptually analogous to the $𝜌 = 1$ case, is much more numerically challenging for generic values of $𝜌$ because of the need to integrate over two disorder fields, $𝑧$ and $𝑤 𝑢$ . For this reason, we focused on the two extreme cases, $𝜌 = 1$ and 0, in which only one of the two disorder fields is present. The results are qualitatively similar, so we expect our conclusions to hold also for intermediate values of $𝜌$ . We confirm it by numerical simulations at $0 < 𝜌 < 1$ .

2. Asymmetric interactions

The mapping to an equilibrium distribution requires symmetry in the interactions: nonsymmetric interactions correspond to nonconservative forces, which explicitly break time reversal and lead to nonequilibrium steady states. To show that our results hold also in this case, at least if the asymmetry is not too strong, we have analyzed the case of small asymmetry in perturbation theory. The analysis of the Martin–Siggia–Rose–De Dominicis–Janssen action allows us to conclude that a small degree of asymmetry ( $𝛾 = 1 - 𝜀, 𝜀 ≪ 1$ ) does not affect qualitatively the results we shall present in the next section, therefore establishing that our findings for the symmetric case also holds for small asymmetry (see Appendix for more details). We have also confirmed this result by numerical simulations for $𝛾 < 1$ .

IV. RESULTS

In the following, we present our analytical results focusing on ecosystems with parameters $𝜎 = 0.5$ and $𝜇 = 1$ , hence a case in which interactions are in average competitive for the pool of species.

A. Characterization of the self-sustained phase

By solving the DMFT equations described in the previous section, one finds that when the diffusion constant is large enough, the system is in a self-sustained phase (active phase in the directed percolation jargon) in which a nonzero abundance is maintained despite the presence of demographic fluctuations. In this regime, although some species go globally extinct on all patches, others survive thanks to the migration from neighboring patches. This mechanism is sufficient to prevent extinctions due to demographic stochasticity and leads to a self-sustained metacommunity.

In the following, we discuss the salient properties of this phase, focusing on two ecologically relevant observables: the average abundance, $ℎ = ‾ ⟨ 𝑁 ⟩$ , and the ecosystem diversity $𝜙$ , defined as the fraction of species that are not globally extinct, i.e., that have nonzero abundance in at least one patch. At stationarity, we can compute the ecosystem diversity as $𝜙 = ‾ 𝜃 (⟨ 𝑁 ⟩)$ .

As expected, demographic noise is detrimental to survival: the fraction of surviving species, or diversity, and the average abundance decrease with the strength of demographic fluctuations; see the bottom panels of Fig. . On the contrary, dispersal is beneficial, as shown in the top panels of Fig. . The behavior of the diversity for species-rich ecosystems with heterogeneous interactions in the presence of demographic noise is a novel result of our approach: in the case of fixed external immigration, previously often considered in the literature, all species are kept alive by the immigration, albeit some at very small abundances, it is therefore not possible to rigorously define the ecosystem diversity . We find that the species that go extinct are those whose growth is on average more affected by the interactions with the rest of the ecosystem, as quantified by the static part of the environmental noise $𝑧 𝜎 \sqrt 𝑞 0$ , which renormalizes the carrying capacity of a species. For $𝜌 = 1$ , if $𝑧$ is larger than a critical value $𝑧 *$ , the corresponding species goes extinct (for $𝑧 > 𝑧 *$ , the renormalized carrying capacity is negative). This is true also for smaller values of $𝜌$ (Appendix ). The case of independent interactions across patches ( $𝜌 = 0$ ) is special, for all species are globally equivalent so that they can only be all surviving or all extinct. In general, all species have some patches in which they are very abundant. Immigrants from these patches can then save them from extinction in the rest of the system. This favorable role of dispersal through which spatial heterogeneity enhances diversity has been discussed in .

FIG. 4.

Average abundance $‾ ⟨ 𝑁 ⟩$ and diversity $𝜙$ as a function of the diffusion constant $𝐷$ for $𝑇 = 0.25$ (top) and as a function of temperature (strength of demographic noise) for $𝐷 = 0.1$ (bottom). The dashed lines represent the $𝑇 = 0$ well-mixed results. $𝜇 = 1, 𝜎 = 0.5$ .

The limits $𝐷 \to \infty$ and $𝑇 \to 0$ can be mapped to the well-mixed case. For $𝐷 \to \infty$ , the timescale of spatial mixing is much smaller than all other timescales, therefore the abundances of each species are equal on all sites. The absence of spatial fluctuations allows one to write an evolution equation involving only the space-averaged abundances, which corresponds to an effective single local community without demographic fluctuations with interactions given by the spatial average of the original ones. The well-mixed result is also recovered (for $𝜌 = 1$ ) in the $𝑇 \to 0$ limit (see the two bottom panels of Fig. ): because the abundances do not fluctuate, there is no migration flux between patches, and the diffusion term plays no role.

As for the distribution of the abundances, we find an exponential decay (see Appendix ), as is the case in other models with random fully connected interactions .

B. Transition to complete extinction: Emergence of a discontinuous transition at low dispersal

When demographic fluctuations are sufficiently strong, decreasing the diffusion constant leads to a continuous phase transition from an active phase in which some species are able to self-sustain, to an inactive phase in which they are all extinct. The critical value of the diffusion constant is the same as would be obtained in the absence of interactions, where the system directly maps to directed percolation, or in the case of constant interactions ; see Fig. and Appendix . This is to be expected: upon approaching the transition, the abundances tend to zero, and therefore the interactions, which have a quadratic dependence on the abundances, become irrelevant. The critical exponents indeed match the ones falling in the Directed Percolation universality class; in particular, the abundance goes to zero linearly [Fig. ]. Interestingly, approaching the transition the diversity does not go to zero and instead tends to a finite value [Fig. ]. The average abundance goes to zero not because more and more species are going extinct, but because all surviving species are simultaneously decreasing their abundances. This homogenization in the behavior of species is yet another consequence of the irrelevance of the interactions, the only trait distinguishing one species from another in our model.

FIG. 5.

(a) The phase diagram for constant interactions across patches ( $𝜌 = 1$ ). The continuous line indicates the continuous transition, and the dotted and dashed lines are the limits of the metastability region, highlighted in gray. At the two limits of the metastability region, one of the two solutions disappears and we have a discontinuous transition. The arrows indicate the parameter range in the right figures. The average abundance $ℎ = ‾ ⟨ 𝑁 ⟩$ and the diversity $𝜙 = ‾ 𝜃 (⟨ 𝑁 ⟩)$ as a function of $𝐷$ across a discontinuous (b), (d) or continuous (c), (e) transition. In subfigure (c) the arrows indicate the direction of the hysteresis cycle: with decreasing $𝐷$ (starting from high values), the ecosystem would follow the finite solution until the discontinuous transition, where the abundances would jump to zero. If we now increase $𝐷$ , it would follow the zero solution until this becomes unstable at $𝐷 0 (𝑇)$ . Gray dashed lines indicate the value of $𝐷$ at which a single species would go (continuously) extinct. Note that we have divided $𝐷$ by the critical value of the diffusion constant for Directed Percolation $𝐷 0 (𝑇)$ in all plots to emphasize the effect of interactions on the already known case. Because $𝐷 0 (𝑇)$ vanishes exponentially for $𝑇 \to 0$ (Appendix ), the metastability region has a vanishing width in this limit and the system is always in the survival phase. $𝜇 = 1, 𝜎 = 0.5$ .

At smaller demographic noise this picture changes drastically and interactions play a major role, as shown in the phase diagram in Fig. . The ecosystem is able to self-sustain at values of the diffusion constant for which in the absence of interactions it would be in the inactive phase. Upon further lowering $𝐷$ , we encounter a discontinuous transition at which all species abruptly go extinct, i.e., species abundances suddenly jump to zero. Before the discontinuous transition, there is an extended region in which the ecosystem is metastable [in gray in Fig. ]: in this regime, the system reaches an equilibrium with high or low abundances depending on the initial conditions. It exhibits hysteresis [Fig. ].

It was recently shown that a metacommunity subject to demographic noise and constant mutualistic interactions exhibits a similar discontinuous phase transition . The authors of also performed numerical simulations with random (patch-independent) interactions, showing that the surviving species have more mutualistic interactions than the total species pool. We find that a similar mechanism is at play in our case: it is an emergent phenomenon due to ecological dynamics which is present even though interactions are not on average mutualistic (in fact they are competitive, $𝜇 = 1$ ). Because of the symmetry in the interaction network, species that interact more competitively are more negatively affected by the interactions with the rest of the ecosystem, and will hence be more easily driven to extinction. This leads to a decrease of the mean of the interaction matrix restricted to surviving species, which we have estimated in the case $𝜌 = 1$ using a result obtained in (Appendix ). Another quantity of interest is the average interaction term for nonextinct species, $‾ Int + = ‾ \sum 𝑗 𝛼 𝑖 𝑗 ⟨ 𝑁 𝑗 ⟩ +$ (the $+$ indicates that the average is carried out only over nonextinct species, $⟨ 𝑁 𝑖 ⟩ > 0$ ), which we find to be negative in the entire metastability region [Fig. ]. For a species to survive in conditions in which without interactions it would go extinct, we need the interaction term (that appears summed to the carrying capacity with a negative sign) to give on average a negative contribution. We indeed find numerically that only species with negative interaction terms manage to survive [Fig. ], thus leading to an enhancement of mutualism between surviving species—see Appendix for details.

FIG. 6.

Thermal averaged interaction term, $Int = ⟨ \sum 𝑗 𝛼 𝑖 𝑗 𝑁 𝑗 ⟩$ , averaged over nonextinct species (indicated by an overline with a $+$ superscript), for two temperatures corresponding to the discontinuous regime. Left: analytical results for $𝑇 = 0.4, 𝜌 = 1$ [as in Figs. –]. $Int +$ is negative in the metastability region; it jumps to zero when all species go extinct at the discontinuous transition. Right: Distribution of the thermal averaged interaction terms in a numerical simulation in the metastability region [ $𝑇 = 0.18, 𝐷 / 𝐷 0 (𝑇) = 0.8, 𝑆 = 200, 𝐿 = 400, 𝑡 max = 500$ , averaged over two runs]. Nonextinct species are highlighted in orange; only species with negative (or close to zero) interaction terms manage to survive. Averaging only over nonextinct species (orange dotted line) leads to a significantly lower (more mutualistic) value than averaging over all species (blue dotted line). $𝜇 = 1, 𝜎 = 0.5$ .

In Fig. we also show the phase diagram in the case of independent ( $𝜌 = 0$ ) interactions across patches, to be compared to the one of Fig. corresponding to constant ( $𝜌 = 1$ ) interactions across patches. In both cases, the upper limit of the metastability region is bounded from below by the critical value of the diffusion constant in the absence of interactions, $𝐷 0 (𝑇)$ . For $𝜌 = 1$ these two lines coincide, whereas for $𝜌 = 0$ the metastability region extends above $𝐷 0 (𝑇)$ in some range of temperature. In the part of the metastability region above $𝐷 0 (𝑇)$ , the two metastable solutions are both finite: one is of order 1 and the other is proportional to the distance from $𝐷 0 (𝑇)$ ; the two solutions coalesce at the tip of the metastability region.

FIG. 7.

The phase diagram for independent interactions across patches ( $𝜌 = 0$ ). The continuous line indicates the continuous transition; the dotted and dashed lines indicate the limits of the metastability region, highlighted in gray. At the two limits of the metastability region, one of the two solutions disappears and we have a discontinuous transition. $𝜇 = 1, 𝜎 = 0.5$ .

One can also analytically show that the phase diagrams remains qualitatively unchanged considering a small asymmetry in the interactions ( $𝛾 = 1 - 𝜀, 𝜀 ≪ 1$ ); see Appendix . Numerical simulations presented in the next sections confirm this result.

V. ASSESSING THE GENERALITY OF THE SCENARIO

To confirm the generality of our results, we now consider different variations of the model studied in the previous section. The aim is to show that our results hold in a broader setting. We shall be particularly interested in considering the case of a large but finite number of species, a large but finite number of patches, a small but finite asymmetry of interactions, as well as intermediate values of $𝜌$ . All these cases could in principle be studied analytically, but they would require very involved (in some cases very challenging) analysis. We therefore turn to direct numerical simulations of the generalized Lotka-Volterra equation , and we show that the results agree with and extend the theory presented in the previous section. The details on the numerical scheme implemented for the simulation can be found in Appendix . These simulations are challenging as we are interested in considering both a large number of species and a large number of patches. Moreover, lowering the temperature results in a strong slowdown of the dynamics (Appendix ), leading to additional computational costs. The slowdown of the dynamics is much stronger in the presence of heterogeneity in the interactions than with zero or constant ones.

A. Finite number of species and finite number of patches

Generically, for moderate system sizes ( $𝑆 < 100$ and $𝐿 < 100$ ) we find strong fluctuations due to the quenched disorder in the interaction matrix, and quantitative finite-size effects compared to the asymptotic $𝑆, 𝐿 \to \infty$ solution, in particular for $𝜌 = 1$ [for $𝜌 = 0$ each patch is characterized by an independent realization of the interaction matrix, thus leading to a faster (self-averaging) convergence of the system to its disorder average]. For larger values of $𝑆$ and $𝐿$ , e.g., $𝑆 = 200, 𝐿 = 400$ , fluctuations and finite-size effects are limited and one finds results that are both qualitative and quantitative in agreement with the analytical solution.

In Fig. we show the behavior of the average abundance as a function of the diffusion constant for two different values of the temperature, starting from two different initial conditions. To probe the existence of hysteresis, and therefore a discontinuous transition and metastability, we numerically simulate systems with different initial conditions. For the green curves, the initial abundances were uniformly sampled between 0 and 1, for the red curves between 0 and 0.1. The former should therefore be more prone to evolve toward the self-sustained solution, if it exists, whereas the latter is prone to evolve toward the “all-extinct” solution.

FIG. 8.

Average abundance $‾ ⟨ 𝑁 ⟩$ as a function of the diffusion constant $𝐷$ for (a) $𝑇 = 0.18$ and (b) $𝑇$ = 0.8. Green and red dots indicate the initial conditions of order 1 and of order 0.1. The dashed line indicates the analytical prediction for the critical value of the diffusion constant for the continuous transition. $𝜇 = 1, 𝜎 = 0.5, 𝑆 = 200, 𝐿 = 400, 𝑡 max = 500$ (left) and 200 (right).

We find that indeed at higher temperatures, $𝑇 = 0.8$ , in agreement with the analytics and the phase diagram in Fig. , the final abundances vary continuously when varying the diffusion constant, and they converge to the same value regardless of the initial condition. The value of $𝐷$ at which the final abundances significantly depart from zero quantitatively matches the analytical result for the critical value of the diffusion constant at the continuous transition.

Instead, at $𝑇 = 0.18$ the final abundances show a strong dependence on the initial condition in an extended interval of diffusion strengths; for a given initial condition, the final abundance exhibits a very abrupt change . Interestingly, the dynamics strongly slows down in this regime, in particular for the decay of the abundances from large initial conditions. In fact, this process occurs via the rare extinctions of species that are asymptotically not able to self-sustain but can persist for very long times, especially in this regime in which demographic fluctuations are weak. The strong dependence on the initial conditions cannot be explained just by the slowdown of the dynamics because the abundances with different initial conditions evolve in opposite directions (see Fig. and Appendix ).

The heterogeneity in the interaction network is essential to allow the ecosystem to self-sustain below the single DP critical point: indeed if we consider the same parameters but take $𝜎 = 0$ , all species go extinct below $𝐷 0 (𝑇)$ , and there is no strong dependence on the initial conditions (Appendix ).

B. Asymmetric interactions and partial correlation between patches

We are now interested in focusing on cases in which the interactions between species are not fully symmetric, and the interaction matrices are partially correlated between patches, i.e., $0 < 𝜌 < 1$ .

As we have already discussed, we have analytically established that a very small asymmetry is not a singular perturbation. Thus, our results should qualitatively hold also for a finite, at least not too large, asymmetry.

To confirm this finding and study intermediate values of $𝜌$ (besides $𝜌 = 0, 1$ considered analytically), we performed simulations with spatial heterogeneity $𝜌 = 0.9$ and asymmetry in the interactions $𝛾 = 0.9$ , and as before for $𝐿 = 400, 𝑆 = 200$ . Also in this case at $𝑇 = 0.8$ we find a continuous transition and no strong dependence on the initial conditions, while at $𝑇 = 0.18$ we find a discontinuous transition and a hysteresis region (Fig. ) . Although the curves quantitatively change with respect to their $𝛾 = 𝜌 = 1$ counterparts, as expected, the results and in particular the existence of a discontinuous transition do remain qualitatively unaltered.

FIG. 9.

Average abundance $‾ ⟨ 𝑁 ⟩$ as a function of the diffusion constant $𝐷$ for (a) $𝑇 = 0.18$ and (b) $𝑇$ = 0.8 with some spatial heterogeneity ( $𝜌 = 0.9$ ) and some asymmetry in the interactions ( $𝛾 = 0.9$ ). Green and red lines indicate the initial conditions of order 1 and of order 0.1, lighter dots show the average abundance at intermediate times ( $50 %$ and $75 %$ of $𝑡 max$ ). $𝜇 = 1, 𝜎 = 0.5, 𝑆 = 200, 𝐿 = 400, 𝑡 max = 500$ (left) and 200 (right).

In conclusion, combining all these numerical tests, we conclude that the scenario obtained from the analytical solution is robust and holds broadly. We will come back to this point in the Conclusion to suggest other extensions and tests.

VI. PRECURSOR OF THE INSTABILITY TOWARD EXTINCTION

In the previous section, we have shown that dispersal can rescue complex and large ecosystems from extinction due to demographic noise. Depending on the strength of the demographic noise, the transition from the self-sustained to the extinct phase can be either continuous or discontinuous. The latter takes place for low demographic noise and low dispersal. In this regime, we have found that the transition is accompanied by a metastable regime and hysteresis. Such a transition is what is called in ecology, and in environmental and social sciences, a tipping point or regime shift , and in physics it is called a spinodal. Tipping points are often catastrophic events, as the abrupt rapid shifts almost always lead to negative consequences and a less favorable state of the system. Our case is no exception, as the system's transition is from a self-sustained state with high diversity to one in which all species are extinct. As was done for several other tipping points , it is therefore important to find early signs or precursors that can allow one to detect the closeness of the system to the tipping point before the catastrophic shift actually takes place.

In our case, following intuition that comes from the physics of spinodal points, we focus on responses to perturbations as a probe of closeness to the tipping point. We can show analytically (see Appendix ) that the instability of the self-sustained state is accompanied by a diverging response to perturbations. This phenomenon is strongly linked to the saddle-node bifurcation of the mean-field equations that governs the transition.

In particular, we have studied the change of the average abundance due to a change in the carrying capacity. Such a response, which can be measured in controlled laboratory experiments, does diverge approaching the discontinuous transition; see Fig. for the $𝜌 = 0$ case. A similar behavior is expected for generic values of $𝜌$ . This probe can therefore be used as an early warning signal of the proximity to the tipping point of the self-sustained phase. In natural ecosystems, where measuring responses to perturbation can be challenging, one could instead monitor the long-term fluctuations of average abundance due to environmental noise affecting the carrying capacity on a long time. This would be a proxy for the response proposed above (it is important to focus on long times as all the processes at play are slow).

FIG. 10.

Average abundance and its response to a perturbation of the carrying capacity $𝑘$ at $𝑇 = 0.153$ for $𝜌 = 0$ approaching the instability of the self-sustained phase. $𝜇 = 1, 𝜎 = 0.5$ .

VII. CONCLUSIONS

We uncovered a rich phase diagram for many-species Lotka-Volterra metacommunities subject to heterogeneous symmetric interactions, demographic noise, and diffusion. If the demographic fluctuations are too strong, they drive all species to extinctions, but when the diffusion constant is large enough these extinctions can be compensated by recolonizations from neighboring sites, and the ecosystem is able to self-sustain at finite abundance and diversity. The system exhibits a phase transition between an extinction and a survival phase. The transition can be either continuous or discontinuous, depending on whether the behavior of the system is dominated by the demographic fluctuations or the heterogeneous interaction network.

When the demographic fluctuations are strong, the transition is continuous and interactions play a secondary role. In fact, the transition is completely analogous to what one would have in the absence of the interactions (even the critical values of the diffusion constant coincide). This is because when the abundances tend to zero, the interactions become subdominant and the system falls in the standard Directed Percolation universality class.

The situation is drastically different at lower demographic noise. In this case, the transition becomes discontinuous and the system exhibits novel features, which are a signature of the complexity of the ecosystem and the major role played by the interactions. There is an extended range of parameters in which without interactions, i.e., for single species, the system would be driven to extinction, but the metacommunity is instead able to self-sustain at finite abundances. This is possible because strongly competing species are eliminated from the community, while surviving species cooperate to self-sustain in such harsh conditions. For small demographic noise and lowering the diffusion constant, the ecosystem reaches a tipping point at which all surviving species go extinct; close to this point the ecosystem is subject to collapses upon small perturbations, and its dynamics exhibits hysteresis. We therefore find that mutualism naturally emerges from an (on average) competitive pool of species when conditions become harsher. This has a double effect: it allows the ecosystem to survive in conditions in which all species in isolation would go extinct, but it also makes it fragile to perturbations. In this regime, it is not possible to predict the vicinity of the catastrophic shift of the ecosystem by looking at the average abundance. As an early warning sign, we propose to monitor the response of the system to perturbations. We have shown that this is a suitable probe, as it diverges approaching the discontinuous transition.

We confirm and complement our analytical approach with numerical simulations, which show that our results are quite robust to modifications of the model, in particular to the introduction of a small asymmetry in the interactions, to various degrees of correlation of the interaction network between different spatial locations, and for a system with a finite number of species and patches.

There are several directions that warrant future investigations. We focused on a fully connected spatial system, which provides a mean-field analysis for generic spatial lattices. On the other hand, our DMFT treatment of the interactions is directly generalizable to any other spatial network, including finite-dimensional ones. It would be very interesting to study cases in which the patches are located in a finite-dimensional lattice or on random structures. In particular, it would be interesting to find out (1) whether the discontinuous transition is also present in this case or if finite-dimensional fluctuations destroy the metastable region, and (2) whether the continuous transition can still be described in terms of directed percolation, or if interactions, although secondary, can alter its universality class. It would also be worth analyzing stronger asymmetries in the interactions, e.g., lowering the value of $𝛾$ . We expect that a significant positive correlation between reciprocal interactions is needed to induce metastability. This ensures that species that interact more competitively are also more negatively affected by the interactions with the rest of the ecosystem and hence go extinct, thus leading to mutualism for the surviving species.

Finally, the species-rich LV model with heterogeneous and strong interactions displays multiple equilibria and chaotic dynamics . The possibility of different patches to converge to different stationary states could strongly modify the behavior of the system, in particular allowing the system to experience higher values of global diversity, possibly violating May's bound .

Note added. During the preparation of this manuscript, we became aware of the work of Denk and Hallatscheck on tipping points in mutualistic Lotka-Volterra communities . Their results are complementary and agree with ours.

ACKNOWLEDGMENTS

We would like to thank J. Denk and O. Hallatscheck for sharing their results and constructive interactions. We also thank Joseph Baron, M. Barbier, J. F. Arnoldi, and L. F. Cugliandolo for stimulating discussions. This work was supported by the Simons Foundation Grant on Cracking the Glass Problem (No. 454935) (G.B.).

Appendices

APPENDIX A: DMFT DERIVATION

Here we outline the derivation, adapted from Ref. , of the Dynamical Mean Field Theory for our system, for a generic value of the spatial correlation of the interactions $𝜌$ .

We consider $𝑆$ species, indexed by $𝑖 = 1, \dots, 𝑆$ , and their Lotka-Volterra dynamics,

̇ 𝑁 𝑖, 𝑢 = 𝑁 𝑖, 𝑢 (1 - 𝑁 𝑖, 𝑢 - \sum 𝑗 𝛼 𝑢 𝑖 𝑗 𝑁 𝑗, 𝑢 + 𝜁 𝑖, 𝑢) + 𝐷 (1 𝑐 \sum 𝑣 \in 𝜕 𝑢 𝑁 𝑖, 𝑣 - 𝑁 𝑖, 𝑢) + 𝜂 𝑢 𝑖 (𝑡) \sqrt 𝑁 𝑖, 𝑢 + 𝜆 𝑖, 𝑢,

(A1)

to which we have added a perturbation to the carrying capacity $𝜁 𝑖, 𝑢$ and an external immigration $𝜆 𝑖, 𝑢$ , which will be taken to be zero at the end of the computation. These equations [for a given value of the $𝜂 𝑢 𝑖 (𝑡)$ ] define the trajectories $𝑁 𝑖, 𝑢 (𝑡)$ . We add a new species, $𝑖 = 0$ , to the system, and we draw its interactions and initial conditions independently from the rest of the system and with the same statistics. At large $𝑆$ , thanks to the scaling of the interactions, the presence of a new species is a small perturbation to the system, so that the trajectories of the other $𝑆$ species will only be slightly modified. We consider their linear response:

𝛿 𝑁 𝑖, 𝑢 (𝑡) = \sum 𝑣 \in 𝜕 𝑢, 𝑗 \int 𝑡 0 𝛿 𝑁 𝑖, 𝑢 (𝑡) 𝛿 𝜁 𝑗, 𝑣 (𝑡') [- 𝛼 𝑣 𝑗 0 𝑁 0, 𝑣 (𝑡')] 𝑑 𝑡' = \sum 𝑣 \in 𝜕 𝑢, 𝑖 \int 𝑡 0 𝑅 𝑢, 𝑣 𝑖, 𝑗 (𝑡, 𝑡') [- 𝛼 𝑣 𝑗 0 𝑁 0, 𝑣 (𝑡')] 𝑑 𝑡' .

(A2)

We have introduced the response function $𝑅 𝑢, 𝑣 𝑖, 𝑗 (𝑡, 𝑡')$ of the abundance of species $𝑖$ in patch $𝑢$ at time $𝑡$ to a variation in the carrying capacity of species $𝑗$ in patch $𝑣$ at time $𝑡'$ .

The dynamics of species 0 will depend on these new trajectories:

̇ 𝑁 0, 𝑢 = 𝑁 0, 𝑢 (1 - 𝑁 0, 𝑢 - \sum 𝑖 𝛼 𝑢 0 𝑖 (𝑁 0 𝑖, 𝑢 + 𝛿 𝑁 𝑖, 𝑢)) + 𝐷 (1 𝑐 \sum 𝑣 \in 𝜕 𝑢 𝑁 0, 𝑣 - 𝑁 0, 𝑢) + 𝜂 0, 𝑢 (𝑡) \sqrt 𝑁 0, 𝑢 .

(A3)

Because the correlations between interaction coefficients in any two patches are the same, these Gaussian variables can generically be decomposed into a common random contribution, identical in all patches and proportional to the correlation $𝜌$ , and one independent in different patches, proportional to $\sqrt 1 - 𝜌 2$ . We thus introduce the matrix $𝑎 𝑖, 𝑗$ and $𝑎 𝑢 𝑖, 𝑗$ such that $𝛼 𝑢 𝑖, 𝑗 = 𝜇 / 𝑆 + 𝜎 (𝜌 𝑎 𝑖, 𝑗 + \sqrt 1 - 𝜌 2 𝑎 𝑢 𝑖, 𝑗), 𝔼 [𝑎 𝑖, 𝑗] = 𝔼 [𝑎 𝑢 𝑖, 𝑗] = 0, 𝔼 [𝑎 2 𝑖, 𝑗] = 𝔼 [𝑎 𝑢 𝑖, 𝑗 2] = 1 / 𝑆, 𝔼 [𝑎 𝑖, 𝑗 𝑎 𝑗, 𝑖] = 𝔼 [𝑎 𝑢 𝑖, 𝑗 𝑎 𝑢 𝑗, 𝑖] = 𝛾 / 𝑆$ , and all other correlations are 0. We can rewrite the interaction term as

- \sum 𝑖 𝛼 0 𝑖 (𝑁 0 𝑖, 𝑢 + 𝛿 𝑁 𝑖, 𝑢) = - \sum 𝑖 [𝜇 / 𝑆 + 𝜎 (𝜌 𝑎 0 𝑖 + \sqrt 1 - 𝜌 2 𝑎 𝑢 0 𝑖)] 𝑁 0 𝑖, 𝑢 + \sum 𝑖, 𝑗 [𝜇 / 𝑆 + 𝜎 (𝜌 𝑎 𝑖 0 + \sqrt 1 - 𝜌 2 𝑎 𝑢 𝑖 0)] \times [𝜇 / 𝑆 + 𝜎 (𝜌 𝑎 0 𝑗 + \sqrt 1 - 𝜌 2 𝑎 𝑢 0 𝑗)] \sum 𝑣 \in 𝜕 𝑢 \int 𝑡 0 𝑅 𝑢, 𝑣 𝑖, 𝑗 (𝑡, 𝑡') 𝑁 0, 𝑣 (𝑡') 𝑑 𝑡' .

(A4)

We want to describe its statistical properties in the limit $𝑆 \to \infty$ . The response function $𝑅 𝑢, 𝑣 𝑖, 𝑗 (𝑡, 𝑡')$ is defined on the unperturbed trajectories, and is therefore uncorrelated from the interactions coefficients with species 0. $𝑅 𝑢, 𝑣 𝑖, 𝑗 (𝑡, 𝑡') \sim 1 / \sqrt 𝑆$ for $𝑖 \neq 𝑗$ , so that the off-diagonal terms can be neglected. Thanks to the central limit theorem, $\sum 𝑗 𝑎 0 𝑗 𝑎 𝑗 0 𝑅 𝑢, 𝑣 𝑗, 𝑗 (𝑡, 𝑡')$ will converge to its average:

\sum 𝑗 𝑎 0 𝑗 𝑎 𝑗 0 𝑅 𝑢, 𝑣 𝑗, 𝑗 (𝑡, 𝑡') \to 𝑆 𝔼 [𝑎 0 𝑗 𝑎 𝑗 0] 𝔼 [𝑅 𝑢, 𝑣 𝑗, 𝑗 (𝑡, 𝑡')] = 𝛾 𝔼 [𝑅 𝑢, 𝑣 𝑗, 𝑗 (𝑡, 𝑡')] .

(A5)

By similarly evaluating all terms in , we obtain

- \sum 𝑗 𝛼 0 𝑗 (𝑁 0 𝑗, 𝑢 + 𝛿 𝑁 𝑗, 𝑢) \to - 𝜇 𝔼 [𝑁 0 𝑗, 𝑢] - 𝜎 𝜌 ̃ 𝜉 𝑢 (𝑡) - 𝜎 \sqrt 1 - 𝜌 2 ˜ 𝜓 𝑢 (𝑡) + 𝜎 2 𝜌 2 𝛾 \sum 𝑣 \in 𝜕 𝑢 \int 𝑡 0 𝔼 [𝑅 𝑢, 𝑣 𝑗, 𝑗 (𝑡, 𝑡')] 𝑁 0, 𝑣 (𝑡') 𝑑 𝑡' + 𝜎 2 (1 - 𝜌 2) 𝛾 \int 𝑡 0 𝔼 [𝑅 𝑢, 𝑢 𝑗, 𝑗 (𝑡, 𝑡')] 𝑁 0, 𝑢 (𝑡') 𝑑 𝑡',

(A6)

where $̃ 𝜉 𝑢 (𝑡)$ and $˜ 𝜓 𝑢 (𝑡)$ are Gaussian fields with zero mean and covariance $𝔼 [̃ 𝜉 𝑢 (𝑡) ̃ 𝜉 𝑣 (𝑡')] = 𝔼 [𝑁 0 𝑗, 𝑢 (𝑡) 𝑁 0 𝑗, 𝑣 (𝑡')], 𝔼 [˜ 𝜓 𝑢 (𝑡) ˜ 𝜓 𝑣 (𝑡')] = 𝛿 𝑢 𝑣 𝔼 [𝑁 0 𝑗, 𝑢 (𝑡) 𝑁 0 𝑗, 𝑢 (𝑡')]$ . Note that $̃ 𝜉 𝑢$ and the first integral of derive from the component of the interactions constant across patches, $𝑎 𝑖 𝑗$ , as we can see from the $𝜌$ -dependent prefactors, and they therefore couple different patches. $˜ 𝜓 𝑢$ and the second integral of derive instead from the component of the interactions independent across patches, $𝑎 𝑢 𝑖 𝑗$ , and therefore they represent diagonal correlations and responses. Plugging this expression in the dynamical equation for species 0, we obtain

̇ 𝑁 0, 𝑢 = 𝑁 0, 𝑢 (1 - 𝑁 0, 𝑢 - 𝜇 𝔼 [𝑁 0 𝑗, 𝑢] - 𝜎 𝜌 ̃ 𝜉 𝑢 (𝑡) - 𝜎 \sqrt 1 - 𝜌 2 ˜ 𝜓 𝑢 (𝑡) + 𝜎 2 𝜌 2 𝛾 \sum 𝑣 \in 𝜕 𝑢 \int 𝑡 0 𝔼 [𝑅 𝑢, 𝑣 𝑗, 𝑗 (𝑡, 𝑡')] 𝑁 0, 𝑣 (𝑡') 𝑑 𝑡' + 𝜎 2 (1 - 𝜌 2) 𝛾 \int 𝑡 0 𝔼 [𝑅 𝑢, 𝑢 𝑗, 𝑗 (𝑡, 𝑡')] 𝑁 0, 𝑢 (𝑡') 𝑑 𝑡') + 𝐷 (1 𝑐 \sum 𝑣 \in 𝜕 𝑢 𝑁 0, 𝑣 - 𝑁 0, 𝑢) + 𝜂 0, 𝑢 (𝑡) \sqrt 𝑁 0, 𝑢 .

(A7)

Species 0 is statistically equivalent to all the others, therefore we can replace the averages over the $𝑆$ original species with averages with respect to this new dynamics for a single species, obtaining some self-consistent equations:

̇ 𝑁 𝑢 = 𝑁 𝑢 (1 - 𝑁 𝑢 - 𝜇 ℎ 𝑢 - 𝜎 𝜌 ̃ 𝜉 𝑢 (𝑡) - 𝜎 \sqrt 1 - 𝜌 2 ˜ 𝜓 𝑢 (𝑡) + 𝜎 2 𝜌 2 𝛾 \int 𝑡 0 \sum 𝑣 \in 𝜕 𝑢 𝑅 𝑢 𝑣 (𝑡, 𝑠) 𝑁 𝑣 (𝑠) 𝑑 𝑠

(A8)

+ 𝜎 2 (1 - 𝜌 2) 𝛾 \int 𝑡 0 𝑅 𝑢 𝑢 (𝑡, 𝑠) 𝑁 𝑢 (𝑠) 𝑑 𝑠) + 𝐷 (1 𝑐 \sum 𝑣 \in 𝜕 𝑢 𝑁 𝑣 - 𝑁 𝑢) + 𝜂 𝑢 (𝑡) \sqrt 𝑁 𝑢,

(A9)

⟨ ̃ 𝜉 𝑢 (𝑡) ̃ 𝜉 𝑣 (𝑠) ⟩ = 𝐶 𝑢 𝑣 (𝑡 - 𝑠) = 𝔼 [𝑁 𝑢 (𝑡) 𝑁 𝑣 (𝑠)],

(A10)

⟨ ˜ 𝜓 𝑢 (𝑡) ˜ 𝜓 𝑣 (𝑠) ⟩ = 𝛿 𝑢 𝑣 𝐶 𝑢 𝑢 (𝑡 - 𝑠),

(A11)

𝑅 𝑢 𝑣 (𝑡, 𝑠) = 𝔼 [𝛿 𝑁 𝑢 (𝑡) 𝛿 𝜁 𝑣 (𝑠) | 𝜁 = 0],

(A12)

ℎ 𝑢 = 𝔼 [𝑁 𝑢] .

(A13)

Since species have been effectively decoupled, we can suppress the species index.

In the single equilibrium phase, we expect the process to reach a time translation invariant regime, in which the one-time averages are time-independent, and two-times observables only depend on the times difference. This was shown in for a single community with demographic noise and fixed immigration, and it is known to be the case for Directed Percolation and in a many-species metacommunity with constant interactions . It is also confirmed by our numerical results, which show a quick relaxation of one-time observables to an asymptotic value (see Appendix ), at least away from phase transitions. Since the autocorrelation of the abundance of one species does not tend to zero at large times, we can decompose $𝜉 𝑢$ and $𝜓 𝑢$ into a constant and a fluctuating component:

̃ 𝜉 𝑢 (𝑡) = ̂ 𝜉 𝑢 + 𝜉 𝑢 (𝑡),

(A14)

˜ 𝜓 𝑢 (𝑡) = ˆ 𝜓 𝑢 + 𝜓 𝑢 (𝑡),

(A15)

where $̂ 𝜉 𝑢$ and $ˆ 𝜓 𝑢$ are (time-independent) Gaussian variables with zero mean and correlations $l i m 𝜏 \to \infty 𝐶 𝑢 𝑣 (𝑡, 𝑡 + 𝜏) = 𝐶 \infty 𝑢 𝑣$ and $𝛿 𝑢 𝑣 𝐶 \infty 𝑢 𝑢$ , and the autocorrelation of $𝜉 𝑢$ and $𝜓 𝑢$ goes to zero at long times. Averaging over $𝜉 𝑢, 𝜓 𝑢$ , and $𝜂$ at fixed $̂ 𝜉 𝑢$ and $ˆ 𝜓 𝑢$ corresponds to performing a time-average for one species in one patch; averaging also over $ˆ 𝜓 𝑢$ and $̂ 𝜉 𝑢$ corresponds to averaging over patches and species. In this sense, $̂ 𝜉 𝑢$ and $ˆ 𝜓 𝑢$ play the role of the quenched disorder, which was previously represented by the interaction matrix $𝛼 𝑢 𝑖 𝑗$ . We will refer to the average over $𝜉$ and $𝜂$ at fixed $̂ 𝜉 𝑢$ and $ˆ 𝜓 𝑢$ as thermal average and indicate it with brackets, and to the average over $̂ 𝜉 𝑢$ and $ˆ 𝜓 𝑢$ as disorder average and indicate it with an overline.

While the derivation is so far valid for any spatial network, we will now restrict ourselves to a fully connected network, in which the empirical average over neighbors can be replaced by its thermal average. In the large- $𝐿$ limit, the connected correlation over thermal fluctuations between $𝑁 𝑢$ and $𝑁 𝑣$ is subdominant, so that $𝜉 𝑢$ and $𝜉 𝑣$ become independent. A perturbation in patch $𝑣$ influences the abundance in patch $𝑢$ through the diffusion term, which in a fully connected network is of order $1 / 𝐿$ , therefore $𝑅 𝑢 𝑣$ for $𝑢 \neq 𝑣$ scales as $1 / 𝐿$ , whereas $𝑅 𝑢 𝑢$ is of order 1. Since all patches are equivalent, the elements of the $𝑅 𝑢 𝑣$ matrix can only take two values:

𝑅 𝑢 𝑢 = 𝑅 𝑑

(A16)

𝑅 𝑢 𝑣 = 𝑅 0 / 𝐿, 𝑢 \neq 𝑣 .

(A17)

The same is true for $𝐶 𝑢 𝑣$ :

𝐶 \infty 𝑢 𝑢 = 𝐶 \infty 𝑑,

(A18)

𝐶 \infty 𝑢 𝑣 = 𝐶 \infty 0, 𝑢 \neq 𝑣 .

(A19)

We separate $̂ 𝜉 𝑢$ in a patch-independent and a patch-dependent part: $̂ 𝜉 𝑢 = 𝑧 \sqrt 𝐶 \infty 0 + 𝑤 𝑢 \sqrt 𝐶 \infty 𝑑 - 𝐶 \infty 0$ . We call patch disorder average the average over $𝑤 𝑢$ and $ˆ 𝜓 𝑢$ , and species disorder average the average over $𝑧$ . $1 𝐿 \sum 𝑣 𝑁 𝑣$ concentrates around its average over thermal fluctuations and patch disorder $𝑁 *$ , which will be a function of the static Gaussian field $𝑧$ . Substituting in the dynamical equation and using time translational

invariance, we obtain

̇ 𝑁 = 𝑁 [𝑘 - 𝑁 - 𝜇 ℎ - 𝜎 (𝜌 \sqrt 𝐶 \infty 0 𝑧 + 𝜌 \sqrt 𝐶 \infty 𝑑 - 𝐶 \infty 0 𝑤 + \sqrt 1 - 𝜌 2 \sqrt 𝐶 \infty 𝑑 ˆ 𝜓 + 𝜌 𝜉 + \sqrt 1 - 𝜌 2 𝜓)] + 𝑁 𝜎 2 𝛾 (𝜌 2 \int 𝑡 0 𝑅 𝑑 (𝑡 - 𝑠) 𝑁 (𝑠) 𝑑 𝑠 + 𝜌 2 \int 𝑡 0 𝑅 0 (𝑡 - 𝑠) 𝑁 * (𝑠) 𝑑 𝑠 + (1 - 𝜌 2) \int 𝑡 0 𝑅 𝑑 (𝑡 - 𝑠) 𝑁 (𝑠) 𝑑 𝑠) + 𝐷 (𝑁 * - 𝑁) + 𝜂 (𝑡) \sqrt 𝑁 = 𝑁 [𝑘 - 𝑁 - 𝜇 ℎ - 𝜎 (𝜌 \sqrt 𝐶 \infty 0 𝑧 + \sqrt 𝐶 \infty 𝑑 - 𝜌 2 𝐶 \infty 0 𝑤 + 𝜉)] + 𝑁 𝜎 2 𝛾 (\int 𝑡 0 𝑅 𝑑 (𝑡 - 𝑠) 𝑁 (𝑠) 𝑑 𝑠 + 𝜌 2 𝑅 int 0 𝑁 *) + 𝐷 (𝑁 * - 𝑁) + 𝜂 (𝑡) \sqrt 𝑁,

(A20)

where we have summed the random variables that had the same behavior of the correlations ( $𝑤$ and $ˆ 𝜓, 𝜉$ and $𝜓$ ), and

𝑅 int 0 = \int \infty 0 𝑑 𝜏 𝑅 0 (𝜏) .

(A21)

The equations simplify in the extreme cases $𝜌 = 1$ and 0, because only one of the components of the static part of the disorder is present, either $𝑤$ or $𝑧$ . For $𝜌 = 1, 𝐶 \infty 𝑑 = 𝐶 \infty 0$ . For $𝜌 = 0, 𝑁 *$ coincides with $ℎ$ , so that we have one less self-consistent equation, and $𝑅 0$ is not present; these two facts greatly simplify the numerical solution of the equations.

APPENDIX B: STATIONARY PROBABILITY DISTRIBUTION IN THE SYMMETRIC CASE

In the case of symmetric interactions ( $𝛾 = 1$ ), in the single equilibrium phase, the system relaxes to equilibrium and it verifies the Fluctuation-Dissipation Theorem (FDT) :

𝑅 𝑑 (𝜏) = - 1 𝑇 𝑑 𝐶 (𝜏) 𝑑 𝜏 .

(B1)

We can integrate by parts the term with the memory kernel:

\int 𝑡 0 𝑅 𝑑 (𝑡 - 𝑠) 𝑁 (𝑠) 𝑑 𝑠 = 1 𝑇 \int 𝑡 0 𝑑 𝐶 𝑑 (𝑡 - 𝑠) 𝑑 𝑠 𝑁 (𝑠) 𝑑 𝑠

(B2)

= 1 𝑇 (𝐶 0 𝑑 𝑁 (𝑡) - 𝐶 (𝑡) 𝑁 (0) - \int 𝑡 0 𝐶 𝑑 (𝑡 - 𝑠) ̇ 𝑁 (𝑠) 𝑑 𝑠)

(B3)

= 1 𝑇 ((𝐶 0 𝑑 - 𝐶 \infty 𝑑) 𝑁 (𝑡) - \int 𝑡 0 [𝐶 𝑑 (𝑡 - 𝑠) - 𝐶 \infty 𝑑] ̇ 𝑁 (𝑠) 𝑑 𝑠) .

(B4)

We have obtained an additional quadratic term in $𝑁 (𝑡)$ , and a friction term. The friction term and the noise $𝜉$ describe the coupling of the system to an effective colored bath at temperature $𝑇$ , which replaces the coupling of one species to all the others.

Using the Martin–Siggia–Rose–De Dominicis–Janssen (MSRDJ) formalism, we can show that the stationary probability distribution associated with the stochastic differential equation

̇ 𝑁 = 𝑁 (𝑘 - 𝐷 (1 - 𝜌 2 𝜎 2 𝑅 int 0 𝑁 *) - 𝜇 ℎ - 𝜌 𝜎 \sqrt 𝐶 \infty 0 𝑧 - 𝜎 \sqrt 𝐶 \infty 𝑑 - 𝜌 2 𝐶 \infty 0 𝑤 - 𝜎 𝜉 (𝑡)) - 𝑁 2 (1 - 𝜎 2 𝑇 (𝐶 0 𝑑 - 𝐶 \infty 𝑑)) + 𝑁 𝜎 2 𝑇 \int 𝑡 0 [𝐶 𝑑 (𝑡 - 𝑠) - 𝐶 \infty 𝑑] ̇ 𝑁 (𝑠) 𝑑 𝑠 + 𝐷 𝑁 * + 𝜂 (𝑡) \sqrt 𝑁

(B5)

is the Boltzmann distribution with the effective Hamiltonian:

𝐻 eff = (1 - 𝜎 2 𝑇 (𝐶 0 𝑑 - 𝐶 \infty 𝑑)) 𝑁 2 / 2 - (𝑘 - 𝐷 (1 - 𝜌 2 𝜎 2 𝑁 * 𝑅 int 0) - 𝜇 ℎ - 𝜌 𝜎 \sqrt 𝐶 \infty 0 𝑧 - 𝜎 \sqrt 𝐶 \infty 𝑑 - 𝜌 2 𝐶 \infty 0 𝑤 + 𝜁) 𝑁 + (𝑇 - 𝐷 𝑁 * + 𝜆) l n 𝑁,

(B6)

where we have reintroduced the perturbations $𝜁$ and $𝜆$ . To show that this is the correct equilibrium distribution, we need to verify that, with this as an initial condition, time reversal is a symmetry of the associated MSRDJ action. We will do it, following Ref. , for a simplified dynamics, which contains all the crucial ingredients:

̇ 𝑁 = 𝑁 (1 - 𝑁 - 𝜎 𝜉 (𝑡) - 𝜎 2 \int 𝑡 0 𝜈 (𝑡, 𝑠) ̇ 𝑁 (𝑠) 𝑑 𝑠) + 𝜂 (𝑡) \sqrt 𝑁 + 𝜆,

(B7)

⟨ 𝜉 (𝑡) 𝜉 (𝑠) ⟩ = 𝑇 𝜈 (𝑡 - 𝑠),

(B8)

⟨ 𝜂 (𝑡) 𝜂 (𝑠) ⟩ = 2 𝑇 𝛿 (𝑡 - 𝑠) .

(B9)

Its equilibrium distribution is given by

𝑃 eq (𝑁) = 𝑒 - 𝛽 𝐻 𝑍,

(B10)

𝐻 = 𝑁 2 / 2 - 𝑁 + (𝑇 - 𝜆) l o g 𝑁,

(B11)

where $𝛽 = 1 / 𝑇$ , the inverse temperature. The white noise should be interpreted according to Ito's discretization. It is convenient to convert it to Stratonovich's discretization, which is left invariant by time reversal. The multiplicative nature of the noise makes the two discretizations not equivalent: we then need to introduce an additional drift term as follows:

𝜂 \sqrt 𝑁 \to 𝜂 \sqrt 𝑁 - 12 \sqrt 2 𝑇 2 \sqrt 𝑁 \sqrt 2 𝑇 𝑁 = 𝜂 \sqrt 𝑁 - 𝑇 2 .

(B12)

The MSRDJ action can be written in terms of a deterministic and a dissipative part ,

𝑆 [𝑁, ˆ 𝑁] = 𝑆 det [𝑁, ˆ 𝑁] + 𝑆 diss [𝑁, ˆ 𝑁],

(B13)

𝑆 det [𝑁, ˆ 𝑁] = l o g 𝑃 eq (𝑁 (- 𝑇)) + \int 𝑇 - 𝑇 𝑑 𝑢 (𝑖 ˆ 𝑁 [𝑁 (1 - 𝑁) + 𝜆 - 𝑇 / 2 - 𝑇 / 2] + 𝑁 - 1 / 2),

(B14)

𝑆 diss [𝑁, ˆ 𝑁] = \int 𝑢 𝑖 ˆ 𝑁 𝑢 \int 𝑣 [𝛿 (𝑢 - 𝑣) + 𝜈 (𝑢 - 𝑣) 𝜃 (𝑢 - 𝑣) 𝑁 𝑢] (𝑖 𝑇 ˆ 𝑁 𝑣 𝑁 𝑣 - ̇ 𝑁 𝑣) .

(B15)

The time reversal transformation for the two fields is given by

𝑁 (𝑡) ⟶ 𝑁 𝑅 (𝑡) = 𝑁 (- 𝑡),

(B16)

𝑖 ˆ 𝑁 (𝑡) ⟶ 𝑖 ˆ 𝑁 𝑅 (𝑡) = 𝑖 ˆ 𝑁 (- 𝑡) + 1 𝑇 𝑁 (- 𝑡) 𝜕 𝜕 𝑡 𝑁 (- 𝑡) .

(B17)

The deterministic and dissipative part of the action are independently invariant under this transformation:

𝑆 det [𝑁 𝑅, ˆ 𝑁 𝑅] = - l o g 𝑍 - 𝛽 𝐻 (𝑁 (𝑇)) + \int 𝑢 [(𝑖 ˆ 𝑁 - 𝑢 + 1 𝑇 𝑁 - 𝑢 𝜕 𝜕 𝑢 𝑁 - 𝑢) (𝑁 - 𝑢 (1 - 𝑁 - 𝑢) + 𝜆 - 𝑇) + 𝑁 - 𝑢 - 1 / 2] = - l o g 𝑍 - 1 𝑇 [𝑁 2 𝑇 / 2 - 𝑁 𝑇 + (𝑇 - 𝜆) l n 𝑁 𝑇] + 1 𝑇 \int 𝑢 𝜕 𝜕 𝑢 (𝑁 2 𝑢 / 2 - 𝑁 𝑢 + (𝑇 - 𝜆) l n 𝑁 𝑢) + \int 𝑢 (𝑖 ˆ 𝑁 𝑢 [𝑁 𝑢 (1 - 𝑁 𝑢) + 𝜆 - 𝑇] + 𝑁 𝑢 - 1 / 2) = - l o g 𝑍 - 𝛽 𝐻 (𝑁 (- 𝑇)) + \int 𝑢 (𝑖 ˆ 𝑁 𝑢 [𝑁 𝑢 (1 - 𝑁 𝑢) + 𝜆 - 𝑇] + 𝑁 𝑢 - 1 / 2) = 𝑆 det [𝑁, ˆ 𝑁],

(B18)

𝑆 diss [𝑁 𝑅, ˆ 𝑁 𝑅] = \int 𝑢 (𝑖 ˆ 𝑁 - 𝑢 + 1 𝑇 𝑁 - 𝑢 𝜕 𝜕 𝑢 𝑁 - 𝑢) \int 𝑣 (𝛿 𝑢 - 𝑣 + 𝜈 𝑢 - 𝑣 𝜃 𝑢 - 𝑣 𝑁 𝑢) 𝑖 𝑇 ˆ 𝑁 - 𝑣 𝑁 - 𝑣 = \int 𝑢 (𝑖 𝑇 ˆ 𝑁 𝑢 𝑁 𝑢 - ̇ 𝑁 𝑢) \int 𝑣 (𝛿 𝑣 - 𝑢 + 𝜈 𝑣 - 𝑢 𝜃 𝑣 - 𝑢 𝑁 𝑣) 𝑖 ˆ 𝑁 𝑣 = 𝑆 diss [𝑁, ˆ 𝑁] .

(B19)

The action is invariant under the time-reversal transformation using $𝑃 eq$ as the initial and final condition, therefore $𝑃 eq$ is the correct equilibrium probability distribution.

APPENDIX C: RESPONSE FUNCTIONS

In the following, we restrict ourselves to the $𝜌 = 1$ case for simplicity, unless specified, and we show how to obtain the self-consistent equation leading to $𝑅 int 0$ .

At equilibrium we can rewrite the integrated disorder-dependent responses to a perturbation of the carrying capacity and of the immigration rate in terms of connected correlation functions of $𝑁$ ,

𝑟 int 𝑑 (𝑧) = \int \infty 0 𝑑 𝜏 ⟨ 𝛿 𝑁 𝑢 (𝜏) 𝛿 𝜁 𝑢 (0) ⟩ = 𝜕 ⟨ 𝑁 𝑢 ⟩ 𝜕 𝜁 𝑢 = 𝛽 (⟨ 𝑁 2 ⟩ - ⟨ 𝑁 ⟩ 2),

(C1)

𝜒 (𝑧) = \int \infty 0 𝑑 𝜏 ⟨ 𝛿 𝑁 𝑢 (𝜏) 𝛿 𝜆 𝑢 (0) ⟩ = 𝜕 ⟨ 𝑁 𝑢 ⟩ 𝜕 𝜆 𝑢 = 𝛽 (⟨ 𝑁 l o g 𝑁 ⟩ - ⟨ 𝑁 ⟩ ⟨ l o g 𝑁 ⟩) .

(C2)

When the time dependence is not present, we are considering a time-independent perturbation.

Adding a perturbation in site $𝑣$ leads to a variation of the abundances in all other sites, because of the coupling by diffusion and the memory term. These variations are of order $1 / 𝐿$ , but since there are $𝐿$ of them they give a significant contribution. When studying $𝜕 ⟨ 𝑁 𝑢 ⟩ 𝜕 𝜁 𝑣$ we need to take into account four contributions: there is an $𝑂 (1)$ variation of $𝑁 𝑣$ that leads to an $𝑂 (1 / 𝐿)$ perturbation of the immigration rate perceived by $𝑁 𝑢$ and an $𝑂 (1 / 𝐿)$ change in its off-diagonal memory term; there are $𝐿 - 2$ variations of $𝑂 (1 / 𝐿)$ of the $𝑁 𝑤$ , with $𝑤 \neq 𝑢, 𝑣$ , each leading to an $𝑂 (1 / 𝐿 2)$ change in both immigration and memory term. Carefully taking into account all these contributions, we can write $𝑟 int 0 (𝑧)$ in terms of $𝑟 int 𝑑 (𝑧), 𝜒 (𝑧)$ , and $𝑟 int 0 (𝑧)$ itself:

𝑟 int 0 (𝑧) = 𝐿 \int \infty 0 𝑑 𝜏 ⟨ 𝛿 𝑁 𝑢 (𝜏) 𝛿 𝜁 𝑣 (0) ⟩ = 𝐿 𝜕 ⟨ 𝑁 𝑢 ⟩ 𝜕 𝜁 𝑣 = 𝐿 ⟨ (𝜕 𝑁 𝑣 𝜕 𝜁 𝑣 (𝐷 𝐿 𝜕 𝑁 𝑢 𝜕 𝜆 𝑢 + 𝜎 2 𝑅 𝑖 𝑛 𝑡 𝑢 𝑣 𝜕 𝑁 𝑢 𝜕 𝜁 𝑢) + \sum 𝑤 \neq 𝑢, 𝑣 𝜕 𝑁 𝑤 𝜕 𝜁 𝑣 (𝐷 𝐿 𝜕 𝑁 𝑢 𝜕 𝜆 𝑢 + 𝜎 2 𝑅 𝑖 𝑛 𝑡 𝑢 𝑤 𝜕 𝑁 𝑢 𝜕 𝜁 𝑢)) ⟩ = (𝐷 𝜒 (𝑧) + 𝜎 2 𝑟 int 𝑑 (𝑧) 𝑅 int 0) (𝑟 int 𝑑 (𝑧) + 𝑟 int 0 (𝑧)) .

(C3)

In the third line we used the fact that the correlations between different patches are subleading to take separately the thermal averages. Solving for $𝑟 int 0 (𝑧)$ , we obtain

𝑟 int 0 (𝑧) = [𝐷 𝜒 (𝑧) + 𝜎 2 𝑟 int 𝑑 (𝑧) 𝑅 int 0] 𝑟 int 𝑑 (𝑧) 1 - [𝐷 𝜒 (𝑧) + 𝜎 2 𝑟 int 𝑑 (𝑧) 𝑅 int 0] .

(C4)

We can then average over $𝑧$ to obtain $𝑅 int 0$ :

𝑅 int 0 = ‾ [𝐷 𝜒 (𝑧) + 𝜎 2 𝑟 int 𝑑 (𝑧) 𝑅 int 0] 𝑟 int 𝑑 (𝑧) 1 - [𝐷 𝜒 (𝑧) + 𝜎 2 𝑟 int 𝑑 (𝑧) 𝑅 int 0] .

(C5)

APPENDIX D: ASYMMETRIC INTERACTIONS

The MSRDJ action with nonsymmetrical interactions is given by

𝑆 [𝑁, ˆ 𝑁] = \int 𝑢 𝑖 (ˆ 𝑁 𝑢 [𝑁 𝑢 (𝑘 - 𝐷 (1 - 𝜎 2 𝛾 𝑁 * 𝑅 int 0) - 𝜇 ℎ + 𝜎 \sqrt 𝐶 \infty 𝑑 𝑧 - 𝑁 𝑢) - 𝑇 + 𝐷 𝑁 *] + 𝑁 - 12)

(D1)

+ \int 𝑢 𝑖 ˆ 𝑁 𝑢 (𝑖 𝑇 ˆ 𝑁 𝑢 𝑁 𝑢 - ̇ 𝑁 𝑢) + 𝜎 2 2 \int 𝑢 𝑖 ˆ 𝑁 𝑢 𝑁 𝑢 \int 𝑣 𝐶 𝑐 (𝑢 - 𝑣) 𝑖 ˆ 𝑁 𝑣 𝑁 𝑣 + 𝛾 𝜎 2 \int 𝑢 𝑖 ˆ 𝑁 𝑢 𝑁 𝑢 \int 𝑣 𝑅 (𝑢 - 𝑣) 𝑁 𝑣 + [l o g 𝑃 (𝑁 (0))],

(D2)

where we have defined $𝐶 𝑐 (𝑢 - 𝑣) = 𝐶 𝑑 (𝑢 - 𝑣) - 𝐶 \infty 𝑑$ . If the introduction of a small asymmetry in the interactions ( $𝜀 = 1 - 𝛾 ≪ 1$ ) is a nonsingular perturbation, all the self-consistently determined quantities in the action ( $ℎ, 𝐶 𝑑, 𝑅 int 0$ , and $𝑅 𝑑$ ) will be close to their counterparts for $𝛾 = 1$ . At first order in $𝜀$ we can neglect their change; therefore, $𝑅 𝑑$ and $𝐶 𝑑$ will still respect a Fluctuation-Dissipation Relation. We can separate the action in a part that would respect FDT and a part that breaks it explicitly:

𝛿 𝑆 = 𝜀 𝜎 2 𝑇 \int 𝑢 > 𝑣 𝐶 𝑐 𝑢 - 𝑣 𝑖 ˆ 𝑁 𝑢 𝑁 𝑢 ̇ 𝑁 𝑣 .

(D3)

An average $⟨ 𝑓 (𝑁 𝑡) ⟩$ can be expanded as

⟨ 𝑓 (𝑁 𝑡) ⟩ = ⟨ 𝑓 (𝑁 𝑡) ⟩ 0 + ⟨ 𝑓 (𝑁 𝑡) 𝛿 𝑆 ⟩ 0 + 𝑂 (𝜀 2),

(D4)

where $⟨ \cdot ⟩ 0$ indicates the average with respect to the action neglecting $𝛿 𝑆$ .

We want to estimate the scaling of

⟨ 𝑓 (𝑁 𝑡) 𝛿 𝑆 ⟩ 0 = 𝜀 𝜎 2 𝑇 \int 𝑢 > 𝑣 𝐶 𝑐 𝑢 - 𝑣 𝑖 ⟨ 𝑓 (𝑁 𝑡) ˆ 𝑁 𝑢 𝑁 𝑢 ̇ 𝑁 𝑣 ⟩ 0 = 𝜀 𝜎 2 𝑇 \int 𝑢 > 𝑣 𝐶 𝑐 𝑢 - 𝑣 𝑖 𝜕 𝜕 𝑣 𝛿 𝛿 𝜁 𝑢 ⟨ 𝑓 (𝑁 𝑡) 𝑁 𝑣 ⟩ 0

(D5)

to show that it is not singular approaching a phase transition. In the simple equilibrium phase, the connected correlation function decays exponentially, with a typical relaxation time $𝜏$ that could diverge at the phase transitions:

𝐶 𝑐 (𝑢 - 𝑣) \sim (‾ ⟨ 𝑁 2 ⟩ - ⟨ 𝑁 ⟩ 2) 𝑒 - | 𝑢 - 𝑣 | / 𝜏 .

(D6)

The correlation function $⟨ 𝑓 (𝑁 𝑡) 𝑁 𝑣 ⟩ 0$ will contain a $𝑣$ independent part (that we can neglect since we will be taking the derivative in $𝑣$ ) and a connected component of order 1 that decays with the same relaxation time $𝜏$ . Perturbing the system with a field $𝜁 𝑢$ , this observable will respond as

𝛿 𝛿 𝜁 𝑢 ⟨ 𝑓 (𝑁 𝑡) 𝑁 𝑣 ⟩ 0 \propto 1 𝑇 𝜏 𝑒 - (𝑡 - 𝑣) / 𝜏 .

(D7)

Inserting these scalings in Eq. , we obtain

⟨ 𝑓 (𝑁 𝑡) 𝛿 𝑆 ⟩ 0 \propto 𝜀 𝜎 2 𝑇 \int 𝑢 > 𝑣 𝑒 - (𝑢 - 𝑣) / 𝜏 𝜕 𝜕 𝑣 (1 𝑇 𝜏 𝑒 - (𝑡 - 𝑣) / 𝜏) = 𝜀 𝜎 2 𝑇 2 𝜏 2 \int 𝑡 - \infty 𝑑 𝑣 𝑒 - (𝑡 - 𝑣) / 𝜏 \int 𝑡 𝑣 𝑑 𝑢 𝑒 - (𝑢 - 𝑣) / 𝜏

(D8)

= 𝜀 𝜎 2 𝑇 2 𝜏 \int 𝑡 - \infty 𝑑 𝑣 𝑒 - (𝑡 - 𝑣) / 𝜏 (1 - 𝑒 - (𝑡 - 𝑣) / 𝜏) = 𝜀 𝜎 2 𝑇 2 𝜏 (𝜏 - 𝜏 2) = 𝜀 𝜎 2 2 𝑇 2 .

(D9)

Considering a small asymmetry in the interactions, observables are shifted by a correction of order $𝜀$ , where the prefactor is of order 1 and has no divergence at the phase transitions. We thus expected the phase diagram to remain qualitatively unchanged.

APPENDIX E: EXTINCTION THRESHOLD AND DIVERSITY (FOR $𝜌 = 1$ )

The self-consistency condition for $𝑁 *$ reads

𝑁 * (𝑧) = ⟨ 𝑁 ⟩ 𝐻 eff (𝑁; ℎ, 𝐶 0 𝑑, 𝐶 \infty 𝑑, 𝑅 int 0, 𝑧, 𝑁 *) = \int \infty 0 𝑑 𝑁 𝑁 𝑒 - 𝛽 𝐻 eff (𝑁; ℎ, 𝐶 0 𝑑, 𝐶 \infty 𝑑, 𝑅 int 0, 𝑧, 𝑁 *) \int \infty 0 𝑑 𝑁 𝑒 - 𝛽 𝐻 eff (𝑁; ℎ, 𝐶 0 𝑑, 𝐶 \infty 𝑑, 𝑅 int 0, 𝑧, 𝑁 *) .

(E1)

$𝑁 * = 0$ is always a solution of this equation; we want to find the value of $𝑧$ at which it becomes unstable.

We can separate the effective Hamiltonian into a quadratic and a logarithmic part:

𝐻 eff (𝑁; ℎ, 𝐶 0 𝑑, 𝐶 \infty 𝑑, 𝑅 int 0, 𝑧, 𝑁 *) = 𝐻 𝑞 (𝑁; ℎ, 𝐶 0 𝑑, 𝐶 \infty 𝑑, 𝑅 int 0, 𝑧, 𝑁 *) + (𝑇 - 𝐷 𝑁 *) l n 𝑁 .

(E2)

For $𝑁 * \to 0$ , the logarithmic part gives rise to a nonintegrable divergence in 0 in the denominator. To improve the numerical stability of our solution at small $𝑁 *$ , we performed an integration by parts of the denominator:

\int \infty 0 𝑑 𝑁 𝑒 - 𝛽 𝐻 eff = \int \infty 0 𝑑 𝑁 𝑒 - 𝛽 𝐻 𝑞 𝑁 - 1 + 𝛽 𝐷 𝑁 * = 1 𝐷 𝑁 * \int \infty 0 𝑑 𝑁 𝑒 - 𝛽 𝐻 𝑞 𝑁 𝛽 𝐷 𝑁 * 𝑑 𝐻 𝑞 𝑑 𝑁 .

(E3)

The integral is now finite for $𝑁 * \to 0$ and we can expand Eq. in powers of $𝑁 *$ :

𝑁 * (𝑧) = 𝐷 𝑁 * (𝑧) \int \infty 0 𝑑 𝑁 𝑒 - 𝛽 𝐻 𝑞 𝑁 𝛽 𝐷 𝑁 * (𝑧) \int \infty 0 𝑑 𝑁 𝑒 - 𝛽 𝐻 𝑞 𝑑 𝐻 𝑞 𝑑 𝑁 𝑁 𝛽 𝐷 𝑁 * (𝑧) = 𝑁 * (𝑧) 𝐷 \int \infty 0 𝑑 𝑁 𝑒 - 𝛽 𝐻 𝑞 \int \infty 0 𝑑 𝑁 𝑒 - 𝛽 𝐻 𝑞 𝑑 𝐻 𝑞 𝑑 𝑁 + 𝑂 ([𝑁 * (𝑧)] 2) .

(E4)

The term of order $𝑁 * (𝑧) 2$ is always negative, therefore the number of solutions depends on the coefficient of the $𝑁 * (𝑧)$ term: if $𝑐 1 (𝑧) < 1$ the only solution is $𝑁 * (𝑧) = 0$ ; if $𝑐 1 (𝑧) > 1$ the $𝑁 * (𝑧) = 0$ solution is unstable and there is a positive stable one. We define the effective growth rate $𝑟 eff = 1 - 𝜎 2 𝛽 (𝐶 0 𝑑 - 𝐶 \infty 𝑑)$ and the effective growth factor $𝑟 eff 𝑔 eff (𝑧) = 𝑘 - 𝜇 ℎ - 𝑧 \sqrt 𝐶 \infty 𝑑 𝜎 + 𝐷 𝜎 2 𝑁 * 𝑅 int 0$ . The extinction threshold $𝑧 *$ (Fig. ) is given by

1 = 𝑐 1 (𝑧 *) = 𝐷 \sqrt 𝛽 𝜋 2 𝑟 eff e x p (𝛽 2 (𝑔 eff (𝑧 *) 𝑟 eff - 𝐷) 2 𝑟 eff) \times [1 + e r f (\sqrt 𝛽 2 (𝑔 eff (𝑧 *) 𝑟 eff - 𝐷 \sqrt 𝑟 eff))] .

(E5)

FIG. 11.

Self-consistent solution for $𝑁 * (𝑧)$ (blue), coefficient of the first-order expansion $𝑐 1 (𝑧)$ (orange), and Gaussian probability distribution $𝑃 (𝑧)$ (green). The highlighted region corresponds to the nonextinct species; its area is the diversity of the ecosystem. $𝑇 = 0.4, 𝐷 = 0.15, 𝜇 = 1, 𝜎 = 0.5$ .

As noted in Ref. , this is the same condition that would determine the criticality of the Directed Percolation process with corresponding growth rate and growth factor:

̇ 𝑁 𝑢 = 𝑟 eff (𝑁 2 𝑢 / 2 - 𝑔 eff 𝑁 𝑢) + 𝐷 (1 𝐿 \sum 𝑣 𝑁 𝑣 - 𝑁 𝑢) + 𝜂 \sqrt 𝑁 𝑢

(E6)

⟨ 𝜂 (𝑡) 𝜂 (𝑡') ⟩ = 2 𝑇 𝛿 (𝑡 - 𝑡') .

(E7)

The diversity, given by the fraction of nonextinct species, can be obtained as

𝜙 = \int \infty 𝑧 * 𝑃 (𝑧) 𝑑 𝑧 = 12 e r f c (𝑧 * \sqrt 2) .

(E8)

APPENDIX F: CONTINUOUS TRANSITION POINT

At the continuous transition, all moments of $𝑁$ tend to zero, and we can expand the extinction condition in powers of these moments. Keeping only the zeroth order, we obtain an equation on the critical value of the diffusion constant:

𝐷 0 \sqrt 𝛽 𝜋 2 𝑒 𝛽 2 (𝑘 - 𝐷 0) 2 [1 + e r f (\sqrt 𝛽 2 (𝑘 - 𝐷 0))] = 1 .

(F1)

This condition has no dependence on the distribution of the interactions; indeed, it is the same one that would be obtained with zero or constant interactions .

For $𝑇 \to 0$ we can expand Eq. and show that $𝐷 0$ vanishes exponentially:

𝐷 0 (𝑇) \sim 1 \sqrt 2 𝜋 𝛽 𝑒 - 𝑘 2 2 𝑇 .

(F2)

The reason for this behavior is that at low demographic noise, the abundances of a species with carrying capacity $𝑘$ undergo a fluctuation toward very low values very rarely. In fact, one needs to wait for a rare fluctuation of the demographic noise that makes the species go against the force due to the logistic growth. This phenomenon is similar to the one encountered in the Kramers' problem for barrier crossing. Using the same line of arguments employed there, one finds that the timescale for this rare event is $𝑒 𝑘 2 2 𝑇$ (the “energy barrier” equals $𝑘 2 / 2$ ). The equation above can therefore be interpreted as a balance between two inverse timescales: the one needed for diffusion to operate and the one over which extinctions take place.

By a careful (and cumbersome) expansion of the self-consistent equations, we can show that approaching the continuous transition $ℎ \propto 𝑞 𝑑 \propto 𝐷 - 𝐷 0, 𝑞 0 \propto (𝐷 - 𝐷 0) 2$ , and $𝑧 *$ (and therefore $𝜙$ ) has a finite limit.

APPENDIX G: ABUNDANCE DISTRIBUTION

As noted before, two types of stochasticity contribute to the distribution of abundances. Each species is subjected to demographic and environmental noise, making their abundance a time-dependent random variable. For each species, the abundance is distributed according to the Boltzmann distribution with Hamiltonian $𝐻 eff$ ; we will call this $𝑃 (𝑁 | 𝑧)$ . On top of this, because of disorder, different species experience different average interactions with the rest of the ecosystem [different values of $𝑧$ , i.e., distributed according to $𝑃 (𝑧)$ Gaussian], leading to species-dependent factors in $𝐻 eff$ . If we want to study the distribution of the abundances of all species at a given time in one site $[𝑃 (𝑁)]$ , we need to take into account both effects. We can compute $𝑃 (𝑁)$ marginalizing over $𝑧$ :

𝑃 (𝑁) = \int 𝑑 𝑧 𝑃 (𝑁 | 𝑧) 𝑃 (𝑧),

(G1)

with

𝑃 (𝑁 | 𝑧) = 𝑒 - 𝛽 𝐻 eff (𝑁; 𝑧, 𝑁 * (𝑧)) \int \infty 0 𝑑 𝑁 𝑒 - 𝛽 𝐻 eff (𝑁; 𝑧, 𝑁 * (𝑧))

(G2)

𝑃 (𝑧) = 𝑒 - 𝑧 2 / 2 \sqrt 2 𝜋 .

(G3)

We could also be interested in the distribution across species of the abundance averaged over patches or time, given by $𝑃 (𝑁 *) = 𝑃 (𝑧) (𝑑 𝑁 * (𝑧) 𝑑 𝑧) - 1$ for $𝑁 * > 0$ . There is also a finite probability $1 - 𝜙$ that $𝑁 * = 0$ , where $𝜙$ is the diversity.

Examples of these abundance probability distributions are shown in Fig. .

FIG. 12.

(a) Probability distribution of the abundance $𝑁$ ; in orange deep in the survival phase $[𝑇 = 0.8, 𝐷 / 𝐷 0 (𝑇) = 1.5]$ , in blue right before the discontinuous transition $[𝑇 = 0.4, 𝐷 / 𝐷 0 (𝑇) = 0.84]$ . Note that the shown distributions do not integrate to 1 because a finite fraction of the species are extinct, leading to a $𝛿$ function in zero with weight $1 - 𝜙$ . (b) Probability distribution of the abundance for a given species, i.e., at fixed $𝑧$ ; in blue for a species close to extinction ( $𝑧 = - 0.45, 𝑁 * = 0.148$ ), in orange for a species far from extinction ( $𝑧 = - 2.45, 𝑁 * = 1.733$ ); $𝑇 = 0.4, 𝐷 / 𝐷 0 (𝑇) = 0.84$ . (c) Probability distribution of the space (or time) averaged abundance, because of extinct species we again have a $𝛿$ function in zero with weight $1 - 𝜙$ . $𝑇 = 0.4, 𝐷 / 𝐷 0 (𝑇) = 0.84, 𝜇 = 1, 𝜎 = 0.5, 𝜌 = 1$ .

APPENDIX H: DIVERGENCE OF RESPONSE TO A VARIATION OF THE CARRYING CAPACITY

The divergence of response functions when approaching a tipping point is a generic feature of saddle node bifurcations . Let us consider a generic dynamical system, described by

𝑑 ⃗ 𝑥 𝑑 𝑡 = 𝐹 (⃗ 𝑥, 𝑘) .

(H1)

$⃗ 𝑥$ contains all the degrees of freedom of the system, whereas $𝑘$ is a control parameter. The zeros of $𝐹$ yield the stationary states $𝑥 *$ :

𝐹 (𝑥 *, 𝑘) = 0 .

(H2)

The stationary point is stable if the Jacobian of $𝐹$ has only negative eigenvalues, ensuring that $𝑥$ returns to $𝑥 *$ upon small perturbations. In a saddle node bifurcation, a stable and an unstable stationary point collide and annihilate each other. Since the Jacobian at the stable stationary point has only negative eigenvalues whereas the one at the unstable stationary point has at least a positive eigenvalue, one of the eigenvalues has to cross zero at the bifurcation. The existence of a zero mode leads to a diverging response to perturbation.

We show below how this mechanism is at play in our case when approaching the stability limit of the self-sustained phase. We study the response of the system to a perturbation in the environmental conditions in the case of independent interaction coefficients ( $𝜌 = 0$ ). We will consider for concreteness a perturbation to the carrying capacity $𝑘$ , but we expect the same qualitative behavior for perturbations to the diffusion constant, the moments of the interactions, or the strength of the demographic fluctuations.

The response of the order parameters to a variation of $𝑘$ involves some connected moments of $𝑁$ and the derivative of $𝐻$ in $𝑘$ :

𝑑 ℎ 𝑑 𝑘 = \int 𝒟 𝑧 (\int 𝑑 𝑁 𝑁 𝑒 - 𝛽 𝐻 (- 𝛽) 𝑑 𝐻 / 𝑑 𝑘 \int 𝑑 𝑁 𝑒 - 𝛽 𝐻 - \int 𝑑 𝑁 𝑁 𝑒 - 𝛽 𝐻 \int 𝑑 𝑁 𝑒 - 𝛽 𝐻 (- 𝛽) 𝑑 𝐻 / 𝑑 𝑘 (\int 𝑑 𝑁 𝑒 - 𝛽 𝐻) 2) = - 𝛽 (‾ ⟨ 𝑁 𝑑 𝐻 𝑑 𝑘 ⟩ - ‾ ⟨ 𝑁 ⟩ ⟨ 𝑑 𝐻 𝑑 𝑘 ⟩),

(H3)

𝑑 𝐶 0 𝑑 𝑑 𝑘 = - 𝛽 (‾ ⟨ 𝑁 2 𝑑 𝐻 𝑑 𝑘 ⟩ - ‾ ⟨ 𝑁 2 ⟩ ⟨ 𝑑 𝐻 𝑑 𝑘 ⟩),

(H4)

𝑑 𝐶 \infty 𝑑 𝑑 𝑘 = - 2 𝛽 (‾ ⟨ 𝑁 ⟩ ⟨ 𝑁 𝑑 𝐻 𝑑 𝑘 ⟩ - ‾ ⟨ 𝑁 ⟩ 2 ⟨ 𝑑 𝐻 𝑑 𝑘 ⟩) .

(H5)

Thanks to the fact that $𝜌 = 0, 𝐻$ has the simplified form

𝐻 = [1 - 𝜎 2 𝛽 (𝐶 0 𝑑 - 𝐶 \infty 𝑑)] 𝑁 2 / 2 + (𝜇 ℎ - 𝑘 + 𝐷 - 𝑧 \sqrt 𝐶 \infty 𝑑 𝜎) 𝑁 + (𝑇 - 𝐷 ℎ) l n 𝑁 .

(H6)

$𝑑 𝐻 / 𝑑 𝑘$ depends on the derivative of the order parameters in $𝑘$ :

𝑑 𝐻 𝑑 𝑘 = - 𝑁 + 𝜕 𝐻 𝜕 ℎ 𝑑 ℎ 𝑑 𝑘 + 𝜕 𝐻 𝜕 𝐶 0 𝑑 𝑑 𝐶 0 𝑑 𝑑 𝑘 + 𝜕 𝐻 𝜕 𝐶 \infty 𝑑 𝑑 𝐶 \infty 𝑑 𝑑 𝑘,

(H7)

𝜕 𝐻 𝜕 ℎ = 𝜇 𝑁 - 𝐷 l n 𝑁,

(H8)

𝜕 𝐻 𝜕 𝐶 0 𝑑 = - 𝜎 2 𝛽 2 𝑁 2,

(H9)

𝜕 𝐻 𝜕 𝐶 \infty 𝑑 = 𝜎 2 𝛽 2 𝑁 2 - 1 2 \sqrt 𝐶 \infty 𝑑 𝑧 𝜎 𝑁 .

(H10)

Substituting $𝑑 𝐻 / 𝑑 𝑘$ in Eqs. , we obtain

𝑑 ℎ 𝑑 𝑘 = - 𝛽 {- (‾ ⟨ 𝑁 2 ⟩ - ‾ ⟨ 𝑁 ⟩ 2) + [𝜇 (‾ ⟨ 𝑁 2 ⟩ - ‾ ⟨ 𝑁 ⟩ 2) - 𝐷 (‾ ⟨ 𝑁 l o g 𝑁 ⟩ - ‾ ⟨ 𝑁 ⟩ ⟨ l o g 𝑁 ⟩)] 𝑑 ℎ 𝑑 𝑘 - 𝜎 2 𝛽 2 (‾ ⟨ 𝑁 3 ⟩ - ‾ ⟨ 𝑁 ⟩ ⟨ 𝑁 2 ⟩) 𝑑 𝐶 0 𝑑 𝑑 𝑘 + [𝜎 2 𝛽 2 (‾ ⟨ 𝑁 3 ⟩ - ‾ ⟨ 𝑁 ⟩ ⟨ 𝑁 2 ⟩) - 1 2 \sqrt 𝐶 \infty 𝑑 𝜎 (‾ ⟨ 𝑁 2 ⟩ 𝑧 - ‾ ⟨ 𝑁 ⟩ 2 𝑧)] 𝑑 𝐶 \infty 𝑑 𝑑 𝑘},

(H11)

𝑑 𝐶 0 𝑑 𝑑 𝑘 = - 𝛽 {- (‾ ⟨ 𝑁 3 ⟩ - ‾ ⟨ 𝑁 2 ⟩ ⟨ 𝑁 ⟩) + [𝜇 (‾ ⟨ 𝑁 3 ⟩ - ‾ ⟨ 𝑁 2 ⟩ ⟨ 𝑁 ⟩) - 𝐷 (‾ ⟨ 𝑁 2 l o g 𝑁 ⟩ - ‾ ⟨ 𝑁 2 ⟩ ⟨ l o g 𝑁 ⟩)] 𝑑 ℎ 𝑑 𝑘 - 𝜎 2 𝛽 2 (‾ ⟨ 𝑁 4 ⟩ - ‾ ⟨ 𝑁 2 ⟩ 2) 𝑑 𝐶 0 𝑑 𝑑 𝑘 + [𝜎 2 𝛽 2 (‾ ⟨ 𝑁 4 ⟩ - ‾ ⟨ 𝑁 2 ⟩ 2) - 1 2 \sqrt 𝐶 \infty 𝑑 𝜎 (‾ ⟨ 𝑁 3 ⟩ 𝑧 - ‾ ⟨ 𝑁 2 ⟩ ⟨ 𝑁 ⟩ 𝑧)] 𝑑 𝐶 \infty 𝑑 𝑑 𝑘},

(H12)

𝑑 𝐶 \infty 𝑑 𝑑 𝑘 = - 2 𝛽 {- (‾ ⟨ 𝑁 ⟩ ⟨ 𝑁 2 ⟩ - ‾ ⟨ 𝑁 ⟩ 3) + [𝜇 (‾ ⟨ 𝑁 ⟩ ⟨ 𝑁 2 ⟩ - ‾ ⟨ 𝑁 ⟩ 3) - 𝐷 (‾ ⟨ 𝑁 ⟩ ⟨ 𝑁 l o g 𝑁 ⟩ - ‾ ⟨ 𝑁 ⟩ 2 ⟨ l o g 𝑁 ⟩)] 𝑑 ℎ 𝑑 𝑘 - 𝜎 2 𝛽 2 (‾ ⟨ 𝑁 ⟩ ⟨ 𝑁 3 ⟩ - ‾ ⟨ 𝑁 ⟩ 2 ⟨ 𝑁 2 ⟩) 𝑑 𝐶 0 𝑑 𝑑 𝑘 + [𝜎 2 𝛽 2 (‾ ⟨ 𝑁 ⟩ ⟨ 𝑁 3 ⟩ - ‾ ⟨ 𝑁 ⟩ 2 ⟨ 𝑁 2 ⟩) - 1 2 \sqrt 𝐶 \infty 𝑑 𝜎 (‾ ⟨ 𝑁 ⟩ ⟨ 𝑁 2 ⟩ 𝑧 - ‾ ⟨ 𝑁 ⟩ 3 𝑧)] 𝑑 𝐶 \infty 𝑑 𝑑 𝑘} .

(H13)

We collect the three order parameters in a vector $⃗ 𝑝 = (ℎ, 𝐶 0 𝑑, 𝑞 0) 𝑇$ . Then $𝑑 ⃗ 𝑝 𝑑 𝑘$ satisfies

𝑑 ⃗ 𝑝 𝑑 𝑘 = ̂ 𝐽 𝑑 ⃗ 𝑝 𝑑 𝑘 + ⃗ 𝑠,

(H14)

where $̂ 𝐽$ is a $3 \times 3$ matrix and $𝑠$ is a vector; their elements are the coefficients of Eqs. . The solution is given by

𝑑 ⃗ 𝑝 𝑑 𝑘 = - (̂ 𝐽 - ̂ 1) - 1 ⃗ 𝑠 .

(H15)

The response to a variation of $𝑘$ diverges if $̂ 𝐽 - ̂ 1$ has a zero eigenvalue, which is found to happen when approaching the discontinuous transition.

We expect the same qualitative behavior of the response to perturbations for generic values of $𝜌$ , but for $𝜌 \neq 0$ we need to take into account also the variations of the function $𝑁 * (𝑧)$ , which leads to the study of an infinite-dimensional matrix.

APPENDIX I: REDUCED INTERACTION MATRIX

In the case of fixed interaction matrices, a finite fraction of the species goes extinct; the interaction matrix restricted to the surviving species has a smaller mean than the starting one. The statistics of the reduced interaction matrix can be computed at 0 temperature :

𝜇' = 𝜙 𝜇 - 2 𝜎 2 ℎ 𝜙 𝑑 𝜙 𝑑 𝜁 = 𝜙 𝜇 - 2 𝜎 ℎ 𝜙 \sqrt 2 𝜋 𝐶 \infty 𝑑 𝑒 - (𝑧 *) 2 / 2 .

(I1)

Since we are not at 0 temperature, in our case this formula is only an approximation, but it provides a useful estimate of the variation of the mean interaction. We find that the interaction mean decreases (more mutualistic) when decreasing the diffusion coefficient [Fig. ]; it is negative in the entire metastability region. In Fig. we show the distribution of the interaction coefficients considering all species or only surviving ones in numerical simulations. The distribution of the interaction coefficients is slightly shifted to more negative values, and indeed $𝜇' = 𝑆 1 𝑆 2 \sum 𝑖 𝑗 𝛼 𝑖 𝑗$ changes from 0.96 to $- 0.28$ .

FIG. 13.

Left: Analytical estimate of the mean of the reduced interaction matrix using Eq. for $𝑇 = 0.4$ . Right: Numerical results for the distribution of the interaction coefficients for all and surviving species. $𝑇 = 0.18, 𝐷 / 𝐷 0 = 0.8$ . In both cases, in the initial species pool $𝜇 = 1, 𝜎 = 0.5$ ; if all species go extinct, we say $𝜇' = 0$ .

To compute the average interaction term, we can again use the cavity method and imagine to add a species (with index 0) to the community. Using Eq. ,

Int 0 = ⟨ \sum 𝑗 𝛼 0 𝑗 𝑁 𝑢 𝑗 ⟩ = 𝜇 ℎ + 𝜎 \sqrt 𝐶 \infty 𝑑 𝑧 - 𝛾 𝜎 2 (𝑅 int 𝑑 + 𝑅 int 0) ⟨ 𝑁 0 ⟩ .

(I2)

We can now average it over all species (all values of $𝑧$ , overline), or over only nonextinct ones ( $𝑧 < 𝑧 *$ , overline with $+$ superscript),

‾ 𝐼 = 𝜇 ℎ - 𝛾 𝜎 2 (𝑅 int 𝑑 + 𝑅 int 0) ℎ,

(I3)

‾ 𝐼 + = 𝜇 ℎ - 𝜎 \sqrt 𝑞 0 𝜙 𝑒 - 𝑧 * 2 / 2 \sqrt 2 𝜋 - 𝛾 𝜎 2 (𝑅 int 𝑑 + 𝑅 int 0) ℎ 𝜙 .

(I4)

Note that we will always find $‾ 𝐼 + < ‾ 𝐼$ ; $‾ 𝐼 +$ is negative in the entire metastability region [Fig. in the main text]. This is also confirmed by numerical simulations: the average interaction term is 0.13 considering all species, and $- 0.46$ considering only nonextinct ones [Fig. in the main text].

In the case of independent interaction matrices, all species survive, so that the interaction matrix is not modified.

APPENDIX J: NUMERICAL SCHEME

The numerical simulation of demographic noise poses some technical challenges. Naively sampling it as a Gaussian variable can result in negative species abundances, an unphysical result that makes the scheme numerically unstable. A clever solution was found in Ref. , and improved in . The idea is to separate the process in a deterministic part:

̇ 𝑁 𝑖, 𝑢 = 𝑁 𝑖, 𝑢 (1 - 𝑁 𝑖, 𝑢 - \sum 𝑗 𝛼 𝑢 𝑖 𝑗 𝑁 𝑗, 𝑢) + 𝐷 (1 𝐿 \sum 𝑣 𝑁 𝑖, 𝑣 - 𝑁 𝑖, 𝑢)

(J1)

and a stochastic one:

̇ 𝑁 𝑖, 𝑢 = \sqrt 𝑁 𝑖, 𝑢 𝜂 𝑖, 𝑢 .

(J2)

At each time step we numerically integrate the two in sequence. For the stochastic part, an exact solution of the associated Fokker-Planck equation is available for any initial condition, and it can be efficiently sampled using Gamma and Poisson variables:

˜ 𝑁 𝑖, 𝑢 (𝑡) = Gamma [Poisson (𝑁 𝑖, 𝑢 (𝑡) 𝑇 𝑑 𝑡)] 𝑇 𝑑 𝑡 .

(J3)

For the deterministic part, we rely on the Euler method,

𝑁 𝑖, 𝑢 (𝑡 + 𝑑 𝑡) = [˜ 𝑁 𝑖, 𝑢 (𝑡) (1 - ˜ 𝑁 𝑖, 𝑢 (𝑡) - \sum 𝑗 𝛼 𝑢 𝑖 𝑗 ˜ 𝑁 𝑗, 𝑢 (𝑡)) + 𝐷 (1 𝐿 \sum 𝑣 ˜ 𝑁 𝑖, 𝑣 (𝑡) - ˜ 𝑁 𝑖, 𝑢 (𝑡))] 𝑑 𝑡 .

(J4)

APPENDIX K: ADDITIONAL NUMERICAL RESULTS

Some of the challenges encountered in numerical simulations become clear examining the time evolution of the average abundances (Figs. and ). At high temperature (top) the average abundance fluctuates significantly even with a large number of species and patches ( $𝑆 = 200, 𝐿 = 400$ ); finite size effects on $𝐿$ lead to an excess of extinctions. At lower temperature (bottom), the dynamics strongly slows down, and at $𝑡 = 200$ some of the abundances (depending on the value of $𝐷$ ) have not yet reached their asymptotic value, leading to a smoothing of the discontinuous transition. Without heterogeneity in the interaction network we observe no dependence on initial conditions, even at $𝑇 = 0.18$ (Fig. ).

FIG. 14.

Time evolution of the average abundance for two different temperatures and two average values of the initial conditions. Note the different time ranges in the top and bottom figures: at high temperature the abundances have converged to their asymptotic values at $𝑡 max = 200$ , while at lower temperature it is necessary to wait much longer ( $𝑡 max = 500$ ). $𝑆 = 200, 𝐿 = 400, 𝜇 = 1, 𝜎 = 0.5$ .

FIG. 15.

Time evolution of the average abundance with partial correlation between patches ( $𝜌 = 0.9$ ) and nonsymmetric interactions ( $𝛾 = 0.9)$ at $𝑇 = 0.18$ and two average values of the initial conditions. At $𝑡 = 500$ the abundances have not yet reached their asymptotic value, leading to an apparent smoothing of the discontinuous transition. Nevertheless, this is ensured by the abrupt change of behavior of the evolution of the average abundance: for one value of the diffusion constant at long times, the abundance is decaying to 0, whereas for the next it shows a (slow) increase. We conclude that the asymptotic values would likewise show an abrupt change. $𝑆 = 200, 𝐿 = 400, 𝜇 = 1, 𝜎 = 0.5$ .

FIG. 16.

Time evolution of the average abundance without heterogeneity in the interaction network ( $𝜎 = 0$ ) at $𝑇 = 0.18$ and two average values of the initial conditions. For $𝐷 < 𝐷 0 (𝑇)$ the abundances converge to 0. $𝑆 = 200, 𝐿 = 400, 𝜇 = 1$ .

Interactions and Migration Rescuing Ecological Diversity

Abstract

Physics Subject Headings (PhySH)

Viewpoint

How Migration May Stabilize the Diversity of Ecosystems

Article Text

I. INTRODUCTION

II. THE MODEL

III. METHODS

A. DMFT and coupled Directed Percolation processes

B. Symmetric interactions, mean-field approximation, and mapping to a system in thermal equilibrium

C. Extensions

1. Spatial heterogeneity

2. Asymmetric interactions

IV. RESULTS

A. Characterization of the self-sustained phase

B. Transition to complete extinction: Emergence of a discontinuous transition at low dispersal

V. ASSESSING THE GENERALITY OF THE SCENARIO

A. Finite number of species and finite number of patches

B. Asymmetric interactions and partial correlation between patches

VI. PRECURSOR OF THE INSTABILITY TOWARD EXTINCTION

VII. CONCLUSIONS

ACKNOWLEDGMENTS

Appendices

APPENDIX A: DMFT DERIVATION

APPENDIX B: STATIONARY PROBABILITY DISTRIBUTION IN THE SYMMETRIC CASE

APPENDIX C: RESPONSE FUNCTIONS

APPENDIX D: ASYMMETRIC INTERACTIONS

APPENDIX E: EXTINCTION THRESHOLD AND DIVERSITY (FOR $𝜌 = 1$ )

APPENDIX F: CONTINUOUS TRANSITION POINT

APPENDIX G: ABUNDANCE DISTRIBUTION

APPENDIX H: DIVERGENCE OF RESPONSE TO A VARIATION OF THE CARRYING CAPACITY

APPENDIX I: REDUCED INTERACTION MATRIX

APPENDIX J: NUMERICAL SCHEME

APPENDIX K: ADDITIONAL NUMERICAL RESULTS

References (80)

Interactions and Migration Rescuing Ecological Diversity

Abstract

Physics Subject Headings (PhySH)

Viewpoint

How Migration May Stabilize the Diversity of Ecosystems

Article Text

I. INTRODUCTION

II. THE MODEL

III. METHODS

A. DMFT and coupled Directed Percolation processes

B. Symmetric interactions, mean-field approximation, and mapping to a system in thermal equilibrium

C. Extensions

1. Spatial heterogeneity

2. Asymmetric interactions

IV. RESULTS

A. Characterization of the self-sustained phase

B. Transition to complete extinction: Emergence of a discontinuous transition at low dispersal

V. ASSESSING THE GENERALITY OF THE SCENARIO

A. Finite number of species and finite number of patches

B. Asymmetric interactions and partial correlation between patches

VI. PRECURSOR OF THE INSTABILITY TOWARD EXTINCTION

VII. CONCLUSIONS

ACKNOWLEDGMENTS

Appendices

APPENDIX A: DMFT DERIVATION

APPENDIX B: STATIONARY PROBABILITY DISTRIBUTION IN THE SYMMETRIC CASE

APPENDIX C: RESPONSE FUNCTIONS

APPENDIX D: ASYMMETRIC INTERACTIONS

APPENDIX E: EXTINCTION THRESHOLD AND DIVERSITY (FOR 𝜌=1)

APPENDIX F: CONTINUOUS TRANSITION POINT

APPENDIX G: ABUNDANCE DISTRIBUTION

APPENDIX H: DIVERGENCE OF RESPONSE TO A VARIATION OF THE CARRYING CAPACITY

APPENDIX I: REDUCED INTERACTION MATRIX

APPENDIX J: NUMERICAL SCHEME

APPENDIX K: ADDITIONAL NUMERICAL RESULTS

References (80)

APPENDIX E: EXTINCTION THRESHOLD AND DIVERSITY (FOR $𝜌 = 1$ )