Statistical physics of community ecology: a cavity solution to MacArthur's consumer resource model

One of the most stunning aspects of the natural world is the diversity of species present in most ecosystems. The community structure of ecosystems are shaped through a complex interplay of the externally supplied resources available in an ecosystem, competition for these resources, as well as stochasticity [1–4]. A fundamental problem in community ecology is to understand how these processes give rise to observed pattern of species abundances. A rich theoretical framework has been developed to address this problem. Niche-based theories have emphasized the role of competition for resources [2, 5–10], while neutral theory has highlighted the role of stochastic effects [4, 11–13], and several works have investigated the interplay between stochasticity and competition [14–18].

Many of these theoretical insights have been synthesized in what is commonly referred to as contemporary niche theory. Contemporary niche theory highlights the role played by equalizing mechanisms, processes that decrease fitness differences between organisms, and stabilizing mechanisms, processes that decrease competition for resources. These basic organizational schema have been successfully applied to understand community structure in a wide range of settings [1–3].

One of the simplest and most influential mathematical models for niche theory is MacArthur's consumer resource model (MCRM) [2, 7, 8, 10]. Most analysis of MCRM—including those that inform contemporary niche theory and modern coexistence theory—have focused on small ecosystems with a few species and and few resources [2, 7, 8, 10]. However, it is unclear to what extent the theoretical insights derived from ecosystems with just a few species can be scaled up to diverse, natural ecosystems. One of the defining features of large complex systems is that they often display new 'emergent behaviors' that cannot be understood or deduced from analyzing small systems with just a few parts [19–22]. For this reason, it is essential to directly analyze large ecosystems with many resources and species and ask how they differ from the few-species ecosystems that have been analyzed previously. Recently, several works suggest that large ecosystems can exhibit unexpected behaviors such as phase transitions, emergent community-level cohesion, and the analogues of critical points [15, 23–27]. This highlights the need for new theoretical frameworks for directly analyzing large, heterogeneous ecosystems.

Perhaps the most successful and ubiquitous approaches for analyzing large systems in statistical physics is mean field theory. We emphasize that what is meant by a mean field theories in statistical physics is distinct from the way it is commonly understood in ecology [28, 29]. Unlike most usages in ecology, mean field theories in physics account for not only the means of various quantities but also fluctuations around the mean. In this paper, whenever we use the term mean field theory, we will mean it in this broader statistical physics definition rather than the narrow usage common in ecology. Mean field models have long history in statistical physics and have played a central role in the study of phase transitions and collective emergent behaviors in physical systems [28, 30]. Most mean field theories in physics focus on homogenous systems with identical components and couplings. However, more sophisticated variants such as the cavity method can be used to analyze heterogeneous 'disordered systems' [31]. Here, we develop a statistical physics inspired mean field theory, based on a generalization of the cavity method, and use it to analyze diverse ecosystems. In this paper, we will refer to this as the cavity method (CM).

Our methods are inspired by and build upon recent work showing the connection between community ecology the physics of disordered systems [15, 23, 24, 27, 32–37]. It is also closely related to the statistical mechanics of interacting socio-economic agents [38]. However, unlike these previous works our analysis explicitly incorporates resource dynamics, including resource heterogeneity and depletion. This allows us to naturally connect our results to contemporary niche theory and modern coexistence theory. One of the most striking aspects of our analysis is that we find that modification of fitness and carrying capacity due to collective phenomena is a generic feature of all diverse ecosystems [39]. In diverse ecosystems, organisms can and do significantly reshape their environments by changing resource abundances and, importantly, depleting resources. Moreover, we show that many of the central theoretical quantities in our novel CM have natural ecological interpretations that generalize many classical quantities and results of niche theory to large ecosystems and quantify the effect of collective phenomena in shaping community structure.

In this work, we will analyze one of the canonical and most influential models in community ecology: MacArthur's consumer resource model (MCRM) [7, 8]. MCRM consists of S species or consumers with abundances N_i ($i=1 \ldots S$) that can consume one of M substitutable resources with abundances $R_\alpha$ ($\alpha=1 \ldots M$). The consumer preferences of species i for resource α are encoded by a $S \times M$ matrix, $c_{i \alpha}$.

In the MCRM, the growth rate $g_i({\bf R})$ of a species depends of the concentration of all the resources. To model the growth rate, following MacArthur, we assume that a species i have some minimum maintenance cost, m_i, that they must meet. The growth rate, $g_i({\bf R})$, is proportional to amount of resources consumed, weighted by a quality factor $w_\alpha$, minus this maintenance cost

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g_i({\bf R})= \sum_\alpha c_{i \alpha} w_\alpha R_\alpha -m_i. \label{Eq:growth_rate} \nonumber \end{align} \tag{ 1 }$$

If g_i  >  0, then this is also the growth rate of species i.

The resources can be considered as 'biotic', i.e. organisms which are themselves being consumed, this is the scenario the original MCRM was designed to describe. These resources have their own internal dynamics which, following MacArthur, we assume can be modeled using logistic growth. Furthermore, when a resource is consumed, it's abundance is reduced. This ecological dynamics is captured by the following coupled, nonlinear differential equations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\rm d}N_i \over {\rm d}t} &= N_i g_i({\bf R}) \nonumber \\ {{\rm d}R_{\alpha} \over {\rm d}t}&= F_\alpha (R_\alpha) - \sum_i N_i c_{i \alpha} R_\alpha, \label{Eq:MCRM-original} \nonumber \end{align} \tag{ 2 }$$

where $F_\alpha(R_\alpha)=R_\alpha(K_\alpha -R_\alpha)$ describes the resource dynamics in the absences of consumption and $K_\alpha$ is the carrying capacity of each resource α. Such dynamics are standard when modelling biotic resources. Note also that both the species and resource abundances N_i and $R_\alpha$ must be strictly non-negative. For our analysis of the MCRM, it will be useful to define an 'effective resource capacity'

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K_\alpha^{\rm eff} ({\bf N}) =K_\alpha - \sum_i N_i c_{i \alpha} \label{Eq:Keff} \nonumber \end{align} \tag{ 3 }$$

that accounts for depletion of resources by consumers [23]. The MCRM can be rewritten in terms of $K_\alpha^{\rm eff} ({\bf N})$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\rm d}N_i \over {\rm d}t} &= N_i g_i({\bf R}), \nonumber \\ {{\rm d}R_{\alpha} \over {\rm d}t}&= R_\alpha(K_{\alpha}^{\rm eff}({\bf N}) -R_\alpha). \label{Eq:MCRM} \nonumber \end{align} \tag{ 4 }$$

A crucial property of these equations is that resources can be completely depleted from the environment. This will play an important role in what follows. Finally, we emphasize that these equations are identical those analyzed by MacArthur, Chesson, and others in deriving modern niche theory.

Previous approaches to analyzing the MCRM have largely been confined to small ecosystems with a few species and resources. Here, we consider the opposite limit of large, diverse ecosystems where both the number of species and number of resources is large, $S, M \gg 1$. In this limit, the number of parameters needed to define the ecosystem dynamics becomes extremely large. To overcome this problem, we follow a long tradition in theoretical ecology pioneered by Robert May of looking at the case where the parameters are drawn from a random distribution [40]. This allows us to ask questions about the behavior of a generic, diverse ecosystem.

