Interactions and Migration Rescuing Ecological Diversity

How diversity is maintained in natural ecosystems is a long-standing question in Theoretical Ecology. By studying a system that combines ecological dynamics, heterogeneous interactions, and spatial structure, we uncover a new mechanism for the survival of diversity-rich ecosystems in the presence of demographic ﬂuc-tuations. For a single species, one ﬁnds a continuous phase transition between an extinction and a survival state, that falls into the universality class of Directed Percolation. Here we show that the case of many species with heterogeneous interactions is different and richer. By merging theory and simulations, we demonstrate that with sufﬁciently strong demographic noise, the system exhibits behavior akin to the single-species case, undergoing a continuous transition. Conversely, at low demographic noise, we observe unique features indicative of the ecosystem’s complexity. The combined effects of the heterogeneity in the interaction network and migration enable the community to thrive, even in situations where demographic noise would lead to the extinction of isolated species. The emergence of mutualism induces the development of global bistability, accompanied by sudden tipping points. We present a way to predict the catastrophic shift from high diversity to extinction by probing responses to perturbations as an early warning signal.


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 [1].
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 [2][3][4][5][6][7][8].Many ecosystems are spatially extended: the ecological dynamics takes place at some local scale, but individuals can then explore different spatial locations through migration [9].This can lead to the appearance of complex ecological phenomena, such as traveling activity fronts, pattern formation, and persistent chaotic dynamics [10][11][12][13][14][15][16][17][18].Ecosystems are subject to noise, in particular environmental and demographic (due to stochasticity in births and deaths).Both noises induce fluctuations which are a key factor in determining abundances distributions, and their time-dependence [19][20][21][22][23][24][25][26][27][28][29].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 speciesrich 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 for example the case in semi-arid ecosystems: the total number of plants is such that global fluctuations are negligible, but at the local level stochasticity can play a fundamental role [26].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 [27].
Previous works, following the classical theory of Island Biogeography by MacArthur and Wilson [30], proposed as a rescuing mechanism the immigration from a static reservoir (or "mainland", when thinking of an islandmainland system) [2,4,27,31,32].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 re-invade, providing an "insurance" (or "storage") effect [33,34].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 stabilisation of high-diversity states by spatial structure is a very general phenomenon: it can arise in the presence of spatial heterogeneity of environmental conditions [9,[33][34][35][36] or when abundances in different spatial locations exhibit unsynchronized fluctuations [14][15][16]37].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 a extinct state.A transition separates the two regimes [11,13,38,39].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 [40].
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 [16].Because the transition is continuous with vanishing abundances, interactions, that are quadratic in the abundances, are subleading at the critical point.In 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 [41][42][43][44], it is important to find early warning signals of these kinds of transition 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 [45].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 [46].This disorder approach, which dates back to May's seminal paper [47], has recently inspired a growing body of work [3-6, 23, 48-51] and also received positive experimental confirmations [31,52].Previous works have explored within this framework the effect of heterogeneous interactions [3][4][5][6]47], demographic fluctuations [5,23] and spatial structure [12][13][14][15][16][17], 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, multi-stability, and chaos [3][4][5][6]31].We expect this complex behavior also in the case of spatially structured ecosystems [12].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 [53][54][55][56].Our results could therefore find applications beyond ecology, notably for the study of global bistability and crashes in economy.We consider a meta-community of S species living on a network of L discrete spatial locations, or patches.A graphical representation of the system is given in Figure 1 in the case of a fully connected network of 3 patches.Each species is characterized by its abundance in each patch, which is modeled by a continuous variable, N i,u , representing the total number of individuals divided by the typical size of the local population Ñtyp .The abun-dance of species i in patch u evolves according to the stochastic differential equation:

