Diverse communities behave like typical random ecosystems

TL;DR

This paper demonstrates that diverse ecological communities often behave like random ecosystems, especially when noise is added to consumer preferences, explaining the success of random models in microbial ecology.

Contribution

The study extends May's random matrix approach to generalized consumer-resource models, showing that diverse communities exhibit typical random ecosystem behavior under certain conditions.

Findings

Adding noise leads to a transition to typicality in community properties.

Random consumer-resource models effectively reproduce microbial community patterns.

Diverse communities are difficult to model with bottom-up approaches.

Abstract

In 1972, Robert May triggered a worldwide research program studying ecological communities using random matrix theory. Yet, it remains unclear if and when we can treat real communities as random ecosystems. Here, we draw on recent progress in random matrix theory and statistical physics to extend May's approach to generalized consumer-resource models. We show that in diverse ecosystems adding even modest amounts of noise to consumer preferences results in a transition to "typicality" where macroscopic ecological properties of communities are indistinguishable from those of random ecosystems, even when resource preferences have prominent designed structures. We test these ideas using numerical simulations on a wide variety of ecological models. Our work offers an explanation for the success of random consumer-resource models in reproducing experimentally observed ecological patterns in microbial communities and highlights the difficulty of scaling up bottom-up approaches in synthetic ecology to diverse communities.

Equations