We consider the case where all the consumption coefficients $c_{i\alpha}$, resource carrying capacities $K_\alpha$, and maintenance costs m_i are drawn from a random distribution. Our analytic calculations depend only on the mean and variances of the probability distributions. Denoting the expectation value of a parameters x over a distribution by $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \<x\>$, then we denote the mean and variances of our parameters by: $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \< c_{i \alpha}\> = \mu_{c} /S$, $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \< (c_{i \alpha} -\<c_{i \alpha}\>){\hspace{0pt}}^2\> = \sigma_{c}^2 /S$, $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \<m_i\>=m$, $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \<(m_i -\<m\>){\hspace{0pt}}^2\>=\sigma_m^2$, $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \<K_\alpha \>=K$, and $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \<(K_\alpha-\<K_\alpha\>){\hspace{0pt}}^2\>=\sigma_K^2$. We can also define a parameter $\gamma=M/S$ that measures the ratio of resources to species. This description is consistent with all species being generalists, that is, feeding from many resources. An examples is a predator which can feed on many different species with each prey making up a small part of its diet. In such a scenario, it is reasonable that the flux from each resources decreases with the number of resources M, through the scaling of $c_{i \alpha}$.

3.1. Invasion, ecological stability, and self-consistency

One of the cornerstones of community ecology is the idea of invasion [6, 41, 42]. In our analysis, we will ask under what circumstances a new species can invade an ecosystem. Denote the growth rate of species i when it tries to invade the ecosystem $g_i^{\rm inv}$. We will call this the invasion growth rate. Since we are interested in statistical properties, we will be primarily concerned with the mean and variances of the invasion growth rate averaged over all species i in the regional species pool: $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \< g_i^{\rm inv}\>=g$ and $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \<(g_i^{\rm inv}){\hspace{0pt}}^2\>-\<g_i^{\rm inv}\>^2=\sigma_g^2$.

The key idea that we will exploit in our analysis is the observation that as S and M get large, both the invasion growth rates, $g_i^{\rm inv}$, and the effective carrying capacities, $K_{\alpha}^{\rm eff}$ are the sum of a large number of small terms. Each individual resource makes only a small contribution (of order 1/M) to the growth of any consumer, and every consumer makes an order 1/S contribution to the effective resource capacity. Thus, from the central limit theorem, the distribution of growth rates $g_i^{\rm inv}$ and the distribution of effective resources $K_{\alpha}^{\rm eff}$ in the ecosystem can be well-approximated by a normal distribution. In the language of the cavity method of statistical physics, this corresponds to the replica symmetric solution. For future reference, denote the means and variance of the effective carrying capacity by $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \<K_{\alpha}^{\rm eff}\>=K^{\rm eff}$ and $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \<(K_{\alpha}^{\rm eff}){\hspace{0pt}}^2\>-\<K_{\alpha}^{\rm eff}\>^2= \sigma_{K^{\rm eff}}^2$ (see figure 1).

Zoom In Zoom Out Reset image size

This suggests the following intuition for thinking about our ecosystem. Each species, i, has a invasion growth rate drawn from a normal distribution. In other words, we can think of $g_i^{\rm inv} \approx g + \sigma_g z_i$, where z_i is a standard normal variable. Similarly, each resource has an effective carrying capacity that is also drawn from a normal distribution, with $K_\alpha^{\rm eff}\approx K^{\rm eff} + \sigma_K \tilde{z}_\alpha$, with $\tilde{z}_\alpha$ a standard, normal variable. In general, the means and variances ($g, K^{\rm eff}, \sigma_g^2, \sigma_{K^{\rm eff}}^2$) depend on the abundances of all other species and resources. Our statistical mechanics inspired mean field approach exploits this observation to self-consistently solve for the means and variances of the invasion growth rate and effective resource carrying capacity. In the physics literature, these is known as cavity method (CM). In general, this is a very subtle calculation but can be done using a generalized cavity equation (see below and in appendix).

In order to derive the CM self-consistency equations, we consider a system with S species and M resources and ask what happens when we add an additional species and resource to the system. We denote the abundances of the additional species and resource by N₀ and R₀ respectively. This two-step cavity where both a resource and species is removed is similar to the procedure employed to analyze the Hopfield model and compressed sensing [43, 44] and is necessary to correctly capture subtle correlations between resource and species dynamics. This approach is intimately related to classic works by MacArthur and Levins that analyzed ecological dynamics by asking if a new species could invade an ecosystem [6]. Whereas their analysis was applicable to small ecosystems with a few species, our analysis is valid for large, diverse ecosystems.

Since the number of species and resources in the original ecosystem is large ($S, M \gg 1$), the addition of the new resource and consumer represent a small perturbation of the original system. For this reason, it is useful to define two susceptibilities, χ and ν that measure the sensitivity of an ecosystem to small perturbations. The resource susceptibility, χ, measures the average change in the mean resource abundance at steady-state if we slightly increase the supply of all the externally supplied resources. Denoting the steady-state value of a quantity X by $\bar{X}$, we can mathematically define χ, as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi = {1 \over M} \sum_\alpha {\partial \bar{R}_\alpha \over \partial K_\alpha}. \nonumber \end{align} \tag{ 5 }$$

The average species-cost susceptibility, ν, measures the change in mean species abundances if we slightly decrease the minimum fitness cost (or equivalently increase the growth rate),

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nu={1 \over S} \sum_i {\partial \bar{N}_i \over \partial g_i} = -{1 \over S} \sum_i {\partial \bar{N}_i \over \partial m_i}. \nonumber \end{align} \tag{ 6 }$$

These susceptibilities characterize the sensitivity of an ecosystem to perturbations and are derivatives which can be directly measured in experiments. To estimate ν, one would alter the environment to put stress on the consumers, for instance in the context of microbes, one might force microbes to express a useless protein such as GFP as done in [45] to reduce m_i. To estimate χ, one would alter the environment to vary the carrying capacity and measure how the resource abundance fixed point changes, this could be done in many ways for one example of an experiment varying the environment to impact resource abundances see [46].

In terms of these quantities, one can derive a simple expression for the steady-state abundances of a newly added consumer and resource (see appendix):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar{N}_0 &= {\max \left[0, {g + \sigma_g {z}_0}\right] \over \gamma \sigma_{c}^2 \chi}, \nonumber \\ \bar{R}_0 &= {\max \left[0, K^{\rm eff}+ \sigma_{K^{\rm eff}}^2 \tilde{z}_0\right] \over 1-\sigma_{c}^2 \nu}, \label{EqNRmain} \nonumber \end{align} \tag{ 7 }$$

where as above z₀ and $\tilde{z}_0$ are independent, unit normal variables. These equations have a beautiful and straightforward interpretation. A new species added to the system will have an invasion growth rate $g_0^{\rm inv}=g + \sigma_g {z}_0$, which is normally distributed. If the growth rate is negative, it will not be able to invade the system and go extinct. If its growth rate is positive when introduced in the ecosystem, then it survives with an abundance proportional to its invasion growth rate. We emphasize that this proportionality constant can differ significantly from what would be expected in a single-species ecosystem and depends on all the other resources and species present in the ecosystem through the susceptibility χ and the variance of the consumption coefficients $\sigma_{c}^2$. For this reason, the invasion growth rate of a species when it invades an ecosystem is positively correlated with its abundance. Similarly, the new resource is depleted if its effective carrying capacity is negative. Otherwise the steady-state abundance of the new resource is proportional to its effective carrying capacity. These equations are similar to the arguments of MacArthur and Levins on the necessary conditions for invasibility to large ecosystems [6]. They also generalize results for species abundances derived in [27] using the Lotka–Volterra equation and the results in [24, 37, 38] which ignored resource depletion and resource fluctuations.

3.2. Comparison with numerics

Unlike small ecosystems, we cannot analytically solve for the all the resource and species abundances. However, we can take a statistical approach that allows us to calculate statistical properties of species and resource abundances at steady-state. We also restrict our analysis to uninvadable steady-states, defined as a steady-state which cannot be invaded by any species. This, both simplifies the mathematics, and allows us to more directly relate our calculations to ecology.