D D
which corresponds to Lotka-Volterra dynamics, with constant growth rate r and carrying capacity k that are set to 1 throughout.The notation ∂u indicates the set of patch neighbors of u (from and to which species in patch u can migrate).The growth of each species is influenced by the abundance of all the others through the interaction coefficients α u ij : if α u ij is positive species j inhibits the growth of species i in patch u and vice versa.Positive α u ij and α u ji correspond to two species competing for resources, whereas α u ij and α u ji both negative correspond to mutualistic behaviour.Predation leads to opposite signs.
To model the heterogeneity in the interactions of species-rich ecosystems, we follow [2,57] and consider the disordered LV model.As already discussed in the introduction, the disorder approach has attracted recently a lot of attention [3-6, 23, 48, 49] and also received positive experimental confirmations [31,52].In this framework, the interaction coefficients are random variables, with mean µ/S and variance σ 2 /S.They are independent in each patch except for α u ij and α u ji , 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...) 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 α u ij and α v ij have a correlation coefficient ρ [14,15].
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. [57]).Without spatial heterogeneity the transition point is not modified by the introduction of a spatial structure [12], and spatial heterogeneity decreases the effective complexity of the interaction network [35], 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 non symmetric interactions, in which a single community with fixed immigration exhibits chaotic dynamics, [3][4][5][6]31] was studied in [14,15] 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.( 1) individuals can migrate on the patches network through diffusion, with a constant diffusion coefficient D/c, where c 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 u to any of its c nearest neighbors v ∈ ∂u.
Each species is subject to a white demographic noise η u i , accounting for the stochasticity in birth and death events in a continuum setting [5,23].We follow Ito's convention, according to which fluctuations in birth and deaths at time t + dt depend on the abundance at the previous time step.The noise is uncorrelated and of constant amplitude across species and patches: The auto-correlation of the demographic noise defines the noise strength T which depends on the birth and death rates and on the typical size of the local population; T scales as T ∝ 1/ Ñtyp [5,23]: the larger the local populations, the more negligible are demographic fluctuations.
In the γ = 1 case T 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 [58,59].Therefore also in the continuous model, extinction is possible over finite times, and not only asymptotically as it would be the case for example with environmental noise.physics and statistical field theory [40].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 Figure 2 for a one-dimensional network.Directed percolation was originally introduced to model spreading phenomena, from forest fires to epidemics [40].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 [40,60].We show in figure 3) 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 [40,60].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.(1) in absence of interactions, one for each species.Each equation corresponds to an independent directed percolation process.
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 pres-ence 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.

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.In order 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, see the appendix) 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 [46].
The large S limit allows for an analytical treatment, as the dynamics of the S interacting degrees of freedom can be replaced by the effective dynamics for a single representative species, through Dynamical Mean Field Theory (DMFT) [61].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 [62], 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 selfconsistent 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 [63,64].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 reference [61] for LV models, and can be found in Appendix A 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 S → ∞ DMFT allows one to replace the interaction term − j α ij N j,u by a stochastic expression that has the same statistical properties: Since this allows us to decouple different species, we will for simplicity omit the species index i in the following.We now discuss the different contributions.Note that in the following empirical averages over species will be denoted as 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, h = E[N u ], that in the steady state does not depend on the patch u 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 auto-correlation of the single species abundances: The noise ξu (t) 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 [61], 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 [61].We decompose the noise by rewriting ξu where ) is the value of the correlation function within the same patch at infinite times, z is a static Gaussian variable with zero mean and unit variance, that now plays the role of quenched disorder, and ξ u (t) is a fluctuating noise whose covariance vanishes at long times.Again z and C ∞ d do not depend on the patch u thanks to translational invariance.
To distinguish the roles of fluctuating and static noises in the GLV equation, we introduce two kinds of averages: ⟨•⟩ 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 • instead stands for the average over the static field z 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 i influences species j, which in turn influences species i.These contributions sum up because of the correlation between α ij and its reciprocal α ji , leading to the factor γ.This feedback mechanism (the famous Onsager reaction in the spinglass literature) generates a memory term, containing the response function of the abundance on patch u to a perturbation in the carrying capacity in patch v: In the S → ∞ limit, there is convergence in law between the statistics of the infinite number of randomly coupled DP processes and the effective one [62,65], i.e. the dynamics of a species satisfying the GLV equation ( 1) is equivalent to the effective one of a single species living on the original spatial network: The DMFT closure consists then in replacing the empirical averages over species E[•] 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.
Eq. ( 6) 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 z is to change the control parameter of the DP process, determining whether this is sub-critical or supercritical.
Interestingly, whereas a system of few species interacting and diffusing on a network was established to boil down to a standard DP problem [11,13,38,39,66], 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 na-ture 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 [40], 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 ( 6) 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 D c v∈∂u N v by its thermal av-erage.This amount to D c v∈∂u N v → DN * , where N * = 1 c v∈∂u ⟨N v ⟩ and, using translation invariance, it simplifies to ⟨N u ⟩ (which is time-independent since we are considering steady states).This procedure corresponds to a mean-field approximation of the spatially dependent DMFT Eqs.(6).Such DMFT 2 approximation becomes exact for a fully connected network.In fact, in this case, taking the L → ∞ limit, the empirical average of the abundances over the patches concentrates around the thermal average N * = ⟨N u ⟩.From now on, we shall focus on this case.
By substituting D c v∈∂u N v with N * in Equation ( 6), one obtains an equation on N u only, with an additional parameter to be determined self-consistently.Note that N * is obtained by averaging only over thermal fluctuations, and not over disorder: therefore, it will have to be determined as a function of z.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 u, we now obtain (for large times, i.e. in the steady state): Since all patches are equivalent on a fully connected lattice, the R uv and C uv matrices (of functions) only have two independent elements: the diagonal ones, R d and C d , and the off-diagonal ones R 0 /L and C 0 (see the appendix for the justification of the scaling with L of R 0 /L and C 0 ).
In the case of symmetric interactions, γ = 1, one can show (see App. B) 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: The memory term and ξ therefore play the role of a friction term and the noise associated with a colored thermal bath at temperature T .The stochastic process maps then to a generalized Langevin equation whose stationary probability distribution is given by the Boltzmann distribution at temperature T and with the effective Hamiltonian: where C 0 d is the equal-time correlation function, namely the second moment of the abundances over disorder and noise, ⟨N 2 ⟩.The long time limit of the correlation function, C ∞ d , represents instead the second moment of the thermal-averaged abundances, ⟨N ⟩ 2 .R int 0 is the integrated off-diagonal response, which is the solution of the self-consistent equation (see Appendix C): χ(z) and r d (z) are the species-dependent response to a perturbation in the immigration rate or the carrying capacity, respectively: The self-consistent equations can be expressed as av-erages with respect to the Boltzmann distribution: and analogously for R int 0 .Dz = dz √ 2π e −z 2 /2 indicates the average over the Gaussian field.
These equations can be solved iteratively: starting from a suitable initial condition for and R int 0 , one updates their values according to equations ( 13)-( 16) until reaching a fixed point.Because very large values of z are exponentially suppressed by the Gaussian distribution, it is sufficient to determine N * (z) for z of O(1).
In conclusion, within the DMFT 2 approximation and for the symmetric case, the formidable self-consistent stochastic equations ( 6) can be analyzed by studying a set of static self-consistent equations on four parameters and one function N * (z).Solving these equations (see next section) allows us to obtain a general picture of the interplay between migration and demographic noise for spatially extended metacommunities.In order 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.

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, that we can decompose as where z is constant and w u independent across locations, and C ∞ d and C ∞ 0 are the infinite time correlation function of the abundance on the same patch and on different patches, that for ρ = 1 coincide.Averaging the abundance across patches to obtain the immigration rate requires an additional step, i.e. averaging also over w u .The solution of the self-consistent equations, albeit conceptually analogous to the ρ = 1 case, is for generic values of ρ much more numerically challenging, because of the need to integrate over two disorder fields, z and w u .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.