$$\{\begin{array}{l}\frac{d{N}_i}{dt}=N_i({\sum }_{\beta }{\bar{C}}_{i\beta }{R}_{\beta }-m_i)\\ \frac{d{R}_{\alpha }}{dt}=R_{\alpha }(K_{\alpha }-R_{\alpha }-{\sum }_j{N}_j{\bar{C}}_{j\alpha }).\end{array}$$

$$\{\begin{array}{l}\frac{d{N}_i}{dt}=N_i({\sum }_{\beta }{\bar{C}}_{i\beta }{R}_{\beta }-m_i)\\ \frac{d{R}_{\alpha }}{dt}=K_{\alpha }-R_{\alpha }-{\sum }_j{N}_j{\bar{C}}_{j\alpha }.\end{array}$$

$$\bar{\mathbf{C}}=\mathbf{B}+\mathbf{C},$$

$$\frac{d{N}_i}{dt}=N_i\left(\sum _{\alpha \in \mathbf{M}}{C}_{i\alpha }{K}_{\alpha }-m_i-\sum _j{A}_{ij}{N}_j\right),$$

$$A_{ij}=\sum _{\alpha \in \mathbf{M}}{\bar{C}}_{i\alpha }{\bar{C}}_{\alpha j}^T$$

$${\lambda }_{\min }={\sigma }_c^2(1-\sqrt{S^{\ast }/M}{)}^2.$$

$$\nu =\frac{1}{M}\sum _i(1/{\lambda }_i)=\frac{1}{M}\mathrm{tr}({\mathbf{A}}^{-1}).$$

$$A_{ij}=\sum _{\alpha \in {\mathbf{M}}^{\ast }}{\bar{C}}_{i\alpha }{\bar{C}}_{\alpha j}^T$$

$$\bar{\mathbf{C}}=\mathbf{B}+\mathbf{C},$$

$$\{\begin{array}{l}\frac{d{N}_i}{dt}=N_i({\sum }_{\beta }{\bar{C}}_{i\beta }{R}_{\beta }-m_i)\\ \frac{d{R}_{\alpha }}{dt}={\kappa }_{\alpha }-R_{\alpha }-{\sum }_j{N}_j{\bar{C}}_{j\alpha }{R}_{\alpha }.\end{array}$$

$${\chi }_{\alpha \beta }^R=\frac{\partial {\bar{R}}_{\alpha }}{\partial K_{\beta }},{\chi }_{i\alpha }^N=\frac{\partial {\bar{N}}_i}{\partial K_{\alpha }},{\nu }_{\alpha i}^R=\frac{\partial {\bar{R}}_{\alpha }}{\partial m_i},{\nu }_{ij}^N=\frac{\partial {\bar{N}}_i}{\partial m_j}$$

$$0$$

$$0=\sum _{\alpha \in {\mathbf{M}}^{\ast }}{\bar{C}}_{i\alpha }\frac{\partial {\bar{R}}_{\alpha }}{\partial K_{\beta }},$$

$${\delta }_{ij}=\sum _{\alpha \in {\mathbf{M}}^{\ast }}{\bar{C}}_{i\alpha }\frac{\partial {\bar{R}}_{\alpha }}{\partial m_j},$$

$$0=\sum _{\alpha \in {\mathbf{M}}^{\ast }}{\bar{C}}_{i\alpha }{\chi }_{\alpha \beta }^R,$$

$${\delta }_{ij}=\sum _{\alpha \in {\mathbf{M}}^{\ast }}{\bar{C}}_{i\alpha }{\nu }_{\alpha j}^R,$$

$$\left(\begin{array}{cc}\bar{\mathbf{C}}& 0\\ \text{mathds}1& {\bar{\mathbf{C}}}^T\end{array}\right)\left(\begin{array}{cc}{\nu }^R& {\chi }^R\\ {\nu }^N& {\chi }^N\end{array}\right)=\text{mathds}1$$

$${\chi }_{\alpha \beta }^R$$

$${\chi }_{i\alpha }^N$$

$${\nu }_{ij}^N$$

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

$$\frac{d{R}_{\alpha }}{dt}$$

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

$$\frac{d{R}_0}{dt}$$

$${\bar{N}}_i$$

$${\bar{R}}_{\alpha }$$

$$z_N$$

$$z_R$$

$${\bar{R}}_0=\text{max}\left[0,\frac{{\sigma }_c^2\chi (K-\mu ⟨N⟩+z_R)-\mu ⟨R⟩+m-z_N}{(1-{\sigma }_c^2\nu ){\sigma }_c^2\chi +1}\right]$$

$${\bar{N}}_0=\text{max}\left[0,\frac{(1-{\sigma }_c^2\nu )(\mu ⟨R⟩-m+z_N)+K-\mu ⟨N⟩+z_R}{(1-{\sigma }_c^2\nu ){\sigma }_c^2\chi +1}\right]$$

$${\bar{R}}_0=0,{\bar{N}}_0=\frac{\mu ⟨R⟩-m+z_N}{{\sigma }_c^2\chi }\text{or,}{\bar{N}}_0=0,{\bar{R}}_0=\frac{K-\mu ⟨N⟩+z_R}{1-{\sigma }_c^2\nu }$$

$${\bar{R}}_0=0\text{ and }{\bar{N}}_0=0.$$

$${\bar{R}}_0=\Theta (R_0)\left[\Theta (N_0)\frac{{\sigma }_c^2\chi (K-\mu ⟨N⟩+z_R)-\mu ⟨R⟩+m-z_N}{(1-{\sigma }_c^2\nu ){\sigma }_c^2\chi +1}+(1-\Theta (N_0))\frac{K-\mu ⟨N⟩+z_R}{1-{\sigma }_c^2\nu }\right]$$

$${\bar{N}}_0=\Theta (N_0)\left[\Theta (R_0)\frac{(1-{\sigma }_c^2\nu )(\mu ⟨R⟩-m+z_N)+K-\mu ⟨N⟩+z_R}{(1-{\sigma }_c^2\nu ){\sigma }_c^2\chi +1}+(1-\Theta (R_0))\frac{\mu ⟨R⟩-m+z_N}{{\sigma }_c^2\chi }\right]$$

Full text

Stability $|$ Complex Ecosystem $|$ Phase Transition $|$

I Introduction

One of the most stunning aspects of the natural world is the immense diversity of ecological communities ranging from rainforests to human microbiomes. Ecological communities are critical for numerous processes ranging from global water cycling processes[1] to animal development and host health[2]. For this reason, understanding the principles governing community assembly and function in diverse communities has wide ranging applications from conservation efforts to pharmaceutical engineering and bioremediation[3].

Many traditional ecological models focus on ecosystems consisting of a few species and resources. In such low dimensional models, it is often possible to characterize the ecological traits of all the species and resources and then use this information to make predictions about community-level properties [4, 5, 6]. However, many natural communities are extremely diverse and the models and parameters are naturally high dimensional. This problem is especially pronounced in in the context of microbial ecology where hundreds of species can coexist in a single location. In this case, a comprehensive parametrization of species and resource traits is no longer feasible, suggesting that new ideas and concepts are required to understand diverse communities.

A similar problem is encountered in statistical physics. For example, an ideal gas is characterized by the unit mole, which has the order of $10^{23}$ particles, making it impossible to simultaneously specify the microscopic state of the system (e.g. the positions and velocities of all particles). Despite this uncertainty, it is still possible to make predictions about macroscopic properties like pressure and the average energy by treating the positions and velocities of particles as independent random variables[7]. The fact that such universal statistical behaviors emerge naturally in large disorder systems composed on many particles suggests that a similar approach maybe possible in ecological systems.

In 1972, Robert May suggested that large complex ecosystems can also be modeled as random systems [8]. May considered a diverse ecosystem composed of $S$ species whose interspecific interactions were sampled randomly and independently from a normal distribution with zero mean and variance $\sigma^{2}$. In particular, May asked when such a diverse random ecosystem would be stable to small perturbations. To answer this question, he examined the largest, i.e., the rightmost eigenvalue $\lambda_{\rm max}$ of the $S\times S$ community interaction matrix $\mathbf{J}$, whose diagonal entries – chosen to be $J_{ii}=-1$ by May– describe intraspecific competition and off-diagonal entries $J_{ij}$ describe how much the growth rate of species $i$ is affected by a small change in the population $N_{j}$ of other species $j$ from its equilibrium value. Using a mathematical formula for the distribution of eigenvalues of large random matrices derived by Ginibre [9], May showed that $\lambda_{\rm max}$ increases with $S$, and derived a stability criterion governing the maximum diversity of an ecosystem: a diverse ecosystem becomes unstable to small perturbations when $\sqrt{S}\sigma>1$[8]. May’s stability criteria has proven to be robust against a wide array of changes in the assumptions, including adding biologically realistic correlation structures to the matrix, or incorporating the dependence of the community matrix on population sizes in the Lotka-Volterra model [10, 11].

In May’s model, all ecosystem properties are encoded in the species-species interaction matrix. A major limitation of these models is that they neglect resource dynamics, making it difficult to understand how ecosystem properties depend on both the external environment and species consumer preferences. For this reason, community assembly is often analyzed using generalized Consumer Resource Models (CRMs)[12, 13]. In these models, species are modeled as consumer that can consume resources, and sometimes also produce resources [14, 15, 16, 17, 18]. Recently, we have shown that such models, initialized with random parameters, can predict lab experiments on complex microbial communities [14, 16] and reproduce large-scale ecological patterns observed in field surveys, including the Earth and Human Microbiome Projects [16]. This suggests that the large-scale, reproducible patterns we see across Microbiomes are emergent features of random ecosystems.

Yet, it remains unclear why random ecosystems can accurately describe real ecological communities. To answer these question, in this paper we exploit ideas from random matrix theory and statistical physics to analyze generalized consumer-resource models in spirit of May’s original analysis. We show that the macroscopic ecological properties of diverse ecosystems can be described using random ecosystems, much like thermodynamic quantities like pressure and average energy of the ideal gas can be described by considering particles to be random and independent.

I.1 Models

To explore these ideas, we devised a more concrete version of May’s original thought experiment describing an ecosystem consisting of $S$ non-interacting species where interactions are gradually turned on. May’s original argument only considered the local dynamics near a pre-specified equilibrium point that eventually becomes unstable. Since we are interested in exploring what happens in consumer resource models, we must make additional modeling assumptions to arrive at a complete set of nonlinear dynamics. We focus on numerous variants of the Consumer Resource Model (CRM)[12], including different choices of resource dynamics, consumer preferences, as well as more dramatic variants such as the Microbial Consumer Resource Model introduced in [14, 15, 16].

The original MacArthur Consumer Resource Model [12] consists of $S$ species or consumers with abundances $N_{i}$ ($i=1...S$) that can consume one of $M$ substitutable resources with abundances $R_{\alpha}$ ($\alpha=1...M$), whose dynamics are described by the equations

[TABLE]

The consumption rate of species $i$ for resource $\alpha$ is encoded by the entry $\bar{C}_{i\alpha}$ in the $S\times M$ consumer preference matrix $\mathbf{\bar{C}}$, $K_{\alpha}$ is the carrying capacity of resource $\alpha$, and $m_{i}$ is a maintenance energy that encodes the minimum amount of energy that a species $i$ must harvest from the environment to survive. When the system is in the steady state, some species and resources can vanish. We denote the numbers of surviving species and resources by $S^{*}$ and $M^{*}$, respectively, and in general at steady state we will have $S^{*}\leq S$ and $M^{*}\leq M$. For this reason, we refer to this model as the CRM with resource extinction and consider its effects analytically and numerically in Section II.3 and Appendix C.1.

In the beginning, we focus primarily on a popular variant of the original CRM introduced by Tilman with slightly different resource dynamics[13]:

[TABLE]

From an ecological perspective, there are significant differences between this model variant and the original CRM. First, the resource supply rate $K_{\alpha}$ is constant instead of following logistic growth rate. Second, the species consume resources at a rate that is independent of the resource concentrations in the environment. This can lead to unphysical, negative resource concentrations. Despite these differences, mathematically the equilibrium solutions of the two models have similar forms. One major difference that does arise is that in the dynamics described by Eq. 2 consumers can no longer cause a resource to go extinct (i.e., $\mathbf{M^{*}}=\mathbf{M}$). This makes this models significantly easier to analyze (especially within the context of Random Matrix Theory) and leads to much simpler analytic expressions. For this reason, we largely focus on this latter model without resource extinction (see Fig. A9, Fig. A10, Fig. A8 for dynamics described by Eq 2 and Appendix J for numerics and Appendix C for analytics on original CRM described by Eq. 1). Despite the unphysical, negative resource concentrations, the CRM without resource extinction captures almost all the qualitative behaviors present in more complicated and physically realistic CRMs (though there are some subtle but important differences discussed below).

Both the models in Eq. 1 and Eq. 2 make very specific assumptions about resource dynamics. To check the generality of our results, we also numerically analyzed generalizations of the CRM including linear resource dynamics where resources are supplied externally, and a model of microbial ecology with trophic feedbacks where organisms can feed each other via metabolic byproducts[14, 19, 16, 15]. This analysis can be found in Appendix A. Furthermore, for simplicity, in most of this work we assume that $S=M$. However, we have numerically checked that our results are robust to breaking on this assumption (see Fig. A6).

In CRMs, the identity of each species is specified by its consumption preferences. In real ecosystems, it is well established that organisms can exhibit strong consumer preferences for particular resources. However, recent work has shown that consumer resource models with random consumer preferences can reproduce experimental observations in field surveys and laboratory experiments [14, 16]. To understand this phenomena, we asked how adding noise to consumer preferences changes macroscopic ecosystem level properties like diversity and average productivity. To do so, we considered a thought experiment where we started with predesigned consumer resource preference, and then added “noise” to the consumer resource preferences. Mathematically, we can decompose the consumer matrix $\bar{\mathbf{C}}$ in Eqs 1 and 2 into two parts:

[TABLE]

where $\mathbf{B}$ encodes pre-designed structures, and $\mathbf{C}$ is a random matrix representing ”noise”.

For simplicity, we started with non- interacting species where each species consumes its own resource. A set of non-interacting species can be constructed by engineering each species to consume a different resource type, with no overlap between consumption preferences. For example, one can imagine designing strains of E. coli where each strain expresses transporters only for a single carbon source with all other transporters edited out of the genome: i.e a strain that can only transport lactose, another strain that can only transport sucrose, etc. An ecosystem with such consumer preference structure is shown in Figure 1(A). In such an experiment, horizontal gene transfer would eventually begin distributing transporter genes from one strain to another, so a realistic model would have to allow for some amount of unintended, “off-target” resource consumption. In line with May, we can model the consumer preferences $\bar{C}_{i\alpha}$ of species $i$ for resource $\alpha$ in such an ecosystem as the sum of the identity matrix $\mathbf{B}=\mathds{1}$ and a random component $C_{i\alpha}$ with variance $\sigma^{2}$ that encodes non-specific preferences (see Figure 1A right). In other words, the full consumer matrix can be written as $\bar{\mathbf{C}}=\mathbf{I}+\mathbf{C}$.

II Results

II.1 Phase transition to random ecosystems

Figure 1(B) shows how the number of surviving species at steady-state changes as one adds more and more non-specific resource preferences to an ecosystem initially composed of non-interacting species. Just as in May’s analysis, the appropriate measure of the importance of the random component is the root-mean-squared off-target consumption $\sigma_{c}=\sqrt{M\sigma^{2}}$ (recall $M=S$). This scaling reflects the fact that two consumer matrices $\bar{\mathbf{C}}$ with the same $\sigma_{c}$ but different system sizes $M$ can have very different amounts of absolute noise as shown Figure 1(C), but exhibit almost identical community-level properties (with all differences coming from finite size effects, see Fig. A4 for the universal behavior at different $M$). Figure 1(B) shows the fraction of surviving species $S^{*}/M$ in the ecosystem as a function of $\sigma_{c}$. At small values of $\sigma_{c}$, all the species survive and $S^{*}=S$. As high as $\sigma_{c}=0.7$, almost all of the original species are still present in the community. But between $\sigma_{c}=0.7$ and $\sigma_{c}=1$, there is a sharp transition in community structure, which results in about half of the original species becoming extinct.

Remarkably, the fraction of surviving species converges to the same value as for a completely random consumer preference matrix and remains finite as $\sigma_{c}\to\infty$ [20]. This means that ecosystems with an arbitrarily large number of species can be stably formed by considering a sufficiently large initial species pool. We also examined two other community-level properties: the mean species abundance $\langle N\rangle$ (i.e., the average productivity), and the second moment of the population size $\langle N^{2}\rangle$, which includes information about the distribution of population sizes of various species. Figure 2 shows that both of these quantities are also well-approximated by the random consumer preference matrix for $\sigma_{c}>1$. These numerical predictions are in excellent agreement with analytic predictions derived in the $S\to\infty$ limit derived in Appendix C using the cavity method [21, 22].

This convergence to random ecosystem behavior is quite robust, and holds for other choices of designed consumer preferences beyond the identity matrix considered above. Figure 2 shows numerical simulations of the diversity $S^{*}/M$, average productivity $\langle N\rangle$, and second moment of the species abundances $\langle N^{2}\rangle$ as a function of the noise $\sigma_{c}$ for two other choices of designed consumer preference matrices: a block structure with pre-defined groups of species exhibiting strong intra-group competition and a unimodal structure where each species is more likely to consume resources similar to its preferred resource. Once again, we see that the ecosystem quickly transitions to a behavior where these macroscopic properties are indistinguishable from those of a random ecosystem. Borrowing terminology from physics, we call systems whose macroscopic properties are well described by random ecosystems as typical. The primary effect of the choice of consumer preference matrix is to adjust the threshold value of $\sigma_{c}$ where the transition to typicality occurs. In all cases, we find that the random behavior takes over when the average total off-target consumption capacity over all $M$ resource types becomes greater than the consumption of the primary resource in the original designed ecosystem in the absence of noise.

The character of the self-organized state is also robust to changes in the sampling scheme for the random component of the consumer preferences. Gaussian noise in consumer preferences simplifies the analytic calculations but also sometimes results in non-physical negative values for consumer preferences. We therefore tested two sampling schemes that always produce positive values for consumer preferences: uniformly sampling the random component of preferences $C_{i\alpha}$ in an interval from 0 to $b$, and binary sampling where $C_{i\alpha}=1$ with probability $p_{c}$ and zero otherwise. Changing $b$ or $p_{c}$ affects both the mean and the variance of the random components of the consumer preferences simultaneously making it difficult to directly compare to the Gaussian case. Nonetheless, as can be seen in the Figure 2, the qualitative behaviors is identical to the Gaussian case, with macroscopic ecological properties becoming indistinguishable from those of a fully random ecosystem when the average off-target resource consumption comparable to the the consumption of the designed resources.

II.2 Sensitivity to perturbations and the transition to typicality

To better understand why mass extinctions happen at $\sigma_{c}^{*}\sim 1$ and allow for comparison with May’s original analysis, we calculated an effective species-species competition matrix $A_{ij}$ between species for an ecosystem whose dynamics are governed by Eq. 2. We exploited the observation by MacArthur and others that if resource abundances always remain close to their steady state values, the steady-states of the CRM coincide with those of an effective generalized Lotka-Volterra model of the form

[TABLE]

with the species-species interaction matrix given by

[TABLE]

(see Figure 3(A) and Appendix D for details). This matrix is related to May’s community matrix governing stability $\mathbf{J}$ discussed in the introduction through the relation $J_{ij}=-\bar{N}_{i}A_{ij}$, where $\bar{N}_{i}$ is the steady-state abundance of species $i$. For symmetric interaction matrices of the form in Eq. 4, it is possible to prove that the largest eigenvalue $\lambda_{\rm max}$ of $\mathbf{J}$ reaches zero from below only when the smallest eigenvalue $\lambda_{\rm min}$ of $\mathbf{A}$ reaches zero from above (see Appendix E).

As shown in Figure 1(B), the behavior broadly falls into one of three different regimes depending on the amount of noise introduced in the consumer preferences: a low-noise regime when $\sigma_{c}\ll 1$, a cross-over regime when $0\ll\sigma_{c}\leq 1$, and a high-noise regime when $\sigma_{c}>1$. Figure 3(B) shows how the eigenvalue spectrum of the corresponding Lotka-Volterra interaction matrix $\mathbf{A}$ change as $\sigma_{c}$ increases.

Low-noise regime ($\sigma_{c}\ll 1$): In the low-noise regime, the engineered structure in the consumer preference controls large scale ecological properties. Furthermore, the eigenvalue spectrum of the LV-interaction matrix $\mathbf{A}$ is centered around $1$ reflecting the fact there is very little competition between species (i.e., species still occupy largely independent niches). For this reason, in this regime all the initial species in the ecosystem survive to steady-state so that ${S^{*}}/{M}=1$.

Crossover regime ($0\ll\sigma_{c}\leq 1$): With increasing $\sigma_{c}$, the eigenvalues due the noise component in $\mathbf{A}$ repel each other like in the Coulomb gas and the spectrum spreads out [23]. $\lambda_{\rm min}$ decreases until it reaches the threshold of stability $\lambda_{\rm min}\cong 0$ at $\sigma_{c}^{*}\approx 1$. Note that $\lambda_{\rm min}$ is close to 0 but not exactly at 0 because the steady-state of the CRM is always stable [24]. In this regime even a small environmental perturbations or small amounts of demographic noise can result in species extinctions [25]. This is closely related to the divergence of structural stability when $\lambda_{\rm min}\sim 0$[26]. In Appendix C we show analytically using the Cavity method [21, 22] that in the limit $M\to\infty$, $\lambda_{\rm min}$ is approaches 0 from above when $\sigma^{*}_{c}=1$. At $\sigma_{c}\sim 1$ the engineered structure and noise have comparable amplitudes. For the case where the consumer preferences are chose to be binary noise, this threshold corresponds to a critical noise level $p_{c}\sim\frac{1}{M}$, meaning on average there is one random nonzero element in the row besides the diagonal one. More generally, our numerics suggest that the threshold to typicality occurs in a wide variety of models when the expected off-target resource consumption rates become comparable to the consumption rate for the designed resources.

Noise-dominated regime ($\sigma_{c}>1$) In this regime, we observe two new phenomena that were not accessible in May’s original framework. First, the spectrum of the species-species interaction matrix $A_{ij}$ approaches the Marchenko-Pastur law [27],

[TABLE]

where $a=\sigma_{c}^{2}(1-\sqrt{c})^{2}$, $b=\sigma_{c}^{2}(1+\sqrt{c})^{2}$, $c={S^{*}}/{M}$ and $\Theta(x)$ represents the Heaviside step function. This differs from May’s analysis where the spectrum of the interaction network follows Girko’s Circular law [28, 29, 30]. The reason for this difference is that species-species interaction matrix obtained from the CRM is the outer product of a random matrix $\mathbf{\bar{C}}$ with itself (i.e., a Wishart matrix, see Eq. 4), reflecting the fact that the CRM has two different kinds of degrees of freedom: resources and species. The Marchenko-Pastur law is the distribution we would expect for an ecosystem with completely random consumer preferences [27]. This helps explain our earlier observations that community-level observables of ecosystems are indistinguishable from the purely random ecosystems when $\sigma_{c}$ is sufficiently large (see Figure 3(B)).

Secondly, as $\sigma_{c}$ increases past 1 and ecosystem properties become typical, the resulting ecosystems once again become insensitive to external perturbation [25]. To see this, we note that we can measure sensitivity to perturbations by examining the minimum eigenvalue of the interaction matrix $A_{ij}$, with larger $\lambda_{\rm min}$ meaning decreased sensitivity to perturabations (see Appendix E). The minimum eigenvalue in the Marchenko-Pastur Distribution is located at

[TABLE]

As one increases $\sigma_{c}$, $S^{*}/M\rightarrow 1/2$ from above since there is increases competition between species for shared resources. Consequently, $\lambda_{\rm min}$ is always much larger than zero once ecosystems crossover to their typical behavior.

The above analysis suggests that $\lambda_{\rm min}$ is an important property that can be used to characterize the three regimes seen in Figure 3(C). In the low-noise regime, species-species interactions are weak and $\lambda_{\rm min}\approx 1$, whereas in the high-noise regime $\lambda_{\rm min}=\sigma_{c}^{2}(1-\sqrt{S^{*}/M})^{2}$. The calculation of $\lambda_{\rm min}$ in Regime B is challenging because of the mixture between the engineered structure and noise. However, we can use techniques from RMT for wireless communication (i.e information-plus-noise models) to analytically estimate $\lambda_{\rm min}$ [31, 32]. The results are shown in the red scatter points in Figure 3(D) (see Appendix F.3). As discussed above, $\lambda_{\rm min}$ approaches zero as $\sigma_{c}$ approaches one.

The spectrum of $\mathbf{A}$ also contains quantitative information about the sensitivity of the ecosystem in the Cavity method. Specifically, as shown in Appendix C, we can define a susceptibility $\nu$ that measures the average response of the steady-state population size $\bar{N}_{i}$ to perturbing of the species maintenance cost $m_{i}$ (see Eq. 2). We further show that $\nu$ is directly related to the the sum of the inverse eigenvalues of $A_{ij}$ through the expression

[TABLE]

Figure 3(D) shows that this quantity is initially constant as $\sigma_{c}$ is increased from 0, then reaches the maximum value at $\sigma_{c}=1$, and finally rapidly decreases to near zero. In Appendix C we provide analytical calculations based on the cavity method confirming these numerical results.

Note that our results are not restricted to Gaussian noise but also apply to the other cases where the noise in consumer preferences is binary or uniform (See Figure A7 and A9). This is because the central limit theorem guarantees that the statistics of eigenvalues of large random matrices converge to the statistics in Gaussian random matrices for many biologically plausible choices of consumer preferences.

II.3 Effect of resource extinction

Thus far we have focused on a CRM without resource extinctions specified by Eqs. 2. As discussed extensively in Appendix C, if we allow for resource extinction (Eqs. 1) and write

[TABLE]

instead of eq. 4, somewhat surprisingly, our cavity method predicts a second-order phase transition to typicality rather than a cross-over as is the case without resource extinction. The signature of such a second order transition is the divergence of the susceptibility matrix $\nu$ discussed above. Figure 4 shows $\nu$ with and without resource extinction, numerically confirming the existence of this second order transition. This second order transition is also reflected in the spectrum of the interaction matrix $\mathbf{A}$ through the the appearance of zero eigenvalue modes for CRMs when resources can go extinct.

The existence of zero modes can be understood by noting that resource extinction and species extinction correspond to the column and row deletion in the consumption matrix (shown in Figure 4(A)). Such deletions can change the engineered component of the effective consumer preferences for surviving species and resources, resulting in large fluctuations in the interaction matrix $\mathbf{A}$. In the presence of these large fluctuations, the interaction matrix no longer self-averages, giving rise to the observed second-order phase transition. This same mechanism also leads to a second-order phase transition to typical behavior when the engineered portion of the consumer resources is block diagonal, even in the absence of resource extinctions (see Figure A10).

III Discussion

It is common practice in theoretical ecology to model ecosystems using random matrices. Yet it remains unclear if and when we can treat real communities as random ecosystems. Here, we investigated this question by generalizing May’s analysis to consumer resource models and asking when the macroscopic, community level properties can be accurately predicted using random parameters. We found that introducing even modest amount of stochasticity into consumer preferences ensures that the macroscopic properties of diverse ecosystems will be indistinguishable from those of a completely random ecosystem. Our calculations and numerics suggest that transition to typicality occurs when the total amount of off-target resource consumption becomes comparable to the consumption rate of targeted resources.

We confirmed our analytic calculations using numerical simulations on CRMs with different types of resource dynamics and different classes of non-specific interactions. We emphasize that despite the fact that random ecosystems can make accurate predictions about macroscopic properties like the average diversity or productivity, they will in general fail to capture species level details. This phenomena is well understood in the context of statistical physics where it is possible to predict thermodynamic quantities such as pressure and temperature even though one cannot accurately predict microstates.

These observations may help explain the surprising success of consumer resource models with random parameters in predicting the behavior of microbial ecosystems in the lab and natural environments [14, 16]. They also suggest that maybe possible to predict macroscopic ecosystem level properties like diversity or total biomass even when ecosystems are poorly characterized or have lots of missing data.

The foregoing analysis has several other interesting implications. First, it suggests that bottom-up engineering of complex ecosystems may prove to be very difficult. As the number of components increases, small uncertainties in each of the interaction parameters may eventually overwhelm the designed interactions, and destabilize the intended steady state. Instead, such system are much more likely end up in a typical state which our theory suggests is much more stable than the intended designed state as ecosystems become more diverse.

Our work also suggests that in ecosystems well described by consumer resource models, crossing a May-like transition generically gives rise to typical random ecosystems rather than a marginal stable phase as was found in a recent analysis of the Generalized Lotka-Volterra model [33, 34, 35] (an important caveat to this statement is that adding non-resource based interactions to consumer resource models can restore complicated behavior reminiscent of the marginally stable phase [25]). This suggests that even when cumulative parameter uncertainties preclude a detailed characterization of an ecosystem, methods from statistical physics and Random Matrix Theory can be employed to predict system-level properties [36, 20]. It will be interesting to explore if and how these insights can be exploited to design top-down control strategies for ecosystems and identify assembly rules for microbial communities with many species [4].

In this paper, we only consider white noise, which is independently and identically added to all interaction components. In the future, it will be interesting to ask how other specialized noises, resulting from demographic stochasticity, phenotypic variation, can affect our results. Based on our experience, we expect that, even in these more complicated ecosystems, our conclusion will hold quite generally in the thermodynamics limit. But much more work needs to be done to confirm if this intuition is really correct.

Acknowledgements.

We thank Josh Goldford, Zhenyu Liao, Jason Rocks, Guangwei Si, Jean Vila and Yu Hu for helpful discussions. We also especially appreciate numerous valuable comments from Stefano Allesina. The work was supported by NIH NIGMS grant 1R35GM119461, Simons Investigator in the Mathematical Modeling of Living Systems (MMLS). The authors are pleased to acknowledge that the computational work reported on in this paper was performed on the Shared Computing Cluster which is administered by Boston University Research Computing Services.

Appendix A Model setup

We primarily analyze CRMs of the form given by Eqs 1 and 2. To do so, we decompose the consumer matrix $\bar{\mathbf{C}}$ into two parts:

[TABLE]

with $\mathbf{B}$ encoding a pre-designed set of resource-mediated interactions, and $\mathbf{C}$ a random matrix encoding “off-target” consumption. We consider three types of $\mathbf{B}$ (see Figure 2): the identity matrix, a square Gaussian-type circulant matrix $B_{i\alpha}=e^{-\text{min}(i,|M-i|)^{2}/r^{2}}$ with $r=7$[37] and a block matrix with identical $10\times 10$ blocks(all elements are 1 inside the $10\times 10$ block). We also consider three types of random matrices $\mathbf{C}$. In all cases, each element in the matrix is sampled independently from an underlying probability distribution. The three distributions we consider are a normal distribution with mean zero and standard deviation $\sigma_{c}/\sqrt{M}$, a uniform distribution where each element is sampled uniformly from $[0,b]$, and a Bernoulli distribution where each element can be $+1$ with probability $p_{c}$ and [math] with probability $1-p_{c}$ (i.e Binary Noise).

For all simulations, unless otherwise specified the default choices for parameters are: $M=100$, $\mu=0$, $K=1$, $\sigma_{K}=0.1$ , $m=0.1$ and $\sigma_{m}=0.01$ and each data point is averaged from 5000 independent realizations. The simulation detail for each figure can be found at Appendix G. All simulations are available on GitHub at https://github.com/Emergent-Behaviors-in-Biology/typical-random-ecosystems.

A.1 Alternative Models used in Appendix

To test the generality of our results, we also simulated more complicated variants of the consumer resource model (see Figure A5 and Appendix I). First, we simulated a consumer resource model with linear resource dynamics [38]:

[TABLE]

In this model resources are supplied externally at a rate rather than described by a logistic growth. This small change in resource dynamics can significantly change the ecosystem properties because it prevents resources from going extinct in the steady state. In the simulations, we set $M=100$, $\mu=1$, $\kappa=1$, $\sigma_{\kappa}=0.1$ , $m=0.1$ and $\sigma_{m}=0.01$ and each data point is averaged from 1000 independent realizations.

Second, we simulated a generalization of the MacArthur’s Consumer Resource model to a model we call the Microbial Consumer Resource Model (MicroCRM). The MicroCRM was introduced in [14] and refined in [19] to simulate microbial communities. In this model, in addition to consuming resources species can produce new resources through cross-feeding. This dramatically changes the resource dynamics through the introduction of trophic feedbacks. Unlike the original CRM and the CRM with linear resource dynamics, the MicroCRM possesses no Lyapunov function. Full details of the model are available in the appendix of [19, 15]. In particular, the dynamics we use are described in equation (17) of [19] with the leakage rate $l=0.4$. The fraction of secretion flux secreted to the same resource type is $f_{s}=0.45$, the fraction of secretion flux to ’waste’ resource is $f_{w}=0.45$ and variability in secretion fluxes among resources is $d_{0}=0.2$. We set $M=100$, $\mu=1$, $K=1$, $\sigma_{K}=0.1$ , $m=0.1$ and $\sigma_{m}=0.01$ and each data point is averaged from 1000 independent realizations.

Appendix B Sensitivity to Parameter Perturbations

We begin by defining four susceptibility matrices that measure how the steady-state resource and species abundances respond to changes in the resource supply and species death(growth) rates:

[TABLE]

where the bar $\bar{X}$ over the variable $X$ denotes the steady-state (equilibrium) solution.

For the extinct species and resources, by definition the susceptibilities are zero. For this reason, we focus only on the surviving resources and species. At steady-state, equation (1) gives:

[TABLE]

where $\mathbf{M^{*}}$ and $\mathbf{S^{*}}$ denote the sets of resources and species, respectively, that survive in the ecosystem at steady-state. Differentiating these equations yields the relations

[TABLE]

Substituting in for the partial derivatives using the susceptibility matrices defined above, we have:

[TABLE]

These two equations can be written as single matrix equation for block matrices:

[TABLE]

To solve this equation, we define a $S^{*}\times S^{*}$ matrix: $A_{ij}=\sum_{\alpha\in M^{*}}\bar{C}_{i\alpha}\bar{C}^{T}_{\alpha j}$. A straightforward calculation yields

[TABLE]

Appendix C Cavity Solution

When the designed component of the consumer preferences is the identity (i.e $\mathbf{B}=\mathds{1}$ in Eq. 9), the effect of random off-target consumption on system-scale properties can be computed analytically in the $M,S\to\infty$ limit using the cavity method [21, 22]. The cavity calculation is straightforward but tedious. For this reason, it is helpful to introduce the notation:

• •

$\frac{M^{*}}{M}=\phi_{R}$, $\left<R\right>=\frac{1}{M}\sum_{\beta}R_{\beta}$ and $q_{R}=\frac{1}{M}\sum_{\beta}R^{2}_{\beta}=\left<R^{2}\right>$ , where $M^{*}$ is the number of surviving resources.

• •

$\frac{S^{*}}{S}=\phi_{N}$, $\left<N\right>=\frac{1}{S}\sum_{j}N_{j}$ and $q_{N}=\frac{1}{S}\sum_{j}N^{2}_{j}=\left<N^{2}\right>$, where $S^{*}$ is the number of surviving species.

• •

$C_{i\alpha}\equiv{\mu\over M}+\sigma_{c}d_{i\alpha}$ assuming $\left<d_{i\alpha}\right>=0$, $\left<d_{i\alpha}d_{j\beta}\right>={\delta_{ij}\delta_{\alpha\beta}\over M}.$ with $\left<c_{i\alpha}\right>={\mu\over M}$, $\left<c_{i\alpha}c_{j\beta}\right>={\sigma_{c}^{2}\over M}\delta_{ij}\delta_{\alpha\beta}+{\mu^{2}\over M^{2}}\approx{\sigma_{c}^{2}\over M}\delta_{ij}\delta_{\alpha\beta}$.

• •

$K_{\alpha}=K+\delta K_{\alpha}$ with $\left<K_{\alpha}\right>=\frac{1}{M}\sum_{\beta}K_{\beta}=K$, $\left<\delta K_{\alpha}\delta K_{\beta}\right>=\delta_{\alpha\beta}\sigma_{K}^{2}.$

• •

$m_{i}=m+\delta m_{i}$ with $\left<m_{i}\right>=m$, $\left<\delta m_{i}\delta m_{j}\right>=\delta_{ij}\sigma_{m}^{2}.$

• •

$\gamma=\frac{M}{S}$ and for the identity matrix $\gamma=1$.

Following similar steps as in [22], we perturb the ecosystem with a new species and resource $N_{0}$ and $R_{0}$. Ignoring $\mathcal{O}(1/M)$ terms yields the following equations:

[TABLE]

Denote by $\bar{N}_{\alpha/0}$, $\bar{R}_{\alpha/0}$ and $\bar{N}_{i}$, $\bar{R}_{\alpha}$ the equilibrium values of the species and resources before and after adding the newcomers, respectively. These can be related to each other using the susceptibilities defined above:

[TABLE]

In what follows we assume Replica Symmetry. In this case, the sums in the equations above can be approximated as Gaussian random variables. For this reason, it is helpful to introduce new auxiliary random variables:

[TABLE]

where $\left<z_{N}\right>=0$, $\sigma_{z_{N}}=\sqrt{\sigma_{c}^{2}q_{R}+\sigma_{m}^{2}}$ and $\left<z_{R}\right>=0$, $\sigma_{z_{R}}=\sqrt{\sigma_{c}^{2}q_{N}+\sigma_{K}^{2}}$.

Case 1: both $R_{0}$ and $N_{0}$ are positive. Following calculations analogous to [22] and noting that $\gamma=\frac{M}{S}=1$ yields:

[TABLE]

[TABLE]

Case 2: either $R_{0}$ or $N_{0}$ is zero. We get exactly the same expression as the random ecosystem we derived in [22].

[TABLE]

Case 3: both $R_{0}$ and $N_{0}$ are zero, namely,

[TABLE]

Combining the cases above, the steady state solution is a Gaussian mixture depending on the positivity of $R_{0}$ and $N_{0}$.

[TABLE]

[TABLE]

Cavity equations for the susceptibilities can be obtained directly by differentiating these equations:

[TABLE]

C.1 With resource extinction

Two solutions are found by solving eq. (34) and eq. (35):

[TABLE]

[TABLE]

C.2 Without resource extinction

In this case, the resource never vanishes so that we can fix $\phi_{R}=1$ and solve eq. (34) and eq. (35). Two solutions are found:

[TABLE]

[TABLE]

Above two solutions are continuous at the transition point: $\chi=0$ i.e. $\phi_{N}=1$. Assume there is a small perturbation near the transition: $\phi_{N}=1-\epsilon$ and $\epsilon\ll 1$ and $\nu$ in eq. (39) can be expanded around $\epsilon$. It is easy to check the $\nu$ in eq. (39) has the same expression as eq. (38) at the first order of $\epsilon$. Therefore, only one solution exists:

[TABLE]

The comparison between cavity solutions and numerical simulations for $\chi$ and $\nu$ are given in Figure A8 and Figure 4 respectively.

C.3 Without resource extinction and species extinction

In this case, both the resource and the species never vanish so that we can fix $\phi_{R}=1$ and $\phi_{N}=1$. Solving eq. (34) and eq. (35), only one solution is found:

[TABLE]

C.4 Behavior in Three Regimes

To understand these solutions and behaviors better, it is helpful to consider three regimes: Regime A where $\chi=\phi_{R}-\phi_{N}=0$, Regime B where $\chi$ becomes nonzero and species start to extinct, and Regime C where $\sigma_{c}\gg 1$ and it becomes a random ecosystem.

In Regime B, resource extinction has a significant effect on the system’s feasibility, shown in Figure 4. With resource extinction, equation (37) shows there is a sudden change for the linear response function $\nu$ from Regime A: $\chi=0$ to Regime B $\chi\neq 0$. As $\nu\sim\frac{1}{\phi_{R}-\phi_{N}}$, even a slightly decrease of the number of surviving species will induce a huge perturbation to the ecosystem, corresponding to a phase transition between Regime A and Regime B at $\sigma_{c}^{*}\sim 0.2$.

Without resource extinction, equation (40) shows the linear response function $\nu$ is continuous from Regime A to Regime B. There is a crossover instead of a phase transition there. The peak for the crossover is a finite value and can be calculated by taking the derivative of equation (40) over $\sigma_{c}$, ignoring the correlation between $\sigma_{c}$ and $\phi_{N}$. It happens approximately at $\sigma_{c}^{*}=\sqrt{4\phi_{N}-2}\sim 1.04$, where $\phi_{N}=0.77$ can be obtained from numerical simulation. The explanation for the difference from random matrix theory are provided in the main text and also the spectrums in Figure 3 and Figure 4.

Without resource and species extinction, as shown in equation (41), $\nu$ diverges at $\sigma_{c}^{*}=1$, corresponding to $\lambda_{\rm min}$ reaching exactly zero. This result is also consistent with equation (59), predicted by random matrix theory, which ignores the effect of row or column deletions in the interaction matrix. This tells there do not exists any feasible solutions for the coexistence of M species and M resources. Therefore species must go extinct before $\sigma_{c}^{*}=1$.

In Regime C, further increasing of $\sigma_{c}$ after $\sigma_{c}>1$, the $\sigma_{c}^{4}$ term in the square root becomes dominating and the the susceptibility $\nu$ behaves like a random ecosystem quickly, which explains the dramatic drop of the species packing shown in Figure 1. It indicates the ecosystem tends to a self-organized random state.

C.5 Solutions in Regime A and C

In Regime A($\sigma_{c}\ll 1$), for eqs. (1) with resource extinction, the solutions for the steady-states become,

[TABLE]

For eqs. (2) without resource extinction, the solutions for the steady-states become,

[TABLE]

For ecosystems without resource and species extinction, the solutions for the steady-states become,

[TABLE]

For Regime C ($\sigma_{c}\gg 1$), for eqs. (1) with resource extinction, the solutions for the steady-states become,

[TABLE]

in agreement with the equations obtained in [22] for purely random interactions. For equations. (2) without resource extinction, the solutions for the steady-states become,

[TABLE]

For ecosystems without resource and species extinction, the solutions for the steady-states become,

[TABLE]

Appendix D Lotka-Volterra Model, Wishart Matrix and Marchenko-Pastur Law

In this section, we show how the generalized Lotka-Volterra model can be related to the CRM, and in particular, the how the steady states of the two models can be made to coincide. Solving for the steady-state values of the non-extinct resources by setting the bottom equation in (1) equal to zero gives:

[TABLE]

Substituting this into the top equation in (1) gives:

[TABLE]

where we have defined an interaction matrix $A_{ij}=\sum_{\alpha\in\mathbf{M^{*}}}\bar{C}_{i\alpha}\bar{C}^{T}_{\alpha j}$ and $\mathbf{M^{*}}$ is the set of surviving resources. We can use this equation to solve for the steady-state (equilibrium) abundances of non-extinct species, and arrive at the expression:

[TABLE]

where $\mathbf{S}^{*}$ is the set of surviving species. In terms of $\bar{N}_{i}$, the Lotka-Volterra equations become:

[TABLE]

with community matrix

[TABLE]

In May’s work, $J_{ij}$ is assumed to be an i.i.d. random matrix and an extension of Wigner’s arguments about Gaussian random matrices is used to compute the leading eigenvalue [8]. Since the $\bar{N}_{i}$ are not known a priori, the stability of Lotka-Volterra type dynamics are more easily studied in terms of the eigenvalues of $A_{ij}$, using the connection between the leading eigenvalues of $\mathbf{J}$ and $\mathbf{A}$ derived below.

Appendix E Relating the eigenvalues of $\mathbf{A}$ and $\mathbf{J}$

In this section, we prove that the largest eigenvalue $\lambda_{\rm max}$ of the community matrix $\mathbf{J}$ (which controls the Lyapunov stability of the fixed point) is negative if and only if the smallest eigenvalue $\lambda_{\rm min}$ of the Lotka-Volterra competition matrix $\mathbf{A}$ is positive. For this stability analysis, we remove the rows and columns corresponding to species that go extinct in the steady state, since allowing $N_{i}=0$ trivially generates zero eigenvalues. $\mathbf{J}$ and $\mathbf{A}$ will always refer to the resulting matrices of dimension $S^{*}\times S^{*}$.

We start by defining the diagonal matrix $\mathbf{\bar{N}}$, whose nonzero elements are the equilibrium population sizes $\bar{N}_{i}$. This lets us write

[TABLE]

where $\mathbf{\bar{N}}^{1/2}$ is the diagonal matrix whose entries are the square roots of the population sizes. This equation says that $\mathbf{J}$ is similar to $-\mathbf{W}\equiv-\mathbf{\bar{N}}^{1/2}\mathbf{A}\mathbf{\bar{N}}^{1/2}$, which implies that they share the same eigenvalues.

Since $\mathbf{W}$ and $\mathbf{A}$ are both symmetric matrices, their eigenvalues are all real, and the positivity of all the eigenvalues is equivalent to the positive-definiteness of the matrix.

Now we note that $\mathbf{W}$ is positive definite if and only if $\mathbf{A}$ is positive definite. For if $\mathbf{A}$ is positive definite, then $\mathbf{x}^{T}\mathbf{A}\mathbf{x}>0$ for all column vectors $\mathbf{x}\neq 0$, including the column vector $\mathbf{x}=\mathbf{\bar{N}}^{1/2}\mathbf{y}$ for any column vector $\mathbf{y}\neq 0$. But this implies that $\mathbf{y}^{T}\mathbf{\bar{N}}^{1/2}\mathbf{A}\mathbf{\bar{N}}^{1/2}\mathbf{y}>0$ for all $\mathbf{y}\neq 0$, i.e., that $\mathbf{W}$ is positive definite. Conversely, if $\mathbf{W}$ is positive definite, then $\mathbf{y}^{T}\mathbf{\bar{N}}^{1/2}\mathbf{A}\mathbf{\bar{N}}^{1/2}\mathbf{y}>0$ for all $\mathbf{y}\neq 0$, including $\mathbf{y}=\mathbf{\bar{N}}^{-1/2}\mathbf{x}$ for any $\mathbf{x}\neq 0$. But this implies that $\mathbf{x}^{T}\mathbf{A}\mathbf{x}>0$ for all $\mathbf{x}\neq 0$, i.e., that $\mathbf{A}$ is positive definite.

We conclude that the eigenvalues of $\mathbf{W}$ are all positive if and only if the eigenvalues of $\mathbf{A}$ are all positive. Therefore the largest eigenvalue of $\mathbf{J}=-\mathbf{\bar{N}}^{1/2}\mathbf{W}\mathbf{\bar{N}}^{-1/2}$ is negative if and only if the smallest eigenvalue of $\mathbf{A}$ is positive, as claimed in the main text.

An alternative but much simpler proof can be provided with the properties of the D-stable matrix[39]. A real square matrix $\mathbf{A}$ is said to be D-stable if the matrix $\mathbf{DA}$ is positive definite for every choice of a positive diagnoal matrix $\mathbf{D}$. A sufficient condition for D-stablility is that $\mathbf{A}+\mathbf{A^{T}}$ is positive definite. The Lotka-Volterra competition matrix $\mathbf{A}$ is symmetric and positive definite, $i.e.$, D-stable. $\mathbf{\bar{N}}$ is all positive. It is obvious that $\mathbf{J}=\mathbf{\bar{N}A}$ is positive definite. Further discussions about its application in ecology can be found in [40, 20]

Appendix F Correspondence between RMT and cavity solution

Our numerical simulations show that after the transition, our ecosystems are well described by purely random interactions. This suggests that we should be able to derive our cavity results using Random Matrix Theory (RMT). We now show that this is indeed the case. Our starting point are the average susceptibilities which are defined as:

[TABLE]

From the cavity calculations, we only care about $\chi^{R}_{\alpha\beta}$ and $\nu^{N}_{ij}$, because the other susceptibilities are lower order in $1/M$.

We can combine these equations with (70) and (19) to obtain

[TABLE]

We now show that the cavity solutions are consistent with results from RMT using equations (70) and(19) in Regime A and Regime C described in the main text.

F.1 Regime A: $\bar{\mathbf{C}}=\mathds{1}$

This regime happens when $\sigma_{c}\ll 1$. Substituting, $\bar{\mathbf{C}}=\mathds{1}$ into equations (70) and (19) yields

[TABLE]

This is consistent with the cavity solution equation (36) with $\sigma_{c}=0$ since in this case $S^{*}=S=M$.

F.2 Regime C: $\bar{C}_{i\alpha}\textit{ i.i.d. }\mathcal{N}(0,\sigma_{c}/\sqrt{M})$

In this regime, $\sigma_{c}\gg 1$. In this case, $A_{ij}=\sum_{\alpha\in\mathbf{S^{*}}}\bar{C}_{i\alpha}\bar{C}^{T}_{\alpha j}$ takes the form of a Wishart Matrix. We will exploit this to calculate $\chi$ and $\nu$. Notice,

[TABLE]

where $\lambda_{i}$ is the eigenvalue of $A_{ij}$. From the Marchenko-Pastur law [27], we know that the eigenvalues of a random Wishart matrix obey the Marchenko-Pastur distribution. Substituting equation (6) into the expression for $\nu$ and replacing the sum with an integral yields:

[TABLE]

The second line of equation (F.2) is obtained by transferring the integral function to a complex analytic function and applying the residue theorem. This result is the same as the cavity solution equation (37) when $\sigma_{c}\gg 1$.

F.3 Regime B using the Stieltjes transformation

In Regime B, it hard to estimate the minimum eigenvalue. We can use Stieltjes transformation of information-plus-noise-type matrices which are well studied in wireless communications[41, 31, 32], where $\mathbf{B}$ represents the information encoded in the signal and $\mathbf{C}$ is the noise in wireless communications. In this case, we have

[TABLE]

[TABLE]

Using Theorem 1.1 in Dozier and Silverstein[41], the Stieltjes transform $m(z)$ of $A_{ij}$ satisfies

[TABLE]

The asymptotic spectrum of $A_{ij}$ can be obtained by $m(z)$, the solution of equation (59) with

[TABLE]

The result is shown in Figure A1. The minimum eigenvalue reaches 0 nearly at $\sigma^{*}_{c}=1$, as predicted by the cavity solution.

We emphasize that the phase transition point, derived from eq. (59), does not change at different $\mu$. In the original paper by Marchenko and Pastur, eq. (5) requires the elements are i.i.d variables with mean 0 and variance $\sigma^{2}$. Recently, it has been shown a nonzero $\mu$ contributes only one eigenvalue $\lambda$. It is either in the domain of MP Law, $\lambda\in[a,b]$ or off the domain $\lambda>b$[42, 43] and thus it does not affect the minimum eigenvalue. We can understand it intuitively with a simple example: $\mathds{1}+\mu\mathds{J}$, where $\mathds{1}$ is the identity matrix, $\mathds{J}$ is a $n\times n$ all-ones matrix. The eigenvalues can be calculated by

[TABLE]

It shows when $n\gg 1$, it only contributes a very large eigenvalue $1+(n-1)\mu$ and the others stay at 1.

Appendix G Parameters in simulations

All simulations are done with the CVXPY package[44]. The code is available on GitHub at https://github.com/Emergent-Behaviors-in-Biology/typical-random-ecosystems.

• •

Figure 1(B), 2(B), 3(C, D): the consumer matrix $\mathbf{C}$ is sampled from the Gaussian distribution $\mathcal{N}(\frac{\mu}{M},\frac{\sigma_{c}}{\sqrt{M}})$. $S=100$, $M=100$, $\mu=0$, $K=1$, $\sigma_{K}=0.1$ , $m=0.1$, $\sigma_{m}=0.01$, and each data point is averaged from 5000 independent realizations. The model is simulated with eqs. (2).

• •

Figure 2(C): the consumer matrix $\mathbf{C}$ is sampled from the uniform distribution $\mathcal{U}(0,b)$. $S=100$, $M=100$, $\mu=0$, $K=1$, $\sigma_{K}=0.1$ , $m=0.1$, $\sigma_{m}=0.01$, and each data point is averaged from 5000 independent realizations. The model is described by eqs. (2).

• •

Figure 2(D): the consumer matrix $\mathbf{C}$ is sampled from the Bernoulli distribution $\mathit{Bernoulli}(p_{c})$. $S=100$, $M=100$, $\mu=0$, $K=1$, $\sigma_{K}=0.1$ , $m=0.1$, $\sigma_{m}=0.01$, and each data point is averaged from 5000 independent realizations. The model is described by eqs. (2).

• •

Figure 3(B), A7: the simulation is the same as Fig. 2(B). Each spectrum is drawn from 10000 independent realizations.

• •

Figure 4: the consumer matrix $\mathbf{C}$ is sampled from the Gaussian distribution $\mathcal{N}(\frac{\mu}{M},\frac{\sigma_{c}}{\sqrt{M}})$. $S=100$, $M=100$, $\mu=0$, $K=1$, $\sigma_{K}=0.1$ , $m=0.1$, $\sigma_{m}=0.01$. The model without resource extinction simulated with eqs. (2), and each data point is averaged from 5000 independent realizations.. The model with resource extinction is simulated with eqs. (1), and each data point is averaged from 4000 independent realizations. Each spectrum is drawn from 1 independent realizations for $S=500$.

• •

Figure A2 : the simulation is the same as Fig. 2(B). Each histogram is drawn from 10000 independent realizations.

• •

Figure A5: the consumer matrix $\mathbf{C}$ is sampled from the Gaussian distribution $\mathcal{N}(\frac{\mu}{M},\frac{\sigma_{c}}{\sqrt{M}})$. $S=100$, $M=100$, $\mu=1$, $K=1$, $\sigma_{K}=0.1$ , $m=0.1$, $\sigma_{m}=0.01$. For (C): $\omega=1$, $\sigma_{\omega}=0$, model details can be found in [38]; For (D), the dynamics is described in equation (17) in Supplementary Information of [19]. The noise is only applied on the consumption matrix and $D$ is kept the same at different $\sigma_{c}$. Each data point is averaged from 4000 independent realizations for (A), from 5000 independent realizations for (B) and 1000 independent realizations for (C, D).

• •

Figure A6: the simulation is the same as Figure 1(B) except $S=200$, $M=100$. For the identity case, the consumer matrix is obtained by concatenating the $M\times M$ identity plus noise matrix and a $(S-M)\times M$ gaussian random matrix. The model without resource extinction simulated with eqs. (2), and each data point is averaged from 5000 independent realizations.

• •

Figure A9, A10, and A8: the simulation is the same as Figure 4(A, B). The model without resource extinction simulated with eqs. (2), and each data point is averaged from 5000 independent realizations.. The model with resource extinction is simulated with eqs. (1), and each data point is averaged from 4000 independent realizations.

Appendix H Distinction between extinct and surviving species

In the main text, we show that the value of species packing $\frac{S^{*}}{M}$ in Fig. 1 and Fig. 2. However, in numerical simulations, even for the extinct species, the abundance is never exactly equal 0 due to numerical errors. As a result, we must choose a threshold to distinguish extinct and surviving species in order to calculate $S^{*}$. Since we are using the equivalence with convex optimization to solve the generalized consumer-resource models[45, 46], we can easily choose a reasonable threshold (e.g. $10^{-10}$ for both species since the surviving species are well separated in two peaks (see Fig. A2).

Appendix I Susceptibility matrix for linear resource dynamics

In the quasi-static limit , eqs. 10 becomes

[TABLE]

which can not be reduced to the Lotka–Volterra model.Therefore, we have to rederive the susceptibility matrix for eqs. 10.

In order to have a well defined susceptibilities, we introduce an auxiliary variable $w_{\alpha}$ and eqs. 10 become

[TABLE]

At the end we can set $w_{\alpha}$ to 1 to recover eqs. 10. Employing the results from [38] and [47], the new susceptibility matrix is

[TABLE]

where the bar $\bar{X}$ over the variable $X$ denotes the steady-state (equilibrium) solution.

For the extinct species and resources, by definition the susceptibilities are zero. For this reason, we focus only on the surviving resources and species. At steady-state, eqs. 61 gives:

[TABLE]

where $\mathbf{S^{*}}$ denote the sets of species, respectively, that survive in the ecosystem at steady-state and $\mathbf{M}$ denotes the full sets of resources as they all are nonzero for the linear resource dynamics. Differentiating these equations yields the relations

[TABLE]

Substituting in for the partial derivatives using the susceptibility matrices defined above, we have:

[TABLE]

These two equations can be written as a single matrix equation for block matrices:

[TABLE]

where $W_{\alpha}=\omega_{\alpha}\!+\!\sum_{j\in\mathbf{S^{*}}}\bar{N}_{j}C_{j\alpha}$, $G_{i\alpha}=C_{i\alpha}\bar{R}_{\alpha}$ and diag is the operator transforming a vector to a diagonal matrix.

To solve this equation, we define two $S^{*}\times S^{*}$ matrices: $A_{ij}=\sum_{\alpha\in\mathbf{M}^{*}}C_{i\alpha}C^{T}_{j\alpha}$ and $H_{ij}=(\sum_{\alpha\in\mathbf{M}^{*}}\frac{\bar{R}_{\alpha}}{W_{\alpha}}C_{i\alpha}C^{T}_{j\alpha})^{-1}$. Employing eq. (3.2) for the square off-diagonal partition, a straightforward calculation yields

[TABLE]

[TABLE]

We can see that the new susceptibilitie: $H_{ij}=(\sum_{\alpha\in\mathbf{M}^{*}}\frac{\bar{R}_{\alpha}}{W_{\alpha}}C_{i\alpha}C^{T}_{j\alpha})^{-1}$ is different with $A_{ij}=\sum_{\alpha\in\mathbf{M}^{*}}C_{i\alpha}C^{T}_{j\alpha}$ in eq. 4. Therefore, it can not behave exactly like Marchenko Pastur distribution, shown in Fig. A3 (B). However, since it is very similar to the Wishart matrix, most of our results are still preserved with eqs. (2).

Appendix J Additional figures