Using (7) is it possible to derive self-consistency equations for the fraction of species in the regional species pool that survive, $\phi_{N}$, the mean abundance of the species $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \<N\>=1/S \sum\nolimits_i N_i$, and variance and second moment of surviving species abundances, $\newcommand{\del}{\partial} \newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \<(\delta N){\hspace{0pt}}^2\>$ and $\newcommand{\del}{\partial} \newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} q_{N}= \<(\delta N){\hspace{0pt}}^2\>+ \<N\>^2= 1/S \sum\nolimits_i N_i^2$ respectively. We can also calculate the analogous equations for resources: the fraction of resources with non-zero abundance, $\phi_{R}$, the mean abundance of resources $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \<R\> = 1/M \sum\nolimits_\alpha R_\alpha$, and variance and second moment of the resource abundances, $\newcommand{\del}{\partial} \newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \<(\delta R){\hspace{0pt}}^2\>$ and $\newcommand{\del}{\partial} \newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} q_{R}= \<(\delta R){\hspace{0pt}}^2\>+ \<R\>^2 = 1/M \sum\nolimits_{\alpha} R_\alpha^2$. The equations are derived in appendix C and can be solved numerically.

The validity of our derivation is dependent on the MCRM having replica symmetry in the correspondence between the closely related cavity and replica methods (see e.g. [47]), and while we do not rigorously prove replica symmetry is not broken in this setting, we expect this assumption to be correct because there is a convex Lyapunov function (see e.g. [8]) and therefore just one fixed-point in the MCRM. Broken replica symmetry occurs when the landscape is complex, see [48]. The other assumption is a large system size $S, M \gg 1$. The theory is expected to be exact for infinite systems, but as we will now discuss, our simulations demonstrate that moderately large ecosystems are well predicted by the theory in the setting of Gaussian distributed consumption coefficients for which it was derived.

To check the accuracy of our CM, we compared our analytic predictions to numerical simulations (see figure 2). We simulated (2) for two different choices of distributions for the $c_{i \alpha}$. In the first set of simulation, the $c_{i \alpha}$ were binary random variables with $c_{i \alpha}=1$ with probability p and $c_{i \alpha}=0$ with probability 1  −  p. The probability p can be viewed as the level of generalism in the regional species pool. As $p \rightarrow 0$, all organisms in the community are specialist and consume a handful of resources. When $p \rightarrow1$, the community consists of generalists who can consume almost all resources. In the second set of simulations, we drew the consumption coefficients from a Gaussian distribution with the same mean and variance as the corresponding Bernoulli distribution with probability p.

Zoom In Zoom Out Reset image size

As shown in figure 2, our analytic results agree remarkably well with numerical simulations. The agreement between theory and numerics is nearly exact when $c_{i \alpha}$ are drawn from a Gaussian and shows qualitative agreement even when the consumption coefficients $c_{i \alpha}$ are binary random variables. This is a result of the Gaussianity assumptions used to derive the cavity equations (see appendix). The discrepancy between the binary case and Gaussian case stems from the fact that the for large S and M the $c_{i \alpha}$ are strictly positive for the binary case but generically contain some negative elements for Gaussian distributions. A negative $c_{i \alpha}$ implies that species i produces resource α at a fitness cost to itself. Thus, all simulations with Gaussian include a small fraction of public good producers that are accounted for in our theoretical calculations but are absent in the simulations with binary variables.

Despite these differences, for both choice of distributions the fraction of surviving species declines with increasing p. This is consistent with the basic idea of niche-theory that as p increases, there is increased competition resulting in greater competitive exclusion. In contrast, the mean abundances of surviving species and resources shows a non-monotonic behavior as a function of p in both numerical simulations and analytics (see appendix and figure E1 for additional simulation results).

The MacArthur consumer resource model has played a central role in the development of niche-based theories of community assembly [2, 7–10]. However, most of these analyses have focused on small ecosystems with just a few species and resources. Here, we discuss the ecological implications of our analysis for understanding community assembly in large ecosystems with many species and resources.

4.1. Relating MCRM parameters to ecology

We begin by relating the parameters of the MCRM to more ecologically meaningful quantities such as the niche overlap, fitness, zero net- growth isoclines (ZNGI), and impact vectors. In ecology, the niche overlap, ρ, measures how much two species compete for the same resources. The larger the niche overlap, the more species compete. For small ecosystems, the niche overlap is bounded between 0 and 1, with a niche overlap of zero meaning the species do not compete for resources and a niche overlap of one indicating the species have identical consumption profiles. In the context of the two species MacArthur resource model, the niche-overlap between species can be thought of as the percentage of variance explained if one performs a regression of the first consumer's consumption vector against the consumption vector of the second species [1, 7, 8]. Using this observation, we can naturally extend the idea of niche overlap to entire ecosystems by defining an ecosystem-level niche overlap ρ in terms of the mean and variances of the consumption coefficients $c_{i \alpha}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho = {\mu_{c} ^2 \over {\mu_{c}^2 +\sigma_{c}^2}}. \label{Eq:nicheoverlap} \nonumber \end{align} \tag{ 8 }$$

One useful way of thinking about ρ is that it measures the niche-overlap between two species randomly drawn from the regional species pool. It is easy to see that when $\sigma_{c}^2 \ll \mu_{c}$, all species have nearly identical consumption preferences and $\rho \rightarrow 1$. In contrast when $\sigma_{c}^2 \gg \mu_{c}$, species will have very distinct consumer preferences and $\rho \rightarrow 0$.

Another fundamental quantity in contemporary niche theory is the ecological fitness of an organism, $f_i= \sum\nolimits_{\alpha} c_{i \alpha} K_\alpha -m_i$ [1, 8]. This fitness is the initial growth rate of organism i in the absence of other species. In general, the actual growth rate of a species will differ significantly from the fitness if the resource abundances differ significantly from the resource carrying capacities $K_\alpha$. For this reason, we will refer to this as the 'naive' fitness.

We show in the appendix that it is also possible to relate our parameters directly to ZGNIs and generalized impact vectors.

4.2. Niche overlap and coexistence

One of the fundamental results of niche-based theories is that as the niche-overlap between species increases, co-existence become more and more difficult [1]. The underlying reason for this is species that have similar consumer preference are more likely to compete with each other, resulting in competitive exclusion. Thus, increasing the niche-overlap in the community should decrease the fraction of species $\phi_{N}$ that can co-exist in a community. On the other hand, stabilizing mechanisms that decrease the fitness differences between species should increase coexistence. We can parameterize the fitness differences in the community by the dimensionless quantity $\sigma_m/m$ equal to the standard deviation over the mean of the maintenance costs m_i over all species in the regional species pool. This choice of parameterization is in line with contemporary niche theory where fitness differences are defined as the difference in growth rates when species have identical consumption preferences [1]. Figure 3 shows $\phi_{N}$ as a function of the niche overlap ρ and $\sigma_m^2/m$. This choice of niche-overlap corresponds to varying the probability p for having a non-zero $c_{i \alpha}$ from 0.1 to 0.9 (see figure 2). As predicted by niche theories, increasing ρ leads to increased competition and a smaller $\phi_{N}$. In constrast decreasing $\sigma_m/m_i$ at a fixed ρ, leads to a larger fraction of species surviving. Thus, in this regard large ecosystems behave quite similarly to predictions made by analyzing smaller models.

Zoom In Zoom Out Reset image size

4.3. Resource depletion, effective fitness, and carrying capacity

One important feature of our analysis, which we are able to analyze statistically with our CM approach, is the large-scale depletion of resources. As shown in figure 2, species can significantly change the resource profile and deplete a large fraction of resources initially present in the environment. This collective behavior can change which species survive and thrive in an environment. One way to measure this effect and the reshaping of the resource profile is to measure the correlation between the naive fitness of an organism, $f_i= \sum\nolimits_{\alpha} c_{i \alpha} \max \{K_\alpha, 0\} -m_i$, and its steady-state abundance in the ecosystem $\overline{N}_i$. The fitness f_i measures the growth rate of organism i if it is introduced into an environment in the absence of other species. For this reason, we expect f_i to be highly predictive of $\overline{N}_i$ when resource abundance profiles are not significantly perturbed by consumption. On the other hand, it is possible for correlation between f_i and $\overline{N}_i$ to decrease significantly.