Asymmetric interactions
The mapping to an equilibrium distribution requires symmetry in the interactions: non-symmetric interactions correspond to non-conservative forces, which explicitly break time reversal and lead to non-equilibrium steady states.In order 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 [67][68][69][70] 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 D 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.), the metastability region has a vanishing width in this limit and the system is always in the survival phase.µ = 1, σ = 0.5.

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 non-zero 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, h = ⟨N ⟩, and the ecosystem diversity ϕ, defined as the fraction of species that are not globally extinct, i.e. that have non zero abundance in at least one patch.At stationarity, we can compute the ecosystem diversity as ϕ = θ(⟨N ⟩).
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 bottom panels of Fig. 4. On the contrary, dispersal is beneficial, as shown in the top panels of Fig. 4. 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 [32].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 zσ √ q 0 , which renormalizes the carrying capacity of a species.For ρ = 1, if z is larger than a critical value z * the corresponding species goes extinct (for z > z * the renormalized carrying capacity is negative).This is true also for smaller values of ρ (Appendix E).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 [9,14,35,36].
The limits D → ∞ and T → 0 can be mapped to the well-mixed case.For D → ∞ 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, that 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 T → 0 limit (see two bottom panels of Figure 4): 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 G), as it is the case in other models with random fully connected interactions [3][4][5]71].

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 that would be obtained in the absence of interactions, where the system directly maps to directed percolation, or in the case of constant interactions [16], see Figure 2 and Appendix F. 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 (Figure 5c).Interestingly, approaching the transition the diversity does not go to zero and instead tends to a finite value (Figure 5e).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.
At smaller demographic noise this picture changes drastically and interactions play a major role, as shown in the phase diagram in Figure 5a.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.Further lowering D 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 meta-stable (in grey in Figure 5a): in this regime, the system reaches an equilibrium with high or low abundances depending on the initial conditions.It exhibits hysteresis (Figure 5b).
It was recently shown that a metacommunity subject to demographic noise and constant mutualistic interactions exhibits a similar discontinuous phase transition Figure 6: Thermal averaged interaction term, Int = ⟨ j α ij N j ⟩, averaged over non extinct species (indicated by an overline with a + superscript), for two temperatures corresponding to the discontinuous regime.Left: analytical results for T = 0.4, ρ = 1 (as in Figures 5b-d).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 (T = 0.18, D/D 0 (T ) = 0.8, S = 200, L = 400, t max = 500, averaged over 2 runs).Non extinct species are highlighted in orange, only species with negative (or close to zero) interaction terms manage to survive.Averaging only over non extinct species (orange dotted line) leads to a significantly lower (more mutualistic) value than averaging over all species (blu dotted line).µ = 1, σ = 0.5 [45].The authors of [45] 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 [72] (Appendix I).Another quantity of interest is the average interaction term for non-extinct species, Int + = j α ij ⟨N j ⟩ + (the + indicates that the average is carried out only over nonextinct species, ⟨N i ⟩ > 0), which we find to be negative in the entire metastability region (Fig. 6a).In order 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. 6b), thus leading to an enhancement of mutualism between surviving species -see Appendix I for details.
In Figure 7 we also show the phase diagram in the case 0.15 0.16 0.17 The phase diagram for independent interactions across patches (ρ = 0).The continuous line indicates the continuous transition, the dotted and dashed lines the limits of the metastability region, highlighted in grey.At the two limits of the metastability region one of the two solutions disappears and we have a discontinuous transition.µ = 1, σ = 0.5 of independent (ρ = 0) interactions across patches, to be compared to the one of Figure 5a 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, D 0 (T ).For ρ = 1 these two lines coincide, whereas for ρ = 0 the metastability region extends above D 0 (T ) in some range of temperature.In the part of the metastability region above D 0 (T ) the two metastable solutions are both finite: one is of order one and the other proportional to the distance from D 0 (T ); the two solutions coalesce at the tip of the metastability region.
One can also analytically show that the phase diagrams remains qualitatively unchanged considering a small asymmetry in the interactions (γ = 1 − ϵ, ϵ ≪ 1), see Appendix D. 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 be in principle 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 ( 1) and 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 J.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 K), 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.Generically, for moderate system sizes (S < 100 and L < 100) we find strong fluctuations due to the quenched disorder in the interaction matrix, and quantitative finite size effects compared to the asymptotic S, L → ∞ 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 S and L, e.g. S = 200, L = 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 figure 8 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.In order 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 selfsustained solution, if it exists, whereas the latter to the "all-extinct" solution.
We find that indeed at higher temperatures, T = 0.8, in agreement with the analytics and the phase diagram in Figure 5a, the final abundances vary continuously when varying the diffusion constant, and they converge to the same value, no matter the initial condition.The value of D 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 T = 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 [73].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 Figure 8 and Appendix K).
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 D 0 (T ), and there is no strong dependence on the initial conditions (Appendix K). Figure 9: Average abundance ⟨N ⟩ as a function of the diffusion constant D for T = 0.18 and T = 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 t max ).µ = 1, σ = 0.5, S = 200, L = 400, t max = 500 (left) and 200 (right).
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 L = 400, S = 200.Also in this case at T = 0.8 we find a continuous transition and no strong dependence on the initial conditions, while at T = 0.18 we find a discontinuous transition and a hysteresis region (Figure 9) [74].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.
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 selfsustained 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 transition is what is called in ecology, in environmental and social sciences a tipping point or regime-shift [75,76] and in physics 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 done for several other tipping points [77,78], 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 probe of closeness to the tipping point.We can show analytically (see Appendix H) 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 response, which can be measured in controlled lab experiments, does diverge approaching the discontinuous transition, see Figure 10 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).

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 behaviour 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 value of the diffusion constant coincide).This is because when the abundances tend to zero the interactions become sub-dominant 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, that 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 selfsustain 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 a (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 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 system with a finite number of species and patches.
There are several directions worth 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 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 finite dimensional fluctuations destroy the metastable region, and (2) whether the continuous transition can still be described in terms of directed percolation or interactions, although secondary, can alter its universality class.It would also be worth analysing 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, species rich LV model with heterogeneous and strong interactions display multiple equilibria and chaotic dynamics [2][3][4][5].The possibility of different patches to converge to different stationary states could strongly modify the behaviour of the system; in particular allowing the system to experience higher values of the global diversity, possibly violating May's bound [47].
to which we have added a perturbation to the carrying capacity ζ i,u and an external immigration λ i,u , that will be taken to zero at the end of the computation.These equations (for a given value of the η u i (t)) define the trajectories N i,u (t).We add a new species, i = 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 S, 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 S species will only be slightly modified.We consider their linear response: We have introduced the response function R u,v i,j (t, t ′ ) of the abundance of species i in patch u at time t to a variation in the carrying capacity of species j in patch v at time t ′ .
The dynamics of species 0 will depend on these new trajectories: 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 1 − ρ 2 .We thus introduce the matrix a i,j and a u i,j such that and all other correlations are 0. We can rewrite the interaction term as: We want to describe its statistical properties in the limit S → ∞.The response function R u,v i,j (t, t ′ ) is defined on the unperturbed trajectories, and is therefore uncorrelated from the interactions coefficients with species 0. R u,v i,j (t, t ′ ) ∼ 1/ √ S for i ̸ = j [61], so that the off-diagonal terms can be neglected.Thanks to the central limit theorem, j a 0j a j0 R u,v j,j (t, t ′ ) will converge to its average: By similarly evaluating all terms in (A4) we obtain: where ξu (t) and ψu (t) are Gaussian fields with zero mean and covariance E ξu (t) ξv (t ′ ) = E N 0 j,u (t)N 0 j,v (t ′ ) , E ψu (t) ψv (t ′ ) = δ uv E N 0 j,u (t)N 0 j,u (t ′ ) .Note that ξu and the first integral of A6 derive from the component of the interactions constant across patches, a ij , as we can see from the ρ-dependent prefactors, and they therefore couple different patches.ψu and the second integral of A6 derive instead from the component of the interactions independent across patches, a u ij , and therefore represent diagonal correlations and responses.Plugging this expression in the is the Boltzmann distribution with the effective Hamiltonian: 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 reference [70], for a simplified dynamics, that contains all the crucial ingredients: Its equilibrium distribution is given by: where β = 1/T , 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 The MSRDJ action can be written in terms of a deterministic and a dissipative part The time reversal transformation for the two fields is given by: The deterministic and dissipative part of the action are independently invariant under this transformation: where we have defined If the introduction of a small asymmetry in the interactions (ϵ = 1 − γ ≪ 1) is a non-singular perturbation, all the self-consistently determined quantities in the action (h, C d , R int 0 and R d ) will be close to their counterparts for γ = 1.At first order in ϵ we can neglect their change; therefore R d and C d 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: An average ⟨f (N t )⟩ can be expanded as: where ⟨•⟩ 0 indicates the average with respect to the action neglecting δS.
We want to estimate the scaling of 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: The correlation function ⟨f (N t )N v ⟩ 0 will contain a v independent part (that we can neglect since we will be taking the derivative in v) and a connected component of order 1 that decays with the same relaxation time τ .Perturbing the system with a field ζ u this observable will respond as: Inserting these scalings in equation D5 we obtain: 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 N * reads: N * = 0 is always a solution of this equation, we want to find the value of z at which it becomes unstable.We can separate the effective Hamiltonian into a quadratic and a logarithmic part: For N * → 0 the logarithmic part gives rise to a non-integrable divergence in 0 in the denominator.To improve the numerical stability of our solution at small N * we performed an integration by parts of the denominator: The integral is now finite for N * → 0 and we can expand eq.(E1) in powers of N * : The term of order N * (z) 2 is always negative, therefore the number of solution depends on the coefficient of the N * (z) term: if c 1 (z) < 1 the only solution is N * (z) = 0, if c 1 (z) > 1 the N * (z) = 0 solution is unstable and there is a positive stable one.We define the effective growth rate r and the effective growth factor The extinction threshold z * (Figure 11) is given by: As noted in reference [16], this is the same condition that would determine the criticality of the Directed Percolation process with corresponding growth rate and growth factor: The diversity, given by the fraction of non extinct species, can be obtained as: Appendix F: Continuous transition point At the continuous transition, all moments of N tend to zero, and we can expand the extinction condition (E5) in powers of these moments.Keeping only the zeroth order we obtain an equation on the critical value of the diffusion constant: This condition has no dependence on the distribution of the interactions, indeed it is the same that would be obtained with zero or constant interactions [16].
For T → 0 we can expand Eq. (F1) and show that D 0 vanishes exponentially: The reason for this behavior is that at low demographic noise, the abundances of a species with carrying capacity k undergoes a fluctuation toward very low values very rarely.In fact, one needs to wait 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 e k 2 2T (the "energy barrier" equals k 2 /2).The equation above can be therefore interpreted as a balance between two inverse time-scales: 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 h ∝ q d ∝ D − D 0 , q 0 ∝ (D − D 0 ) 2 , and z * (and therefore ϕ) has a finite limit.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 H ef f , we will call this P (N |z).On top of this, because of disorder, different species experience different average interactions with the rest of the ecosystem (different values of z, that is distributed according to P (z) Gaussian), leading to species-dependent factors in H ef f .If we want to study the distribution of the abundances of all species at a given time in one site (P (N )) we need to take into account both effects.We can compute P (N ) marginalizing over z: with We could also be interested in the distribution across species of the abundance averaged over patches or time, given by P (N * ) = P (z) dN * (z) dz −1 for N * > 0. There is also a finite probability 1−ϕ that N * = 0, where ϕ is the diversity.Examples of these abundance probability distributions are shown in Fig. 12.

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 [77].Let us consider a generic dynamical system, described by ⃗ x contains all the degrees of freedom of the system, whereas k is a control parameter.The zeros of F yield the stationary states x * : The stationary point is stable if the Jacobian of F has only negative eigenvalues, ensuring that x returns to x * 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 k, 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 k involves some connected moments of N and the derivative of H in k: (H13) We collect the three order parameters in a vector ⃗ p = (h, C 0 d , q 0 ) T .Then d⃗ p dk satisfies: where Ĵ is a 3 × 3 matrix and s a vector; their elements are the coefficients of Equations (H11-H13).The solution is given by: The response to a variation of k diverges if Ĵ − 1 has a zero eigenvalue, this is found to happen when approaching the discontinuous transition.We expect the same qualitative behaviour of the response to perturbations for generic values of ρ, but for ρ ̸ = 0 we need to take into account also the variations of the function N * (z), 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 [72]: Figure 13: Left: Analytical estimate of the mean of the reduced interaction matrix using Eq.(I1) for T = 0.4.Right: Numerical results for the distribution of the interaction coefficients for all and surviving species.T = 0.18, D/D 0 = 0.8.In both cases in the initial species pool µ = 1, σ = 0.5; if all species go extinct we say µ ′ = 0.
Since we are not at 0 temperature, in our case this formula is only an approximation, but it provides an useful estimate of the variation of the mean interaction.We find that the interaction mean decreases (more mutualistic) when decreasing the diffusion coefficient (Figure 13a); it is negative in the entire metastability region.In Figure 13b 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 µ ′ = S 1 S 2 ij α ij changes from 0.96 to -0.28.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 Equation A6 We can now average it over all species (all values of z, overline), or over only non extinct ones (z < z * , overline with + superscript).
Note that we will always find I + < I; I + is negative in the entire metastability region (Figure 6a 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 non extinct ones (Figure 6b 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 reference [59], and improved in [5,81].The idea is to separate the process in a deterministic part: and a stochastic one: Ṅi,u = N i,u η i,u .(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: Ñi,u (t) = Gamma P oisson N i,u (t) T dt T dt .

(J3)
For the deterministic part we rely on Euler method.

Figure 1 :
Figure 1: A metacommunity of 7 species living on 3 patches.Each individual interacts with the local community to which it belongs possibly migrating to neighboring patches with diffusion coefficient D. time

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

Figure 4 : 1 0 8 Figure 5 :
Figure 4: Average abundance ⟨N ⟩ and diversity ϕ as a function of the diffusion constant D for T = 0.25 (top) and as a function of temperature (strength of demographic noise) for D = 0.1 (bottom).The dashed lines represent the T = 0 well-mixed results.µ = 1, σ = 0.5.

Figure 7 :
Figure 7: The phase diagram for independent interactions across patches (ρ = 0).The continuous line indicates the continuous transition, the dotted and dashed lines the limits of the metastability region, highlighted in grey.At the two limits of the metastability region one of the two solutions disappears and we have a discontinuous transition.µ = 1, σ = 0.5

8 Figure 8 :
Figure 8: Average abundance ⟨N ⟩ as a function of the diffusion constant D for T = 0.18 and T = 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, S = 200, L = 400, t max = 500 (left) and 200 (right).

Figure 10 :
Figure 10: Average abundance and its response to a perturbation of the carrying capacity k at T = 0.153 for ρ = 0 approaching the instability of the self-sustained phase.µ = 1, σ = 0.5.

Figure 12 :
Figure 12: (a) Probability distribution of the abundance N ; in orange deep in the survival phase (T = 0.8, D/D 0 (T ) = 1.5), in blue right before the discontinuous transition (T = 0.4, D/D 0 (T ) = 0.84).Note that the shown distributions do not integrate to 1 because a finite fraction of the species are extinct, leading to a delta function in zero with weight 1 − ϕ.(b) Probability distribution of the abundance for a given species, i.e. at fixed z; in blue for a species close to extinction (z = −0.45,N * = 0.148), in orange for a species far from extinction (z = −2.45,N * = 1.733);T = 0.4, D/D 0 (T ) = 0.84.(c) Probability distribution of the space (or time) averaged abundance, because of extinct species we again have a delta function in zero with weight 1 − ϕ.T = 0.4, D/D 0 (T ) = 0.84.µ = 1, σ = 0.5, ρ = 1.