Figure 4 shows f_i versus $\bar{N}_i$ for numerical simulations where the $c_{i \alpha}$ drawn from a binomial distribution with p  =  0.1 and $\sigma_m/m=0.1$, as well as the case where parameters are Gaussian random variables with mean an variance matching the binomial setting. From the figure, it is clear there is a significant correlation between f_i and $\bar{N}_i$. Organisms with higher fitness disproportionately survive in the ecosystem. However, a significant number of organisms that have a high naive fitness f_i can still go extinct in the ecosystem (black points). The difference between plots (A), (B) and (C), (D) is that in the former $K_\alpha$ and m_i are kept positive by ensuring they are drawn from a gamma distribution and the consumer preferences $c_{i \alpha}$ are always positive since they are binary (1 with probability p or 0 otherwise). In (C), (D), each of these parameters is drawn from a Gaussian distribution, but with the same mean and variance as in (A), (B). This allows $K_\alpha$, m_i, and $c_{i \alpha}$ to be negative. We will be agnostic about the physical meaning of negative $c_{i \alpha}$: it could for instance be related to a common good production by a species at an energetic cost to itself or simply that species i is preyed on by resource α. A more careful breakdown of relevant behaviors which could lead to the ecological engineering phenomena we will discuss such as poisons and common good production, and this is a direction of future work, but not considered in the present manuscript.

Zoom In Zoom Out Reset image size

We find that allowing these negative consumer preference, does not dramatically alter the model predictions in terms of number/abundance or species/resources: we demonstrate a relatively good agreement between models which allow (Gaussian) and do not allow (binomial) negative consumer preferences in figure 2. However, we also find that the two cases lead to different ecological engineering phenomena when one observes individual species in the environment as in figure 4. Here, in panel D the red points correspond to resources α with a negative carrying capacity ($K_\alpha <0$) which end up in the environment due to 'help' from a set of fit species (i₁,...) because of the negative consumer preference $(c_{i_1 \alpha}, ...)$.

In figure 4(C), red points indicate species with negative fitness that can stably exist in the community (f_i  <  0, $\bar{N}_i>0$) due to the effect of negative consumer preferences, which can boost resources above their carrying capacity. Thus, we see species that cannot survive in the environment in the absence of other species can fixate (red points). Importantly, this emergent phenomenon is a collective property of the whole ecosystem and results from a complex interplay between organisms and environment and is an interesting direction for future study.

Additionally, figure 4 shows predictions from our CM for the correlation between f_i and $\overline{N}_i$. Within our replica-symmetric ansatz, these correlations are described by normal distributions whose variances and covariances can be calculated using our self-consistent equations. The contour lines represent half a standard deviation spread of our normal distribution. Our theory qualitatively captures the shape of the correlation between f_i and $\overline{N}_i$. We give explicit expression for these correlations as well as the mutual information between species abundances and naive fitness in appendix E.

Niche-based theories have played a fundamental role in shaping our understanding of community assembly and community ecology. In this work, we use ideas and methods from statistical physics to analyze a canonical model in community ecology, MacArthur's Consumer Resource Model (MCRM). Unlike previous works, our statistical physics inspired approach allows us to analyze large ecosystems with many species and resources. Our results suggest that organisms can significantly perturb their environments. The abundance of resources can be significantly altered and resource can even be completely depleted. We find that such niche-construction is a generic feature of MCRM. This suggests that in complex ecosystems, organisms actively construct their environment. To quote Levins and Lewontin, 'they are not the passive objects of external forces, but creators and modulators of these forces' [22]. These effects are even more dramatic when consumers can produce public goods at a fitness cost to themselves. In this case, species and resources that could not survive in isolation can fixate in the ecosystem.

To carry out our analysis, we developed a sophisticated theory based on the cavity method. One of the most striking things about our analysis is that many physical quantities that appear in the 'cavity equations' have natural ecological interpretations in terms of invasion growth rates and effective carrying capacities. The underlying reason for this is that the cavity methods is based on asking how ecosystems are perturbed when a new species and a new resource are introduced into the ecosystem. Conceptually, this is very similar in spirit to many classical arguments in community ecology pioneered by Levins and MacArthur that ask whether a new species can invade [2, 6]. This naturally allows us to generalize many of the results from niche-based theories to large, diverse ecosystems. However, the price we pay for using our cavity approach is that we are limited to making statistical predictions.

An important question for future investigation is to ask how our results change if we make the model more realistic. In the MCRM, all species are assumed to have a linear, Type I functional response. It will be interesting to generalize our model to non-linear functional responses. We have also neglected the effects of environmental and demographic stochasticity. Stochasticity can induce phases transitions in ecosystems from a niche-like phase where competitive effects dominate community assembly to an ecologically neutral-like phase where stochasticity is the primary determinant of community structure [15, 23]. It will be interesting to see if the techniques developed here can be generalized to this more complicated setting. Finally, we have assumed that our population can be modeled as a well-mixed community. However, spatial effects can qualitatively change the behavior of cellular populations [49, 50] and are likely to play an important role in community assembly.

We would like to thank Josh Goldford, Kirill Korolev, Seppe Kuehn, Alvaro Sanchez, Daniel Segrè, Cui Wenping for many useful discussions. PM was supported by NIH NIGMS grant 1R35GM119461, a Simons Investigator award in the Mathematical Modeling of Living Systems (MMLS), and a Scialog grant from the Simons Foundation and Research Corporation. MA was supported by the Swartz Program in Theoretical Neuroscience at Harvard.

We briefly summarize MacArthur's classical consumer resource model (MCRM). Species $i =1 \ldots S$ grows at a rate proportional to its utilization of resources, $R_\alpha$, $\alpha=1\ldots M$, in the environment. This is described by the equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {1 \over N_i} {{\rm d}N_i \over {\rm d}t} = \sum_{\alpha} c_{i \alpha} R_{\alpha} -m_i +h_i, \label{CRM_N} \nonumber \end{align} \tag{ A.1 }$$

where $w_\alpha$ is the value of one unit of resource to a species (e.g. ATPs that can be extracted); $c_{i \alpha}$ is the rate at which species i consumes resource α and converts that into a 'growth rate', m_i; m_i is the minimum amount of resources that must be consumed in order to have a positive growth rate. We have also added a small perturbation h_i to the system that will do a linear expansion in. The original MCRM corresponds to the choice h_i  =  0. We define the growth rate to be

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g_i({\bf R})= \sum_\alpha c_{i \alpha} w_\alpha R_\alpha -m_i. \nonumber \end{align} \tag{ A.2 }$$

In consumer resource model, resources satisfy their own dynamical equations:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\rm d}R_{\alpha} \over {\rm d}t}= R_\alpha (K_\alpha +b_\alpha -R_\alpha)- \sum_j c_{j \alpha} N_j R_\alpha, \label{CRM_R} \nonumber \end{align} \tag{ A.3 }$$

where the first term (with $b_\alpha =0$) describes the resource dynamics in the absence of any species and the second term models the consumption of resource by species in the environment, and $b_\alpha$ is small perturbation. The original MCRM corresponds to the choose $b_\alpha=0$. Furthermore, define the effective carrying capacity

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K_{\alpha}^{\rm eff}=K_\alpha -\sum_j c_{j \alpha} N_j. \nonumber \end{align} \tag{ A.4 }$$

We will consider the case when the consumer preferences $c_{i \alpha}$ are random variables that can be characterized by their means and variances. In particular,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \< c_{i \alpha}\>= {\mu_{c} \over S}, \nonumber \end{align} \tag{ A.5 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle \<c_{i \alpha} c_{j \beta}\>= {\sigma_{c}^2 \over S} \delta_{ij} \delta_{\alpha \beta} + {\mu_{c}^2 \over S^2} \approx {\sigma_{c}^2 \over S}. \nonumber \end{align} \tag{ A.6 }$$

To perform the cavity equations, it is useful to define several other quantities. Let us a define the fluctuating part of the consumer preferences $d_{i \alpha}$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c_{i \alpha} \equiv {\mu_{c} \over S} + \sigma_{c} d_{i \alpha}. \nonumber \end{align} \tag{ A.7 }$$

Then, we have that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \< d_{i \alpha}\>=0, \nonumber \end{align} \tag{ A.8 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle \<d_{i \alpha} d_{j \beta} \>={\delta_{ij} \delta_{\alpha \beta} \over S}. \nonumber \end{align} \tag{ A.9 }$$

We will also assume that the carrying capacities are drawn from a Gaussian distribution with

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \< K_\alpha \>=K, \nonumber \end{align} \tag{ A.10 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle \<\delta K_\alpha \delta K_\beta\>=\< (K_\alpha -K)(K_\beta-K)\>= \delta_{\alpha \beta} \sigma_K^2. \nonumber \end{align} \tag{ A.11 }$$

Finally, we assume that the minimum survival coefficients are also drawn from Gaussian distribution with

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<m_i\>=m, \nonumber \end{align} \tag{ A.12 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle \< \delta m_i \delta m_j\>= \<(m_i-m)(m_j-m)\>=\delta_{ij}\sigma_m^2. \nonumber \end{align} \tag{ A.13 }$$

For future reference, it will also be helpful to define the ratio

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \gamma={M \over S}, \nonumber \end{align} \tag{ A.14 }$$

the average resource abundance,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<R\>= {1 \over M} \sum_{\alpha} R_\alpha, \nonumber \end{align} \tag{ A.15 }$$

and the average species abundance

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<N\> = {1 \over S} \sum_j N_j. \nonumber \end{align} \tag{ A.16 }$$

With these definitions, notice that we can rewrite (A.1) as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle {1 \over N_i} {{\rm d}N_i \over {\rm d}t} = \mu \gamma \<R\> -m +\sigma \sum_{\alpha} d_{i \alpha} R_\alpha - \delta m_i, \label{CRM_N2} \nonumber \end{align} \tag{ A.17 }$$

and rewrite (A.3) as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle {1 \over R_\alpha}{{\rm d}R_{\alpha} \over {\rm d}t}= [K -\mu_{c} \<N\>] -R_\alpha -\sigma_{c} \sum_j d_{j \alpha} N_j + \delta K_\alpha. \nonumber \end{align} \tag{ A.18 }$$

We can define the mean growth rate of the population

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle g= \mu_{c} \gamma \<R\> -m \nonumber \end{align} \tag{ A.19 }$$

and the mean effective capacity of resources in the ecosystem to be

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle K^{\rm eff}=K -\mu_{c} \<N\> \nonumber \end{align} \tag{ A.20 }$$

as in the main text. In terms of these quantities, we can rewrite these equations as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle {1 \over N_i} {{\rm d}N_i \over {\rm d}t} &= g+\sigma \sum_{\alpha} d_{i \alpha} R_\alpha - \delta m_i, \nonumber \\ {1 \over R_\alpha}{{\rm d}R_{\alpha} \over {\rm d}t} &= K^{\rm eff}-R_\alpha -\sigma_{c} \sum_j d_{j \alpha} N_j + \delta K_\alpha. \nonumber \end{align} \tag{ A.21 }$$

The terms on the right hand sides of the equation above have a natural interpretation as the 'fluctuating parts' of the growth rate and effective carrying capacity. In particular, we have rewritten the growth rate for species i as the sum of the mean growth rate g and a fluctuating component $\newcommand{\del}{\partial} \delta g_i$ defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle \delta g_i =\sigma \sum_{\alpha} d_{i \alpha} R_\alpha - \delta m_i. \nonumber \end{align} \tag{ A.22 }$$

We have split the effective carrying capacity of resource α is divided into a mean $K^{\rm eff}$ and fluctuating component $\newcommand{\del}{\partial} \delta K^{\rm eff}_\alpha$ defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle \delta K_\alpha ^{\rm eff} = - \sum_j d_{j \alpha} N_j + \delta K_\alpha. \nonumber \end{align} \tag{ A.23 }$$

To derive the cavity equations, we will relate a system with S species and M resources to a new system where we add an additional resource R₀ and and additional species N₀. Thus, the cavity equations relate a ecosystem with S  +  1 and M  +  1 resources to a ecosystem with S and M resources.

Then we can write equations for this new ecosystem (to leading order in S):

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle {1 \over N_i} {{\rm d}N_i \over {\rm d}t} =g +\sigma_{c} \sum_{\alpha} d_{i \alpha} R_\alpha - \delta m_i + \sigma_{c} d_{i0} R_0, \label{CRM_N0} \nonumber \end{align} \tag{ B.1 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle {1 \over R_\alpha}{{\rm d}R_{\alpha} \over {\rm d}t}= K^{\rm eff} -R_\alpha -\sigma_{c} \sum_j d_{j \alpha} N_j + \delta K_\alpha - \sigma_{c} d_{0 \alpha}N_0. \label{CRM_R0} \nonumber \end{align} \tag{ B.2 }$$

We can also write down the corresponding equations for the new resource and species:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle {1 \over N_0} {{\rm d}N_0 \over {\rm d}t}= g +\sigma_{c} \sum_{\alpha} d_{0 \alpha} R_\alpha - \delta m_0 + \sigma_{c} d_{0 0} R_0, \label{CEN} \nonumber \end{align} \tag{ B.3 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle {1 \over R_0}{{\rm d}R_0 \over {\rm d}t}=K^{\rm eff} -R_0 -\sigma_{c} \sum_j d_{j 0} N_j + \delta K_0 - \sigma_{c} d_{0 0}N_0. \label{CER} \nonumber \end{align} \tag{ B.4 }$$

We now focus on steady-state. Let us denote the steady-state value of a quantity X by $\bar{X}$. Then, we can define some susceptibilities that are extremely useful for what follow:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_{i\beta}^{(N)} &= {\partial \bar{N}_i \over \partial K_\beta} \nonumber \\ \chi_{\alpha \beta}^{(R)} &= {\partial \bar{R}_\alpha \over \partial K_\beta} \nonumber \end{align} \tag{ B.5 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nu_{i j}^{(N)} &= {\partial \bar{N}_i \over \partial m_j} \nonumber \\ \nu_{\alpha j}^{(R)} &= {\partial \bar{R}_\alpha \over \partial m_j}. \nonumber \end{align} \tag{ B.6 }$$

Now we are in a position to perform the cavity calculation. Let us denote the steady-state value of a quantity X in the absence of the new resource and species as $\bar{X}_{/0}$. Then, since the addition of a resource and species represents a small perturbation (order 1/S), we can write:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar{N}_i= \bar{N}_{i /0} -\sigma_{c} \sum_{\beta} \chi_{i\beta}^{(N)}d_{0 \beta} \bar{N}_0 -\sigma_{c} \sum_{j} \nu_{ij}^{(N)}d_{j0}\bar{R}_0, \label{Nlinearresponse} \nonumber \end{align} \tag{ B.7 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar{R}_\alpha = \bar{R}_{\alpha /0} - \sigma_{c} \sum_{\beta} \chi_{\alpha \beta}^{(R)} d_{0 \beta} \bar{N}_0 -\sigma_{c} \sum_{j} \nu_{\alpha j}^{(R)}d_{j0}\bar{R}_0. \label{Rlinearresponse} \nonumber \end{align} \tag{ B.8 }$$

We can now plug in these expressions into the steady-state equations for N₀ and R₀. This gives:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle 0=\bar{N}_0 \left[g + \sigma d_{0 0} R_0 + \delta m_0 +\sigma_{c} \sum_{\alpha} d_{0 \alpha} R_{\alpha /0} -\sigma_{c}^2 \sum_{\alpha \beta} \chi_{\alpha \beta}^{(R)} d_{0 \alpha} d_{0 \beta} N_0 -\sigma_{c}^2 \sum_{\alpha , j} \nu_{\alpha j}^{(R)} d_{0 \alpha} d_{j 0} R_0 \right]. \nonumber \end{align} \tag{ B.9 }$$

If we now take leading order contributions to S in this expression, and take expectation value over expressions this reduces to

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle 0=\bar{N}_0 \left[g+ \sigma d_{00} R_0 -\gamma \sigma^2 {1 \over M} \sum_{\alpha} \chi_{\alpha \alpha}^{(R)} N_0 + \delta m_0 +\sigma \sum_{\alpha} d_{0 \alpha} R_{\alpha /0}\right]. \nonumber \end{align} \tag{ B.10 }$$

Notice that, to leading order in S, we can model the term $\newcommand{\del}{\partial} \delta m_0 +\sigma \sum\nolimits_{\alpha} d_{0 \alpha} R_{\alpha /0}$, which is just the invasion growth rate minus the mean growth rate $g_0^{\rm inv}-g$, as a Gaussian random field with mean 0 and variance:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_g^2= \sigma_{c}^2 \gamma {1 \over M} \sum_{\alpha} \bar{R}^2_{\alpha / 0} + \sigma_m^2 = \sigma_{c}^2 \gamma q_{R} + \sigma_m^2, \label{defsigmaN} \nonumber \end{align} \tag{ B.11 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle q_{R} ={1 \over M} \sum_{\alpha} R_{\alpha /0}^2 . \nonumber \end{align} \tag{ B.12 }$$

If we let z_N be random field with mean zero and unit variance, and define the average suspectibility

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi= {1 \over M} \sum_{\alpha} \chi_{\alpha \alpha}^{(R)} \nonumber \end{align} \tag{ B.13 }$$

then we can write the equation for N₀ as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle 0=\bar{N}_0 \left[\mu \gamma \<R\> -m -\gamma \sigma^2 \chi \bar{N}_0 + \sigma_g z_{N} \right]. \nonumber \end{align} \tag{ B.14 }$$

We can also derive a similar equation for R₀. This is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle 0= R_0 \left[K^{\rm eff} -R_0 -\sigma_{c} \sum_j d_{j 0} N_{j /0} -\sigma_{c}^2 \sum_{j \beta} \chi_{i\beta}^{(N)} d_{j 0} d_{0 \beta} \bar{N}_0 -\sigma_{c}^2 \sum_{j, k} \nu_{j k}^{(N)} d_{j 0} d_{k 0} \bar{R}_0+ \delta K_0 - \sigma_{c} d_{00}N_0 \right]. \nonumber \end{align} \tag{ B.15 }$$

Using the same logic as above, to leading order we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 0= R_0 \left[K^{\rm eff} -R_0 +\sigma_{c}^2 \nu \bar{R}_0 +\sigma_{K^{\rm eff}} z_{R} \right], \nonumber \end{align} \tag{ B.16 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nu = {1 \over S} \sum_j \nu_{j j}^{(N)}, \nonumber \end{align} \tag{ B.17 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{K^{\rm eff}}^2 = \sigma_K^2+ \sigma_{c}^2 q_{N}, \label{defsigmaR} \nonumber \end{align} \tag{ B.18 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle q_{N} ={1 \over S} \sum_{j} N_{j /0}^2 . \nonumber \end{align} \tag{ B.19 }$$

We can solve these equations and get

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle N_0 = {\max \left[0, {g + \sigma_g z_{N}}\right] \over \gamma \sigma^2 \chi} \nonumber \end{align} \tag{ B.20 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R_0 = {\max \left[0, K^{\rm eff} + \sigma_{K^{\rm eff}} z_{R} \right] \over 1-\sigma^2 \nu}. \label{EqNR} \nonumber \end{align} \tag{ B.21 }$$

Thus, the distributions for N and R are given by truncated Gaussians.

Let us now write some self-consistency equations in the replica symmetric phase. Let us define the number of non-zero species and resources as S^* and M^* respectively. Furthermore, define

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \phi_{S} = {S^* \over S}, \nonumber \end{align} \tag{ C.1 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \phi_M = {M^* \over M}. \nonumber \end{align} \tag{ C.2 }$$

Our goal is, given some parameters $\{K, \sigma_K, m, \sigma_m, \mu, \sigma, S, M\}$, to find the values for $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \{\phi_{S}, \phi_{N}, \<N\>, \<R\>, q_{R}, q_{N}, \chi, \nu\}$. Since there are eight unknowns we will need eight equations. It will also be useful to define:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \Delta_g = {\mu \gamma \<R\>-m \over \sigma_g}, \nonumber \end{align} \tag{ C.3 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \Delta_{K^{\rm eff}} = {K-\mu\<N\> \over \sigma_{K^{\rm eff}}}, \label{DefDelta} \nonumber \end{align} \tag{ C.4 }$$

and the function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle w_j(\Delta)= \int_{-\Delta}^\infty {{\rm d}z \over \sqrt{2 \pi}} {\rm e}^{-z^2 \over 2} (z+\Delta)^{\,j}. \nonumber \end{align} \tag{ C.5 }$$

First, let us write self-consistency equations for the susceptibilities. Taking derivatives with respect to m and K of (B.21) and noting that the fraction of non-zero species and non-zero resources is $\phi_{N}$ and $\phi_{R}$ respectively gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nu = -{\phi_{N} \over \gamma \sigma_{c}^2 \chi}, \nonumber \end{align} \tag{ C.6 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi= {\phi_{R} \over 1- \sigma_{c}^2 \nu}. \label{SupSCE} \nonumber \end{align} \tag{ C.7 }$$

Notice now that if we define $y= \max[0, {a \over b} + {c \over b} z]$ with z a gaussian random variable we have that:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<y^{\,j}\>= \left({b \over c}\right)^{\,j} \int_{-{b \over a}}^ \infty {\rm d}z {1 \over \sqrt{2 \pi}} {\rm e}^{-z^2 \over 2}(z+{b \over a})^{\,j}. \nonumber \end{align} \tag{ C.8 }$$

We can now use the fact that (B.21) implies that the species distribution and resource distribution is given by a truncated Gaussian to write self consistency equations for the fraction of nonzero resources and species as well as the the moments of their abundances:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \phi_{N} = w_0(\Delta_g), \label{eq:phi_N} \nonumber \end{align} \tag{ C.9 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \phi_{R} = w_0(\Delta_{K^{\rm eff}}),\label{eq:phi_R} \nonumber \end{align} \tag{ C.10 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<N\>= \left({\sigma_g \over \gamma \sigma_{c}^2 \chi}\right)w_1(\Delta_g), \nonumber \end{align} \tag{ C.11 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<R\>=\left({\sigma_{K^{\rm eff}} \over 1-\sigma_{c}^2 \nu}\right)w_1(\Delta_{K^{\rm eff}}), \label{eq:avg_R} \nonumber \end{align} \tag{ C.12 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle q_{N}=\<N^2\>=\left({\sigma_g \over \gamma \sigma_{c}^2 \chi}\right)^2 w_2 (\Delta_g), \nonumber \end{align} \tag{ C.13 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle q_{R}=\<R^2\>=\left({\sigma_{K^{\rm eff}} \over 1-\sigma_{c}^2 \nu}\right)^2 w_2(\Delta_{K^{\rm eff}}). \label{ExSCE} \nonumber \end{align} \tag{ C.14 }$$

Together (C.4), (C.7) and (C.14) define the 10 self-consistency equations we need, along with the definitions (B.11) and (B.18).

We solve these mean field equations numerically using the sum of squared differences between the left and right sides of equations (C.9)–(C.14) as an energy function which we minimize using the basinhopping optimization algorithm from the scipy.optimize. The algorithm uses random perturbations, local minimization, and an accept or reject criterion to attempt to minimize function which may be non-convex. The parameters we used were a temperature of 1, a step size of 0.5, and 5 iterations or initializations. Note that the equations can also be solved iteratively, but we found these solutions were stable for a smaller set of parameter values using this approach.

We can also easily relate our mean-field quantities to ZGNIs. Recall, that ZGNI's delineate range of resource conditions in which a species maintains a positive growth rate [9, 10]. Each species i defines a ZGNI in the resource space ${\bf R}_i^{\rm ZGNI}$ defined by the equation $g_i ({\bf R}_i^{\rm ZGNI})=0$. Geometrically, we can view ${\bf R}_i^{\rm ZGNI}$ as a hyperplane in resource space whose dot product with the consumption coefficients $c_{i \alpha}$ of species i equals m_i (see equation (1)). If $m_i \ll 1$, the ZGNI is well-approximated by the plane perpendicular to $c_{i \alpha}$. We can calculate some statistical properties of these ZGNI. Notice that the mean value of each component is just

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<R_{i \alpha}^{\rm ZGNI}\>= {m \over \gamma \mu_{c}}, \nonumber \end{align} \tag{ D.1 }$$

and the expected value of the square is just

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<(R_{i \alpha}^{\rm ZGNI})^2 \>= {\sigma_m^2 + m^2 \over \gamma \sigma_{c}^2}. \nonumber \end{align} \tag{ D.2 }$$

In Modern Niche Theory, another important quantity is the impact vector of a species i. The impact vector describes how resources are depleted by the addition of another individual. Here, we introduce the idea of generalized impact vectors that measure how the steady state concentration of a resource α changes due to the introduction of a species j Alternatively, we can consider a system without resource β and then ask how it's addition changes the species abundance of a species i. These define the generalized impact vectors (GIVs).

These are of course the leading order contribution (in S) to the cavity equations under the replica symmetric assumption, namely (B.8) and (B.7). Thus, the components of the two 'generalized' impact vectors are given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar{R}_\alpha - \bar{R}_{\alpha /j}= - \sigma_{c} \sum_{\beta} \chi_{\alpha \beta}^{(R)} d_{j\beta} \bar{N}_j, \nonumber \end{align} \tag{ D.3 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar{N_i}- \bar{N}_{i /\beta} = -\sigma_{c} \sum_{j} \nu_{ij}^{(N)}d_{j\beta}\bar{R}_\beta. \nonumber \end{align} \tag{ D.4 }$$

We want to quantify how much the naive fitness (the growth rate of an organism without other species present) is correlated to the invasion growth rate $g_i^{\rm inv}$, which in turn is closely related to the steady state abundance of each of the species in the community. From the definition of the invasion growth rate:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g_i^{\rm inv} = \sum_\alpha{c_{i\alpha}R_{\alpha / i}} - m_i, \nonumber \end{align} \tag{ E.1 }$$

and naive fitness:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_i = \sum_\alpha{c_{i\alpha}K_\alpha} - m_i, \nonumber \end{align} \tag{ E.2 }$$

we can compute the level of correlation between the two using the CM. To begin, we compute the means of each of these distributions:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<f_i\> = \sum_\alpha{\<c_{i \alpha}\> \<K_\alpha\>} - \<m_i\> = \gamma \mu_{c} K - m, \nonumber \end{align} \tag{ E.3 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<g^{\rm inv}_i\> = \sum_\alpha{\<c_{i \alpha}\> \<R_{\alpha/i}\>} - \<m_i\> = \gamma \mu_{c} \< R \> - m. \nonumber \end{align} \tag{ E.4 }$$

The final equality in the preceding equation holds in the asymptotic limit of large S and M. Note that the consumer preferences $c_{i \alpha}$ of an individual species are independent of the resources steady state levels when that species i is not included in the community as in the previous equation. We can also compute the correlation between fluctuations from the mean naive fitness and the mean invasion growth rate: i.e. if a species has a higher or lower individual fitness, how should we expect this to impact its growth rate in the community? To understand this correlation, we define $\newcommand{\del}{\partial} \newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \delta f_i = f_i - \< f_i\>$ and $\newcommand{\del}{\partial} \newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \delta g^{\rm inv}_i = g^{\rm inv}_i - \< g^{\rm inv}_i\>$ and compute:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle \<\delta f_i \delta g^{\rm inv}_i\> &= \< \left(\sum_\alpha{\left[\sigma_{c} d_{i \alpha} K_\alpha + \frac{\mu_{c}}{S} (K_\alpha -\<K_\alpha\>)\right]} + \delta m_i \right) \nonumber \\ & \cdot \left(\sum_\alpha{\left[\sigma_{c} d_{i \alpha} R_{\alpha/i} + \frac{\mu_{c}}{S} (R_{\alpha/i} -\<R_{\alpha/i}\>)\right]} + \delta m_i \right) \>. \nonumber \end{align*}$$

In the large S limit, the important terms remaining in the average above are:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle \<\delta f_i \delta g^{\rm inv}_i\> = \gamma \frac{\sigma_{c}^2}{M} \sum_\alpha{R_{\alpha/i} K_\alpha} + \sigma_m^2, \nonumber \end{align} \tag{ E.5 }$$

thus in the asymptotic limit we can write the correlation between the two forms of fitness as:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle \<\delta f_i \delta g^{\rm inv}_i\> = \gamma \sigma_{c}^2 \< R_{\alpha/i} K_\alpha \> + \sigma_m^2. \nonumber \end{align} \tag{ E.6 }$$

To compute the correlation between carrying capacity and resource level we modify (B.21) from our capacity calculation, which yields such a relationship:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle R_\alpha(K_\alpha) = {\max \left[0, K_\alpha-\mu_{c}\<N\> + \sigma_{c} \sqrt{q_{N}} z_{N}\right] \over 1-\sigma_{c}^2 \nu}. \nonumber \end{align} \tag{ E.7 }$$

Using this relationship and letting k be drawn from the same distribution as $K_\alpha$, where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle R(k) = {\max \left[0, k-\mu_{c}\<N\> + \sigma_{c} \sqrt{q_{N}} z_{N}\right] \over 1-\sigma_{c}^2 \nu}, \nonumber \end{align} \tag{ E.8 }$$

we compute

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \< k R(k) \>_{k, z_{N}}. \nonumber \end{align} \tag{ E.9 }$$

The full form of this integral is thus:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<\int{{\rm d}k k R(k) {\rm e}^{-\frac{(k - K)^2}{2 \sigma_K^2}}}\>_{z_{N}}. \nonumber \end{align} \tag{ E.10 }$$

By rewriting this as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<\int{{\rm d}k (k - K)R(k) {\rm e}^{-\frac{(k - K)^2}{2 \sigma_K^2}}} + \int{{\rm d}k K R(k) {\rm e}^{-\frac{(k - K)^2}{2 \sigma_K^2}}}\>_{z_{N}}, \nonumber \end{align} \tag{ E.11 }$$

it may be simplified via integration by parts on the first term, yielding:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \< k R(k) \>_{k, z_{N}} &= \sigma_K^2 \<R'(k)\> + K \< R(k)\> \nonumber \\ &= \frac{1}{1 - \sigma^2 \nu} \left(\sigma_K^2 w_0(\Delta_{K^{\rm eff}}) + K w_1(\Delta_{K^{\rm eff}})\right), \nonumber \end{align*}$$

where $\newcommand{\ra}{\rightarrow} \renewcommand{\>}{\rangle} \newcommand{\<}{\langle} \Delta_{K^{\rm eff}} = \frac{K - \<N\>}{\sigma_{K^{\rm eff}}}$, and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle w_j(\Delta)= \int_{-\Delta}^\infty {{\rm d}z \over \sqrt{2 \pi}} {\rm e}^{-z^2 \over 2} (z+\Delta)^{\,j}. \nonumber \end{align} \tag{ E.12 }$$

Note, we can also compute Pearson's correlation coefficient for these two fitness metrics:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle c_\rho = \frac{\<\delta f \delta g^{\rm inv}\>}{\sqrt{\<(\delta f)^2 \> \<(\delta g^{\rm inv})^2 \>}} \nonumber \end{align} \tag{ E.13 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle &= \frac{\gamma \sigma_{c}^2 \< k R(k)\>_{k, z_{N}} + \sigma_m^2}{\sqrt{(\sigma_m^2+\gamma \sigma_{c}^2 (\sigma_K^2 + K^2))(\gamma \sigma_{c}^2 \<R^2\> + \sigma_m^2)}} \nonumber \\ &= \frac{\frac{\gamma \sigma_{c}^2}{1 - \sigma_{c}^2 \nu} \left(\sigma_K^2 w_0(\Delta_{K^{\rm eff}}) + K w_1(\Delta_{K^{\rm eff}})\right) + \sigma_m^2}{\sqrt{(\sigma_m^2+\gamma \sigma_{c}^2 (\sigma_K^2 + K^2))(\gamma \sigma_{c}^2 \<R^2\> + \sigma_m^2)}}\nonumber .\label{eq:crho} \nonumber \end{align} \tag{ E.14 }$$

Using (C.10) and (C.12) we can write the preceding expression as:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \frac{\frac{\gamma \sigma_{c}^2}{1 - \sigma_{c}^2 \nu} \left(\sigma_K^2 \phi_{R} + K \left(\frac{1 - \sigma_{c}^2 \nu}{\sigma_{K^{\rm eff}}}\right)\<R\>\right) + \sigma_m^2}{\sqrt{(\sigma_m^2+\gamma \sigma_{c}^2 (\sigma_K^2 + K^2))(\gamma \sigma_{c}^2 \<R^2\> + \sigma_m^2)}}. \nonumber \end{align} \tag{ E.15 }$$

E.1. Abundance versus naive fitness

Note that given the relationship that N is a scaled version of $g^{\rm inv}$ where all negative values are truncated to zero (B.21), it follows that we can compute the correlation between $\tilde{N}_i$ and f_i, where $\tilde{N}_i = \frac{g^{\rm inv}_i}{\gamma \sigma_{c}^2 \chi}$. Where we let $N_i = \tilde{N}_i$ if $\tilde{N}_i>0$ and N_i  =  0 otherwise. To better understand that the correlation between the abundance and fitness of a species, we compute the correlation between $\tilde{N}$ and f:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle \<\delta \tilde{N} \delta f \>= \frac{1}{\gamma \sigma_{c}^2 \chi}\left(\frac{\gamma \sigma_{c}^2}{1 - \sigma_{c}^2 \nu} \left(\sigma_K^2 \phi_{R} + K \left(\frac{1 - \sigma_{c}^2 \nu}{\sigma_{K^{\rm eff}}}\right)\<R\>\right) + \sigma_m^2 \right), \nonumber \end{align} \tag{ E.16 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle \<(\delta \tilde{N})^2\>= \frac{(\gamma \sigma_{c}^2 \<R^2\> + \sigma_m^2)}{(\gamma \sigma_{c}^2 \chi)^2}, \nonumber \end{align} \tag{ E.17 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle \< (\delta f)^2\> = \sigma_m^2+\gamma \sigma_{c}^2 (\sigma_K^2 + K^2). \nonumber \end{align} \tag{ E.18 }$$

Using these covariances, along with the means:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<f\> = \gamma \mu_{c} K - m, \nonumber \end{align} \tag{ E.19 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<\tilde{N}\> = \frac{\gamma \mu_{c} \< R \> - m}{\gamma \sigma_{c}^2 \chi}, \nonumber \end{align} \tag{ E.20 }$$

we are able to generate theoretical predictions for the distribution of $f, N$. See figure 4 where the theoretical plot of f, $\tilde{N}$ is compared with values of fitness f and abundance N for all species in a network over many realizations.

E.2. Resource capacity versus resource abundance

In the same figure 4, we additionally plot a theoretical prediction overlayed with numerics of how the resource capacities $K_\alpha$ are related to the resource abundances $R_\alpha$.

If we define $\tilde{R}_\alpha = \frac{K^{\rm eff}_\alpha}{1 - \sigma_{c}^2 \nu}$, which is equal to the prediction of $R_\alpha$ when the resources abundance is positive, and use the fact

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle \delta K^{\rm eff}_\alpha = - \sum_j d_{j \alpha} N_{j} + \delta K_\alpha, \nonumber \end{align} \tag{ E.21 }$$

then we may combine these two relations to yield:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle \delta \tilde{R}_\alpha = \frac{1}{1 - \sigma_{c}^2 \nu}\left(-\sum_j d_{j \alpha} N_{j/\alpha} + \delta K_\alpha\right). \nonumber \end{align} \tag{ E.22 }$$

We can thus compute the correlations:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle \frac{1}{M}\sum_\alpha{\delta K_\alpha \delta K_\alpha} \ra \sigma_K^2, \nonumber \end{align} \tag{ E.23 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle \frac{1}{M}\sum_\alpha{\delta K_\alpha \delta \tilde{R}_\alpha} \ra \frac{\sigma_K^2}{1 - \sigma_{c}^2 \nu}, \nonumber \end{align} \tag{ E.24 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ra}{\rightarrow} \newcommand{\del}{\partial} \displaystyle \frac{1}{M}\sum_\alpha{\delta \tilde{R}_\alpha \delta \tilde{R}_\alpha} \ra \frac{\sigma_{c}^2 q_{N} + \sigma_K^2}{(1 - \sigma_{c}^2\nu)^2}. \nonumber \end{align} \tag{ E.25 }$$

Also the means are easy to compute:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<K_\alpha\> \ra K, \nonumber \end{align} \tag{ E.26 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\<}{\langle} \renewcommand{\>}{\rangle} \newcommand{\ra}{\rightarrow} \displaystyle \<\tilde{R}_\alpha\> \ra \frac{K - \mu \<N\>}{1 - \sigma_{c}^2\nu}. \nonumber \end{align} \tag{ E.27 }$$

This allows us to make theoretical predictions for how the resource abundances are correlated to resource capacities.

E.3. Mutual information between individual fitness and true growth rate

The mutual information between two Gaussian variables x and z is simply (note the means of these random variables do not contribute so we will assume them are zero mean random variables):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle I(x,z) = \frac{1}{2}\ln \left(\frac{\sigma_x^2 \sigma_z^2}{\sigma_x^2 \sigma_z^2 - \sigma_{xz}^2} \right). \nonumber \end{align} \tag{ E.28 }$$

Thus,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle I(\,f, g^{\rm inv}) = -\frac{1}{2}\ln \left(1 - c_\rho^2\right). \nonumber \end{align} \tag{ E.29 }$$

This gives us a theoretical prediction using the predicted form for the correlation coefficient (E.14).

We discussed how the theoretical curves were generated in appendix C. The numerical simulations were performed by solving the corresponding ODEs (4) and integrating numerically until time 50 000 with 1000 steps. Although it is not always needed, we improved the accuracy by additionally including a small amount of migration noise which we lowered linearly to a negligible roundoff error over the course of the integration to help ensure that a species that was favored to survive would not go extinct.

We also ran simulations in other regimes, such as the one shown in figure E1 where we consider fixing $\mu_{c} = 1$ while varying $\sigma_{c}$ to study the setting when we are less interested in comparing specialists to generalists and more interested in the effect of niche overlap and how a high overlap in the generating distribution can reduce the number of surviving species.

Zoom In Zoom Out Reset image size