Cross-feeding shapes both competition and cooperation in microbial ecosystems

Pankaj Mehta [email protected] Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, MA 02139 Robert Marsland III current address: Pontifical University of the Holy Cross, Rome, Italy Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, MA 02139

Abstract

Recent work suggests that cross-feeding – the secretion and consumption of metabolic biproducts by microbes – is essential for understanding microbial ecology. Yet how cross-feeding and competition combine to give rise to ecosystem-level properties remains poorly understood. To address this question, we analytically analyze the Microbial Consumer Resource Model (MiCRM), a prominent ecological model commonly used to study microbial communities. Our mean-field solution exploits the fact that unlike replicas, the cavity method does not require the existence of a Lyapunov function. We use our solution to derive new species-packing bounds for diverse ecosystems in the presence of cross-feeding, as well as simple expressions for species richness and the abundance of secreted resources as a function of cross-feeding (metabolic leakage) and competition. Our results show how a complex interplay between competition for resources and cooperation resulting from metabolic exchange combine to shape the properties of microbial ecosystems.

Species Packing

|

Ecology

|

Resource Dynamics

|

Microbial communities are found everywhere on the globe from hot springs to human bodies Thompson et al. (2017); Integrative et al. (2019). Large-scale surveys suggest that a these microbial communities are extremely diverse, with hundreds to thousands of distinct microbial species all coexisting in a single environment. Recent lab experiments suggest that that diverse communities can also stably coexist even in simple environments with just a handful of externally supplied resources suggesting that classical ecological theories of competition must be modified to describe microbial ecosystems Goldford et al. (2018); Enke et al. (2019); Gralka et al. (2020); Dal Bello et al. (2021); Estrela et al. (2021).

It is now believed that metabolic interactions facilitated by cross-feeding – the secretion and consumption of metabolic biproducts by microbes – play a fundamental role in shaping the structure and function of microbial ecosystems. Goldford et al. (2018); Pacheco et al. (2019); Amarnath et al. (2021); Dal Bello et al. (2021); Pacheco et al. (2021); Lechón et al. (2021).This has led to a substantial body of work trying to extend classical ecological model to incorporate cross-feeding interactions Goldford et al. (2018); Butler and O’Dwyer (2018); Marsland III et al. (2019); Goyal et al. (2018, 2021); Dal Bello et al. (2021). One prominent example of this is the Microbial Consumer Resource Model (MiCRM) which extends classical ecological consumer resource models Mac Arthur (1969); MacArthur and Levins (1967); Tilman (1982); Chesson (1990); Chase and Leibold (2009) by including cross-feeding interactions between consumers Marsland III et al. (2019); Marsland et al. (2020a); Fant et al. (2021). The MiCRM has been successfully used to explain both laboratory experiments and the natural patterns found in large-scale surveys Marsland et al. (2020a, b); Gralka et al. (2020); Fant et al. (2021).

Whereas classical consumer resource models emphasize the importance of competition in shaping ecosystems Tilman (1982); Chesson (1990); Chase and Leibold (2009), the MiCRM includes both competition for resources as well as cooperation between microbes facilitated by cross-feeding interactions. A fundamental question in microbial ecology is to understand how this interplay between competition and cooperation shapes the properties of diverse communities Lechón et al. (2021). Answering this question has proven difficult due to a lack of analytic results. The MiCRM, in contrast with many other ecological models, cannot be cast as an optimization problem Mac Arthur (1969); MacArthur and Levins (1967); Mehta et al. (2019); Marsland III et al. (2020) and for this reason is not amenable to analysis using the replica method Tikhonov and Monasson (2017); Biroli et al. (2018); Altieri et al. (2021). Here, we use the cavity method, another technique from the statistical physics of disordered systems Advani et al. (2018); Bunin (2017); Barbier et al. (2018), to derive a mean-field theory for the steady-state solutions for the MiCRM. Unlike replicas, the cavity method does not require the existence of an optimization functional making it ideally suited for analyzing the MiCRM Pearce et al. (2020). Since this calculation is quite technically involved, we limit ourselves in the main text to presenting the results and implications. A detailed derivation is provided in the appendix.

Model. We briefly summarize Microbial consumer resource model (MiCRM) (see Marsland III et al. (2019); Marsland et al. (2020a) and appendix for details). Species $i=1\ldots S$ with abundances $N_{i}$ can consume $M$ potential resources $R_{\beta}$ ( $\beta=1\ldots M$ ) that can exist in the environment. A fraction $(1-l)$ of the energy in consumed resources is used for growth while the remaining energy fraction $0\leq l\leq 1$ is released back into the environment as metabolic biproducts (see Fig. 1). This is described by the equation:

	$\displaystyle{dN_{i}\over dt}$	$\displaystyle=N_{i}\left(\Phi^{g}_{i}-m_{i}\right)$			(1)
	$\displaystyle\Phi^{g}_{i}$	$\displaystyle=\sum_{\alpha}(1-l)w_{\alpha}g_{i\alpha}c_{i\alpha}R_{\alpha}$

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 $\alpha$ , $m_{i}$ is the minimum amount of resources that must be consumed in order to have a positive growth rate, and $g_{i\alpha}$ is a factor that converts consumption into a growth rate for species $i$ when consuming resource $\alpha$ . To derive our analytic solution, we consider the case when the consumer preferences $c_{i\alpha}$ are drawn randomly with $\langle c_{i\alpha}\rangle={\mu_{c}\over M}$ and $\langle c_{i\alpha}c_{j\beta}\rangle={\sigma_{c}^{2}\over M}\delta_{ij}\delta_{\alpha\beta}+{\mu_{c}^{2}\over M^{2}}\approx{\sigma_{c}^{2}\over M}\delta_{ij}\delta_{\alpha\beta}$ . We have numerically checked that in the thermodynamic limit where $M,S\rightarrow\infty$ with $\gamma=M/S$ fixed that our analytic expressions hold even when $c_{i\alpha}$ are drawn from distributions where the $c_{i\alpha}$ are strictly positive (see Figures S1 and 2).

In the MiCRM, the resources $R_{\alpha}$ satisfy their own dynamical equations:

	$\displaystyle{dR_{\alpha}\over dt}$	$\displaystyle=$	$\displaystyle\kappa_{\alpha}-\omega_{\alpha}R_{\alpha}-\Phi_{\alpha}^{d}+\Phi_{\alpha}^{p}$		(2)
	$\displaystyle\Phi_{\alpha}^{d}$	$\displaystyle=$	$\displaystyle\sum_{j}c_{j\alpha}N_{j}R_{\alpha},\hskip 7.22743pt\Phi_{\alpha}^{p}=\sum_{j,\beta}l_{\beta}{w_{\beta}\over w_{\alpha}}D_{\alpha\beta}c_{j\beta}R_{\beta}N_{j}.$

The first two terms on the right hand side of the top equation describe the dynamics of $R_{\alpha}$ in the absence of microbes, $\Phi_{\alpha}^{d}$ describes resource consumption, and $\Phi_{\alpha}^{p}$ the production of resources due to cross feeding. The cross-feeding matrix $D_{\alpha\beta}$ encodes the fraction of energy leaked from resource $\beta$ that goes into resource $\alpha$ . Motivated by the universality of metabolism, we assume that the $D_{\alpha\beta}$ is the same for all species $i$ and each row can be described using a Dirichlet distribution. To ensure a good thermodynamic limit, we further assume that elements of $D_{\alpha\beta}$ for a fixed $\beta$ follow a Dirichlet distribution with shape parameters $\{\alpha_{j}\}$ which scale like $\alpha_{j}\sim\sqrt{M}$ . Under this assumption, to leading order in $M$ we can ignore the anti-correlations between the elements of $D_{\alpha\beta}$ when deriving our mean-field equations and instead approximate the the elements of $D_{\alpha\beta}$ as independent, normal random variables with $\langle D_{\alpha\beta}\rangle={1\over M}$ and $\langle D_{\alpha\beta}D_{\gamma\eta}\rangle=\delta_{\alpha\gamma}\delta_{\beta\eta}{\sigma^{2}_{D}\over M}$ (see appendix). We emphasize that our numerical simulations simulate the full dynamics and do not rely on this approximation.

In what follows we set $w_{\alpha}=w$ and $g_{i\alpha}=g$ for all $\alpha$ and $i$ . We also assume that the $m_{i}$ are independent normal variables with mean $m$ and variance $\sigma_{m}^{2}$ . Finally, inspired by recent experiments Goldford et al. (2018); Enke et al. (2019); Gralka et al. (2020); Dal Bello et al. (2021); Estrela et al. (2021) we focus on the case where just one of the resources $R_{E}$ is supplied externally at a rate $\kappa_{E}=\tilde{\kappa}_{E}M$ , and all the other resources result from cross-feeding (i.e $\kappa_{\alpha}=0$ for $\alpha\neq E$ ). This corresponds to the high-energy regime that supports diverse communities in Marsland III et al. (2019).

Refer to caption — Figure 1: Microbial Consumer Resource Model (MiCRM). Resources $\{R_{\alpha}\}$ are supplied externally into the environment at a rate $\kappa_{\alpha}$ and degraded at a rate $\omega$ . Microbes (also called consumers) with abundances $\{N_{i}\}$ consume resources from the environment ( $\Phi_{\alpha}^{d}$ ), using some of the energy to grow ( $\Phi_{i}^{g}$ ), while the remaining energy is secreted back into the environment as metabolic biproducts ( $\Phi_{\beta}^{p}$ ). We assume there is a single externally supplied resource ( $\alpha=E$ ) so that $\kappa_{\alpha}=0$ if $\alpha\neq E$ .

Parameters, competition, and cross-feeding. We are especially interested in understanding the interplay between competition and cooperation. In all consumer resource models, consumers/microbes that have similar preference for resources compete more than those who have very different consumer preferences since they occupy similar metabolic niches. This intuition can be generalized to diverse communities by looking at the quantity $\sigma_{c}/\mu_{c}$ which characterizes the diversity of consumer resource preferences in the regional species pool. When $\sigma_{c}\ll\mu_{c}$ , all species have the similar consumer preferences and there is lots of competition, whereas when $\sigma_{c}\gg\mu_{c}$ there is much less competition between species since species have distinct consumer preferences. For this reason, we interpret the ratio $\sigma_{c}/\mu_{c}$ as a proxy for (the inverse of) the amount of competition in the regional species pool.

The amount of cross-feeding in the community is controlled by the leakage parameter $l$ which governs the fraction of energy consumed by microbes that is released as metabolic biproducts. When $l=0$ , all the energy microbes consume is used for growth and $\Phi_{p}=0$ . In this case, our equations for the MiCRM reduce to the usual equations for a consumer resource model without cross-feeding Tikhonov and Monasson (2017); Cui et al. (2020). In contrast, if $l=1$ all the energy is released as metatabolic biproducts and species cannot grow. More generally, increasing $l$ increases the fraction of energy released as metabolic biproducts and hence serves as a proxy for the amount of cross-feeding in the community Marsland III et al. (2019).

Mean-field equations. In the appendix, we derive mean-field equations for the MiCRM using the cavity method. These calculations are extremely long and technical and here we confine ourselves to simply stating the final results. As in the consumer resource model with linear resource dynamics Cui et al. (2020), to derive our expressions we assume replica symmetry and solve the cavity-equations using a Taylor expansion that approximates the fluxes $\Phi_{i}^{g}$ , $\Phi_{\alpha}^{d}$ , and $\Phi_{\alpha}^{p}$ using their means and variances. We note that the the expectation value $\langle X\rangle$ denotes averages over both random realizations of the parameters (quenched averages) as well as species/resource abundance distributions. The end result of the cavity procedure is a set of self-consistency equations relating the fraction of surviving species $\phi_{N}$ , the mean $\langle N\rangle$ and second-moment $\langle N^{2}\rangle$ of the species abundance distribution, the mean $\langle R\rangle$ and second moment $\langle R^{2}\rangle$ of the resources secreted through cross-feeding, and the cavity susceptibilities $\chi$ and $\nu$ (see Appendix for full definitions).

These equations can be most succinctly written in terms of the average effective degradation rate of resources in the environment $\langle\omega^{\rm eff}\rangle=\omega+\mu_{c}\langle N\rangle$ , which is the sum of the dilution rate $\omega$ and the average rate at which microbes consume resources $\mu_{c}\langle N\rangle$ . In addition, we define the parameters

	$\displaystyle q$	$\displaystyle=l\left(1-{\omega\over\langle\omega_{eff}\rangle}\right)$		(3)
	$\displaystyle\epsilon_{d}^{2}$	$\displaystyle=\gamma{\sigma_{c}^{2}\over\mu_{c}^{2}}{\langle N^{2}\rangle\over\langle N\rangle^{2}},\hskip 7.22743pt\epsilon_{p}^{2}={\sigma_{D}^{2}\over D^{2}}q^{2},\hskip 7.22743pt\epsilon_{k}^{2}=\gamma^{-1}\phi_{N}+{1\over 2}\epsilon_{d}^{2}$

which characterize fluctuations in the fluxes appearing in Eq. Cross-feeding shapes both competition and cooperation in microbial ecosystems. In terms of these quantities, the species distributions satisfy the self-consistency equations

$\displaystyle\phi_{N}$	$\displaystyle=w_{0}(\Delta_{N})$	(4)
$\displaystyle\langle N\rangle$	$\displaystyle=\frac{\sigma_{N}}{(1-l)g\sigma_{c}^{2}\langle\chi R\rangle}w_{1}(\Delta_{N})$
$\displaystyle\langle N^{2}\rangle$	$\displaystyle=\left(\frac{\sigma_{N}}{(1-l)g\sigma_{c}^{2}\langle\chi R\rangle}\right)^{2}w_{2}(\Delta_{N})$
$\displaystyle w_{n}(\Delta)$	$\displaystyle=\int_{-\Delta}^{\infty}{1\over\sqrt{2\pi}}e^{-z^{2}\over 2}(z+\Delta)^{n}$
$\displaystyle\Delta_{N}$	$\displaystyle=\frac{(1-l)\mu_{c}g\langle R\rangle-m}{\sigma_{N}}$
$\displaystyle\sigma_{N}^{2}$	$\displaystyle=\sigma_{m}^{2}+(1-l)^{2}\sigma_{c}^{2}g^{2}\langle R\rangle^{2}(1+\epsilon_{d}^{2}+\epsilon_{p}^{2}),$

The full species abundace distribution can be calculated by noting that probability of observing a species with abundance $N_{0}$ in the microbial community is given by a truncated Gaussian described of the form

N_{0}=\mathrm{max}\,\left[0,{(1-l)\mu_{c}g\langle R\rangle-m+\sigma_{N}z_{0}^{g}\over(1-l)g\sigma_{c}^{2}\langle\chi R\rangle}\right],

(5)

where $z_{g}$ is a Gaussian random normal variable describing fluctuations in the growth flux $\Phi_{0}^{g}$ and the equation for the expectation $\langle\chi R\rangle$ is given below.

We also find that the moments of the resource abundances satisfy the self-consistency equations

	$\displaystyle\langle R\rangle$	$\displaystyle={\tilde{k}_{E}\over\langle\omega^{eff}\rangle}{q\over 1-q}\left(1+{\epsilon_{d}^{2}\over 1-q}-{(1+\gamma)\epsilon_{k}^{2}\over 1-q}\right)+\mathcal{O}(\epsilon^{3})$
	$\displaystyle\langle R^{2}\rangle$	$\displaystyle=\langle R\rangle^{2}(1+\epsilon_{d}^{2}+\epsilon_{p}^{2})+\mathcal{O}(\epsilon^{3}),$		(6)

, where $\mathcal{O}(\epsilon^{3})$ indicates that we have only kept terms to second order in $\epsilon_{d}$ , $\epsilon_{p}$ , and $\epsilon_{k}$ . The full resource distribution is given by

\begin{split}R_{\alpha}&=\frac{\langle R\rangle}{2\gamma^{-1}\phi_{N}+\epsilon_{d}^{2}}\times\\ &\left[\sqrt{(1+\epsilon_{d}z_{\alpha}^{d})^{2}+2(2\gamma^{-1}\phi_{N}+\epsilon_{d}^{2})(1+\epsilon_{p}z_{\alpha}^{p})}-1-\epsilon_{d}z_{\alpha}^{d}\right],\end{split}

(7)

where $z_{\alpha}^{d}$ and $z_{\alpha}^{p}$ are standard random normal variables with mean zero and variance $1$ describing the fluctuations in the depletion and production fluxes.

Finally, these equations must be supplemented by equations for expectations involving the cavity susceptibilities of the form

	$\displaystyle\tilde{\nu}\langle R\rangle$	$\displaystyle={m\over g\mu_{c}}\left(\gamma^{-1}\phi_{N}+\frac{1}{2}\epsilon_{d}^{2}\right)+\mathcal{O}(\epsilon^{3})$
	$\displaystyle\langle\chi R\rangle$	$\displaystyle=\frac{\gamma^{-1}\phi_{N}}{\tilde{\nu}}=\frac{\gamma^{-1}\phi_{N}\langle R\rangle}{\langle\omega^{\rm eff}\rangle\left(\gamma^{-1}\phi_{N}+\frac{1}{2}\epsilon_{d}^{2}\right)}.$		(8)

To check our expressions, we performed numerical simulations of the MiCRM using the Community Simulator package (see appendix and associated Github repository) Marsland et al. (2020a). In these numerical simulations a single externally supplied resource is supplied to the environment, but microbes can produce M=300 additional resources through cross-feeding (leakage rate $l=0.8$ ). The consumer preference matrix $c_{i\alpha}$ were chose from a binomial distribution. Figure S1 shows that the analytic expressions given by Eqs. 5 and 7 for the species and resource distributions (solid lines) agree well with distributions obtained numerically suggesting the cavity equations are capturing the essential ecology of the MiCRM.

Average metabolite abundance depends is controlled by leakage. One striking result of our analytic expressions is that to leading order in $\epsilon=\{\epsilon_{d},\epsilon_{p},\epsilon_{k}\}$ the average abundance of metabolic biproducts $\langle R\rangle$ depends only on the leakage rate $l$ and the effective degradation rate $\langle\omega_{eff}\rangle$ but not on the amount of competition in the regional species pool as measured by $\sigma_{c}/\mu_{c}$ (see Eq. 6). When $\langle\omega_{eff}\rangle\gg\langle\omega\rangle$ , we can use energy conservation to derive a particularly simple approximate expression for the average metabolite abundances

\langle R\rangle\approx{\langle N\rangle ml/(1-l)(\langle\omega_{eff}\rangle-\omega)}.

(9)

A comparison of this expression with numerics is shown in Fig. 2a. From this, we conclude that metabolite concentrations are affected by the regional species pool only through the average species abundance $\langle N\rangle$ and average resource degradation rate $\langle\omega_{eff}\rangle$ . This shows the resource abundances reflect the functional structure of the community rather than taxonomy (i.e. exactly what microbes are present).

Species distribution shaped by competition in regional species pool. The cavity solution also shows that the species abundance distribution is shaped primarily by competition ( $\sigma_{c}/\mu_{c}$ ) and depends much more weakly on the amount of cross-feeding as measured by the leakage rate $l$ . As shown in the appendix, we find that to leading order in $\epsilon=\{\epsilon_{d},\epsilon_{p},\epsilon_{k}\}$ , that the quantity $\Delta_{N}$ depends only on ${\sigma_{c}^{2}/\mu_{c}}$ . Since $\Delta_{N}$ is the primary determinant of species abundances, our analysis implies that both the fraction of surviving species $\phi_{N}$ as well as the ratio

\langle N^{2}\rangle/\langle N\rangle^{2}=w_{2}(\Delta_{N})/[w_{1}(\Delta_{N})]^{2}

(10)

should depend strongly on the amount of competition $\sigma_{c}/\mu_{c}$ but be largely agnostic to $l$ . Fig.2b shows numerical simulations confirming that this is indeed the case.

Species packing bound. One final consequence of our analytic solution is that it allows us to derive a species packing bound for microbial ecosystems Mac Arthur (1969); MacArthur and Levins (1967). In particular, we can ask how many species can co-exist in a microbial ecosystem with cross-feeding? We show in the appendix using a very similar argument to Cui et al. (2020) that the ratio of surviving species $S^{*}$ to metabolic biproducts $M$ must be less than half (i.e $S^{*}/M\leq 1/2$ ) in the thermodynamic limit. Fig. 3 show numerical simulations confirming this result. These results indicate that the species packing bound found for linear resource dynamics without cross-feeding in Cui et al. (2020) also holds for the MiCRM where resources are primarily produced by cross-feeding.

Discussion. In this paper, we have used the cavity method to analytically calculate self-consistent mean-field equations for the MiCRM. Our solution illustrates that an intricate interplay between competition for resources and cooperation resulting from metabolic exchange shape the properties of microbial ecosystems. We found that metabolite abundances are primarily controlled by the leakage and community-level consumption rates of resources. This finding is consistent with recent experimental and theoretical work suggesting that community-level functional structure is much more conserved in similar environments than taxonomy Louca et al. (2016a, b). Indeed, our analytic solutions suggest that the metabolite abundances are largely agnostic to which microbes are producing or consuming resources as long as the community is sufficiently diverse.

In contrast, we find that the species abundance distribution depends primarily on the amount of competition in the regional species pool and is largely agnostic to whether metabolites were supplied externally or internally generated by the community through cross-feeding. This suggests that competition is the primary force shaping taxonomic structure. However, we emphasize that cross-feeding interactions modify the environment by creating new metabolic niches and thus fundamentally shape the terms on which microbes comepete. This basic picture is given further support by our calculation showing that species packing in the MiCRM is analogous to the bound derived for ecosystems without cross-feeding where all resources are supplied externally: the number of surviving species is always less than equal to the number of metabolites in the ecosystem. In other words, species packing in diverse communities is agnostic to how resources are generated.

There are number of promising directions worth pursuing further. In our calculations, both the consumer preferences and cross-feeding matrix were unstructured and we focused on steady-states. It should be possible to extend the calculations presented here to include taxonomic and metabolic structure in the consumption and cross-feeding matrices, as well as metabolic constraints Posfai et al. (2017); Caetano et al. (2021); Pacciani-Mori et al. (2020). More ambitiously, if will be interesting to extend our calculation to derive dynamical mean-field equations Pearce et al. (2020); Roy et al. (2019). Another promising direction is to understand how space and migration modify the intuitions found here Weiner et al. (2019); Gupta et al. (2021); Dukovski et al. (2020). Finally, we would like to explore if we can exploit the mean-field equations derived here to understand eco-evolutionary dynamics in the presence of cross-feeding Good et al. (2018). Recent simulations and analytic arguments suggest that such eco-evolutionary dynamics may depend largely on the functional structure of the community, opening a potential avenue of extending these calculations to include evolutionary dynamics Good et al. (2018); Moran et al. (2021); Tikhonov et al. (2020); Fant et al. (2021).

Acknowledgments: We thank Wenping Cui, Jason Rocks, and Josh Goldford for useful discussions. This work was supported by NIH NIGMS 1R35GM119461 and a Simons Investigator in the Mathematical Modeling of Living Systems (MMLS) award to PM. We also acknowledge support from the SSC computing cluster at BU. P.M. and R.M. contributed equally to this work.

References

Thompson et al. (2017) L. R. Thompson, J. G. Sanders, D. McDonald, A. Amir, J. Ladau, K. J. Locey, R. J. Prill, A. Tripathi, S. M. Gibbons, G. Ackermann, et al., Nature 551 (2017).
Integrative et al. (2019) H. Integrative, L. M. Proctor, H. H. Creasy, J. M. Fettweis, J. Lloyd-Price, A. Mahurkar, W. Zhou, G. A. Buck, M. P. Snyder, J. F. Strauss III, et al., Nature 569, 641 (2019).
Goldford et al. (2018) J. E. Goldford, N. Lu, D. Bajić, S. Estrela, M. Tikhonov, A. Sanchez-Gorostiaga, D. Segrè, P. Mehta, and A. Sanchez, Science 361, 469 (2018).
Enke et al. (2019) T. N. Enke, M. S. Datta, J. Schwartzman, N. Cermak, D. Schmitz, J. Barrere, A. Pascual-García, and O. X. Cordero, Current Biology 29, 1528 (2019).
Gralka et al. (2020) M. Gralka, R. Szabo, R. Stocker, and O. X. Cordero, Current Biology 30, R1176 (2020).
Dal Bello et al. (2021) M. Dal Bello, H. Lee, A. Goyal, and J. Gore, Nat Ecol Evol (2021).
Estrela et al. (2021) S. Estrela, A. Sanchez-Gorostiaga, J. C. Vila, and A. Sanchez, Elife 10, e65948 (2021).
Pacheco et al. (2019) A. R. Pacheco, M. Moel, and D. Segrè, Nature communications 10, 1 (2019).
Amarnath et al. (2021) K. Amarnath, A. V. Narla, S. Pontrelli, J. Dong, T. Caglar, B. R. Taylor, J. Schwartzman, U. Sauer, O. X. Cordero, and T. Hwa, bioRxiv (2021).
Pacheco et al. (2021) A. R. Pacheco, M. L. Osborne, and D. Segrè, Nature communications 12, 1 (2021).
Lechón et al. (2021) P. Lechón, T. Clegg, J. Cook, T. P. Smith, and S. Pawar, bioRxiv (2021).
Butler and O’Dwyer (2018) S. Butler and J. P. O’Dwyer, Nature communications 9, 1 (2018).
Marsland III et al. (2019) R. Marsland III, W. Cui, J. Goldford, A. Sanchez, K. Korolev, and P. Mehta, PLoS computational biology 15, e1006793 (2019).
Goyal et al. (2018) A. Goyal, V. Dubinkina, and S. Maslov, The ISME journal 12, 2823 (2018).
Goyal et al. (2021) A. Goyal, T. Wang, V. Dubinkina, and S. Maslov, Nature communications 12, 1 (2021).
Mac Arthur (1969) R. Mac Arthur, Proceedings of the National Academy of Sciences 64, 1369 (1969).
MacArthur and Levins (1967) R. MacArthur and R. Levins, The American Naturalist 101, 377 (1967).
Tilman (1982) D. Tilman, Resource competition and community structure (Princeton university press, 1982).
Chesson (1990) P. Chesson, Theoretical Population Biology 37, 26 (1990).
Chase and Leibold (2009) J. M. Chase and M. A. Leibold, Ecological niches (University of Chicago Press, 2009).
Marsland et al. (2020a) R. Marsland, W. Cui, J. Goldford, and P. Mehta, Plos one 15, e0230430 (2020a).
Fant et al. (2021) L. Fant, I. Macocco, and J. Grilli, bioRxiv (2021).
Marsland et al. (2020b) R. Marsland, W. Cui, and P. Mehta, Scientific Reports 10, 1 (2020b).
Mehta et al. (2019) P. Mehta, W. Cui, C.-H. Wang, and R. Marsland III, Physical Review E 99, 052111 (2019).
Marsland III et al. (2020) R. Marsland III, W. Cui, and P. Mehta, The American Naturalist 196, 291 (2020).
Tikhonov and Monasson (2017) M. Tikhonov and R. Monasson, Physical review letters 118, 048103 (2017).
Biroli et al. (2018) G. Biroli, G. Bunin, and C. Cammarota, New Journal of Physics 20, 083051 (2018).
Altieri et al. (2021) A. Altieri, F. Roy, C. Cammarota, and G. Biroli, Physical Review Letters 126, 258301 (2021).
Advani et al. (2018) M. Advani, G. Bunin, and P. Mehta, Journal of Statistical Mechanics: Theory and Experiment 2018, 033406 (2018).
Bunin (2017) G. Bunin, Physical Review E 95, 042414 (2017).
Barbier et al. (2018) M. Barbier, J.-F. Arnoldi, G. Bunin, and M. Loreau, Proceedings of the National Academy of Sciences 115, 2156 (2018).
Pearce et al. (2020) M. T. Pearce, A. Agarwala, and D. S. Fisher, Proceedings of the National Academy of Sciences 117, 14572 (2020).
Cui et al. (2020) W. Cui, R. Marsland III, and P. Mehta, Physical Review Letters 125, 048101 (2020).
Louca et al. (2016a) S. Louca, S. M. Jacques, A. P. Pires, J. S. Leal, D. S. Srivastava, L. W. Parfrey, V. F. Farjalla, and M. Doebeli, Nature ecology & evolution 1, 1 (2016a).
Louca et al. (2016b) S. Louca, L. W. Parfrey, and M. Doebeli, Science 353, 1272 (2016b).
Posfai et al. (2017) A. Posfai, T. Taillefumier, and N. S. Wingreen, Physical review letters 118, 028103 (2017).
Caetano et al. (2021) R. Caetano, Y. Ispolatov, and M. Doebeli, Elife 10, e67764 (2021).
Pacciani-Mori et al. (2020) L. Pacciani-Mori, A. Giometto, S. Suweis, and A. Maritan, PLoS computational biology 16, e1007896 (2020).
Roy et al. (2019) F. Roy, G. Biroli, G. Bunin, and C. Cammarota, Journal of Physics A: Mathematical and Theoretical 52, 484001 (2019).
Weiner et al. (2019) B. G. Weiner, A. Posfai, and N. S. Wingreen, Proceedings of the National Academy of Sciences 116, 17874 (2019).
Gupta et al. (2021) D. Gupta, S. Garlaschi, S. Suweis, S. Azaele, and A. Maritan, arXiv preprint arXiv:2104.01256 (2021).
Dukovski et al. (2020) I. Dukovski, D. Bajić, J. M. Chacón, M. Quintin, J. C. Vila, S. Sulheim, A. R. Pacheco, D. B. Bernstein, W. J. Rieh, K. S. Korolev, et al., arXiv preprint arXiv:2009.01734 (2020).
Good et al. (2018) B. H. Good, S. Martis, and O. Hallatschek, Proceedings of the National Academy of Sciences 115, E10407 (2018).
Moran et al. (2021) J. Moran, D. Finlay, and M. Tikhonov, Physical Review E 103, 062402 (2021).
Tikhonov et al. (2020) M. Tikhonov, S. Kachru, and D. S. Fisher, Proceedings of the National Academy of Sciences 117, 8934 (2020).

Supplemental Appendix

Appendix A Basic Setup

We start by briefly summarize Microbial consumer resource model (MiCRM) introduced and developed in Goldford et al. (2018); Marsland III et al. (2019); Marsland et al. (2020b, a). Species, indexed $i=1\ldots S$ , grow by consuming $M$ potential resources, $R_{\beta}$ ( $\beta=1\ldots M$ ), that can be found in the the environment. In the MiCRM, the consumption of resources is never one-hundred percent efficient and a fraction $l_{\alpha}$ of the energy in the consumed resources “leak” back into the environment as metabolic biproducts.

In this MiCRM, this is described by the equation:

	$\displaystyle{dN_{i}\over dt}$	$\displaystyle=$	$\displaystyle N_{i}\left(\Phi^{g}_{i}-m_{i}\right)$
	$\displaystyle\Phi^{g}_{i}$	$\displaystyle=$	$\displaystyle\sum_{\alpha}(1-l_{\alpha})\mathrm{w}_{\alpha}g_{i\alpha}c_{i\alpha}R_{\alpha}$		(11)

where $\mathrm{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 $\alpha$ , $m_{i}$ is the minimum amount of resources that must be consumed in order to have a positive growth rate. Finally, $g_{i\alpha}$ is a factor that converts consumption into a growth rate for species $i$ when consuming resource $\alpha$ .

In the MiCRM, resources satisfy their own dynamical equations

$\displaystyle{dR_{\alpha}\over dt}$	$\displaystyle=$	$\displaystyle\kappa_{\alpha}-\omega_{\alpha}R_{\alpha}-\Phi_{\alpha}^{d}+\Phi_{\alpha}^{p}$
$\displaystyle\Phi_{\alpha}^{d}$	$\displaystyle=$	$\displaystyle\sum_{j}c_{j\alpha}N_{j}R_{\alpha}$
$\displaystyle\Phi_{\alpha}^{p}$	$\displaystyle=$	$\displaystyle\sum_{j,\beta}l_{\beta}{\mathrm{w}_{\beta}\over\mathrm{w}_{\alpha}}D_{\alpha\beta}c_{j\beta}R_{\beta}N_{j}$	(12)

where the first two terms describe resource dynamics in the absence of microbes, $\Phi_{\alpha}^{d}$ describes the depletion of resource $\alpha$ due to consumption, and the production of resources due to cross feeding from metabolic biproducts is encoded in $\Phi_{\alpha}^{p}$ . The matrix $D_{\alpha\beta}$ appearing in $\Phi_{\alpha}^{p}$ encodes the fraction of the energy released in metabolic biproduct $\alpha$ when a microbe consumes a unit of resource $\beta$ (see below for more details). We note that these equations reduce to the equations for the original Consumer Resource Model of Levins and MacArthur in the absence of leakage (i.e. $l_{\alpha}=0$ )Mac Arthur (1969); MacArthur and Levins (1967); Tilman (1982) .

To proceed, it is helpful to introduce some additional notation and make some simplifications:

•

Define the ratio of potential resources to species by

$\gamma=M/S$ (13)
•

For simplicity assume that both the growth rates and resource weights can be assumed to be constant $g_{i\alpha}=g$ and $\mathrm{w}_{\alpha}=\mathrm{w}=1$ . This is equivalent to the physical assumption that all metabolic pathways are basically equally efficient.
•

We will also assume that all resources are leaked at the same rate so that $l_{\alpha}=l$ and that the resource dilution rates are all equal with $\omega_{\alpha}=\omega$ .
•

We focus on the case where just one resource $R_{E}$ is supplied externally, with influx $\kappa_{E}$ , and all the others are internally generated through biproduct secretion. In other words, $\kappa_{\alpha}=0$ for $\alpha\neq E$ .
•

Finally, all our calculations will be performed in the thermodynamic limit with $M,S\rightarrow\infty$ and the ratio $\gamma$ fixed. For this reason, we will only keep terms that are order $1$ in both $M$ and $S$ .
•

We assume that there is no family structure in either the species consumer preferences or the cross-feeding matrix. Even this simple case has been shown to describe many real experiments on microbial communities Marsland et al. (2020b).

We will consider the case when the consumer preferences $c_{i\alpha}$ are drawn randomly with

\langle c_{i\alpha}\rangle={\mu_{c}\over M}

(14)

and

\langle c_{i\alpha}c_{j\beta}\rangle={\sigma_{c}^{2}\over M}\delta_{ij}\delta_{\alpha\beta}+{\mu_{c}^{2}\over M^{2}}\approx{\sigma_{c}^{2}\over M}\delta_{ij}\delta_{\alpha\beta},

(15)

where in writing the second line we have ignored terms that are higher order in $M$ . Since there are no correlations between the consumer preferences for different species, we have made the ecological assumption that there is no family structure in the regional species pool.

We also assume no mebaolic structure in the cross feeding matrix $D_{\alpha\beta}$ . In general, this matrix encodes the fraction of energy leaked from resource $\beta$ that goes into resource $\alpha$ . Following earlier works, we assume that each row in the matrix follows a Dirichlet distribution. To ensure a good thermodynamic limit, we assume that the shape parameters $\{\alpha_{j}\}$ of the Dirichlet distribution scale like $\alpha_{j}\sim\sqrt{M}$ and that the sum of the shape parameters $\alpha_{0}=\sum_{j=1}^{M}\alpha_{j}$ scales like $M^{3/2}$ . With these assumptions, the expectation, variance, and covariance of a variables $\{X_{i}\}$ draw from a Dirichlet distribution scale like

$\displaystyle E[X_{i}]$	$\displaystyle=$	$\displaystyle{\alpha_{i}\over\alpha_{0}}\sim{1\over M}$	(16)
$\displaystyle Var[X_{i}]$	$\displaystyle=$	$\displaystyle{\alpha_{i}(\alpha_{0}-\alpha_{i})\over\alpha_{0}^{2}(\alpha_{0}+1)}\sim{1\over M}$	(17)
$\displaystyle Cov[X_{i},X_{j}]$	$\displaystyle=$	$\displaystyle{-\alpha_{i}\alpha_{j}\over\alpha_{0}^{2}(\alpha_{0}+1)}\sim{1\over M^{2}}\approx 0.$	(18)

Thus, to leading order in $M$ we can ignore the anti-correlations between elements of the Dirichlet distribution in the thermodynamic limit. We can translate these observation to the cross-feeding matrix $D_{\alpha\beta}$ to get

\langle D_{\alpha\beta}\rangle={D\over M}={1\over M}

(19)

and

\langle D_{\alpha\beta}D_{\gamma\eta}\rangle=\delta_{\alpha\gamma}\delta_{\beta\eta}{\sigma^{2}_{D}\over M}.

(20)

Importantly, we can ignore correlations between the elements of $D_{\alpha\beta}$ . Note also that by assumption $D=1$ , since $D_{\alpha\beta}$ is defined as a matrix that partitions a given total quantity of output flux over possible metabolites, with $\sum_{\alpha}D_{\alpha\beta}=1$ .

Appendix B Energy Conservation

One important feature of the MiCRM is that it explicitly accounts for energy conservation. The total power (energy per unit time) that flows into the system is just

P_{input}=w\kappa_{E},

(21)

where $\kappa_{E}$ is the rate at which units of the external resource is added to the system and $w$ is the energy contained in each unit of resource. At steady-state, this energy is either consumed by the microbes at a rate $P_{microbe}$ or lost to the environment $P_{environment}$ due to resources disappearing at a rate $\omega_{\alpha}$ (i.e. we are considering an open system where resource $\alpha$ is lost at a rate $\omega_{\alpha}$ ). We know that by construction

P_{microbe}=\sum_{i}N_{i}(1-l_{\alpha})\mathrm{w}c_{i\alpha}R_{\alpha}=\sum_{i}g^{-1}N_{i}m_{i},

(22)

where in going to the second line we have Eq. 11, noting that at steady-state ${dN_{i}\over dt}=0$ for all microbes. We also have by definition that

P_{environment}=\sum_{\alpha\neq E}\mathrm{w}\omega_{\alpha}R_{\alpha}+\mathrm{w}\omega_{E}R_{E},

(23)

where for future convenience we have explicitly separated the externally supplied resource. Energy conservation is the statement that

P_{input}=P_{microbe}+P_{environment}

(24)

Plugging in the expression above yields the relationship

\displaystyle\kappa_{E}=\omega_{E}R_{E}+\sum_{\alpha\neq E}\omega_{\alpha}R_{\alpha}+\sum_{i}\frac{N_{i}m_{i}}{\mathrm{w}g}.

(25)

In the limit, where $P_{microbe}\gg P_{environment}$ this becomes

\kappa_{E}\approx\sum_{i}\frac{N_{i}m_{i}}{\mathrm{w}g}.

(26)

Appendix C Externally Supplied Resource

In our set-up, there is a single distinguished resource $R_{E}$ that is supplied externally. The dynamics of this resource are given by Eq. 12 with $\kappa_{E}\neq 0$ . In fact, to ensure a good thermodynamic limit we must scale $\kappa_{E}$ with the number of resources $M$ Marsland III et al. (2019). In other words, we have

\kappa_{E}=\tilde{\kappa}_{E}M.

(27)

With this choice, we see that the steady-state abundance of the externally supplied resource will also scale extensively $R_{E}\sim M$ . For this reason, we can solve for the steady-state abundance of $R_{E}$ using “naive mean-field theory”. In particular, within naive MFT, the production and depletion flux in Eq. 12 can be replaced by their mean-field averages

	$\displaystyle\langle\Phi_{E}^{d}\rangle$	$\displaystyle=$	$\displaystyle\langle\sum_{j}c_{jE}N_{j}R_{\alpha}\rangle\approx\gamma\mu_{c}\langle N\rangle R_{\alpha}$
	$\displaystyle\langle\Phi_{E}^{p}\rangle$	$\displaystyle=$	$\displaystyle\langle\sum_{j,\beta}lD_{E\beta}c_{j\beta}R_{\beta}N_{j}\rangle\approx\gamma^{-1}lD\mu_{c}\langle N\rangle\langle R\rangle,$		(28)

where we have defined the averages

	$\displaystyle\langle N\rangle={1\over S}\sum_{j=1}^{S}N_{j}$
	$\displaystyle\langle R\rangle={1\over M}\sum_{\alpha=1}^{M}R_{\alpha}.$		(29)

Plugging these equations into Eq. 12 and solving for the steady-state gives

R_{E}={\kappa_{E}+D\mu_{c}\langle N\rangle\langle R\rangle\over\omega+\gamma^{-1}\mu_{c}\langle N\rangle}\approx M{\tilde{\kappa}_{E}\over\langle\omega^{\mathrm{eff}}\rangle},

(30)

where in going to the second line we have dropped sub-leading order term in $M$ and defined the average effective degradation rate

\langle\omega^{\mathrm{eff}}\rangle=\omega+\gamma^{-1}\mu_{c}\langle N\rangle.

(31)

Thus, as expected the steady-state abundance of the externally supplied resource $R_{E}$ is extensive in $M$ .

Appendix D Means, variances, and covariances for fluxes of self-generated resources

In the thermodynamic limit, we treat the production flux $\Phi_{\alpha}^{p}$ and depletion flux $\Phi_{\alpha}^{d}$ in Eq. 12 for self-generated resources as random variables. Under our assumption of replica-symmetry, such random variables are fully characterized by their means, variances, and covariances. To calculate these quantities, we make use of the cavity method. In particular, we will consider a system in the absence of resource $\alpha$ . We will denote the steady-state abundances reached by species $N_{j}$ and resources $R_{\beta}$ in the absence of resource $\alpha$ by $N_{j/\alpha}$ and $R_{\beta/\alpha}$ . Since through out this paper we will only be considering steady-states, our notation makes no distinction between steady-state quantities and dynamical quantities. We will then relate the means and variances of the fluxes $\Phi_{\alpha}^{p}$ and $\Phi_{\alpha}^{d}$ to $N_{j/\alpha}$ and $R_{\beta/\alpha}$ using perturbation theory. The logic behind this somewhat involved cavity method is that unlike naive-MFT this method properly accounts for all correlations in the problem.

D.1 Depletion flux

Let us first write

c_{i\alpha}={\mu_{c}\over M}+{\sigma_{c}\over\sqrt{M}}x_{i\alpha}

(32)

where the $x_{i\alpha}$ are uncorrelated, standard normal variables. We first find

\Phi_{\alpha}^{d}=R_{\alpha}\gamma^{-1}\mu_{c}\langle N\rangle+R_{\alpha}{\sigma_{c}\over\sqrt{M}}\sum_{j}x_{j\alpha}N_{j/\alpha}

(33)

where through an abuse of notation we have defined the expectation value of $\langle N\rangle$ as

\langle N\rangle={1\over S}\sum_{j}N_{j/\alpha}

(34)

and $N_{j/\alpha}$ is the population of species $j$ in the equilibrium state with resource $\alpha$ absent. This gives us that the mean of the depletion flux is just

\langle\Phi_{\alpha}^{d}\rangle=\gamma^{-1}\mu_{c}\langle N\rangle R_{\alpha}

(35)

We can define

\delta\Phi_{\alpha}^{d}=\Phi_{\alpha}^{d}-\langle\Phi_{\alpha}^{d}\rangle=R_{\alpha}{\sigma_{c}\over\sqrt{M}}\sum_{j}x_{j\alpha}N_{j/\alpha}

(36)

We can also calculate the variance of this flux

	$\displaystyle\sigma_{d,\alpha}^{2}\equiv\tilde{\sigma}_{d\alpha}^{2}R_{\alpha}^{2}\equiv\langle(\delta\Phi_{\alpha}^{d})^{2}\rangle$	$\displaystyle=$	$\displaystyle R_{\alpha}^{2}{\sigma_{c}^{2}\over M}\sum_{j,k}\langle x_{j\alpha}x_{k\alpha}\rangle N_{j/\alpha}N_{k/\alpha}$		(37)
		$\displaystyle=$	$\displaystyle\gamma^{-1}\sigma_{c}^{2}\langle N^{2}\rangle R_{\alpha}^{2}$		(37)

where we have defined

\langle N^{2}\rangle={1\over S}\sum_{j}N_{j/\alpha}^{2}

(38)

and have defined our shorthand symbol $\tilde{\sigma}_{d}^{2}$ for this variance in a way that keeps the $R_{\alpha}^{2}$ visible, since it will be important later on.

This allows us to approximate the production flux as a Gaussian random variable with

\Phi_{\alpha}^{d}=\langle\Phi_{\alpha}^{d}\rangle+\tilde{\sigma}_{d,\alpha}R_{\alpha}z_{\alpha}^{d},

(39)

where $z_{\alpha}^{d}$ is a unit normal random variable and $\langle\Phi_{\alpha}^{d}\rangle$ and $\tilde{\sigma}_{d,\alpha}$ are defined in Eq. 35 and Eq. 37 respectively.

D.2 Production flux

To calculate the mean and variance of the production flux, we once again start by making the dependence of $D_{\alpha\beta}$ of the $M$ explicit by writing

D_{\alpha\beta}={D\over M}+{\sigma_{D}\over\sqrt{M}}y_{\alpha\beta},

(40)

where the $y_{\alpha\beta}$ are uncorrelated, standard normal variables. This lack of correlation follows from the discussion of the Dirichlet distribution above. In terms of these variables we have

	$\displaystyle\Phi_{\alpha}^{p}$	$\displaystyle=$	$\displaystyle l\sum_{j,\beta\neq E,\alpha}\left({D\over M}+{\sigma_{D}\over\sqrt{M}}y_{\alpha\beta}\right)\left({\mu_{c}\over M}+{\sigma_{c}\over\sqrt{M}}x_{j\beta}\right)R_{\beta/\alpha}N_{j/\alpha}$		(41)
			$\displaystyle+\gamma^{-1}lD\mu_{c}\langle N\rangle{R_{E}\over M}$		(41)

where once again we have treated the external resource in naive-MFT since it is extensive in $M$ .

We have to take care in computing the mean flux, because there is an important correlations among $N_{j}$ , $x_{j\beta}$ and $R_{\beta}$ :

\displaystyle\langle\Phi_{\alpha}^{p}\rangle

\displaystyle=\gamma^{-1}lD\left(\mu_{c}\langle N\rangle(\langle R\rangle+{R_{E}\over M})+\sqrt{M}\sigma_{c}\langle NxR\rangle\right)

(42)

where in the last line we have used that fact this second term scales like $\sqrt{M}^{-1/2}$ (see below for self-consistency proof) and we have defined

\langle R\rangle={1\over M}\sum_{\beta}R_{\beta/\alpha}

(43)

and

\langle NxR\rangle={1\over MS}\sum_{\beta j}N_{j/\alpha}x_{j\beta}R_{\beta/\alpha}.

(44)

In writing this we have used the fact that $y_{\alpha\beta}$ is not correlated with $R_{\beta/\alpha}$ and $N_{j\alpha}$ and hence all averages involving this quantity can be set to zero.

In addition to the mean value, to apply the cavity method we will also have to think about the fluctuations around the mean:

	$\displaystyle\delta\Phi_{\alpha}^{p}$	$\displaystyle\equiv$	$\displaystyle\Phi_{\alpha}^{p}-\langle\Phi_{\alpha}^{p}\rangle$		(45)
		$\displaystyle=$	$\displaystyle\frac{l\sigma_{D}}{\sqrt{M}}\sum_{j,\beta\neq\alpha}y_{\alpha\beta}\left({\mu_{c}\over M}+{\sigma_{c}\over\sqrt{M}}x_{j\beta}\right)R_{\beta/\alpha}N_{j/\alpha}-\gamma^{-1}lD\sqrt{M}\sigma_{c}\langle NxR\rangle$		(45)

A straightforward but slightly tedious calculation shows that

\sigma_{p,\alpha}^{2}\equiv\langle(\delta\Phi_{\alpha}^{p})^{2}\rangle=\gamma^{-2}l^{2}\sigma_{D}^{2}\left[\mu_{c}^{2}\langle R^{2}\rangle\langle N\rangle^{2}+\sigma_{c}^{2}M(\langle NxR^{2}xN\rangle+\langle NxR\rangle^{2})+2\mu_{c}\sigma_{c}\langle N\rangle\sqrt{M}\langle xNR^{2}\rangle\right].

(46)

where again in the approximation we have ignored terms that scale as $\sqrt{M}^{-1/2}$ and defined

	$\displaystyle\langle NxR^{2}xN\rangle$	$\displaystyle=$	$\displaystyle{1\over MS^{2}}\sum_{\beta jk}N_{j/\alpha}x_{j\beta}R_{\beta/\alpha}^{2}x_{k\beta}N_{k/\alpha}$		(47)
	$\displaystyle\langle NxR^{2}\rangle$	$\displaystyle=$	$\displaystyle{1\over MS}\sum_{\beta j}N_{j/\alpha}x_{j\beta}R_{\beta/\alpha}^{2}.$		(48)

This allows us to approximate the production flux as a Gaussian random variable with

\Phi_{\alpha}^{p}=\langle\Phi_{\alpha}^{p}\rangle+\sigma_{p,\alpha}z_{\alpha}^{p},

(49)

where $z_{\alpha}^{p}$ is a unit normal random variable and $\langle\Phi_{\alpha}^{p}\rangle$ and $\sigma_{p,\alpha}$ are defined in Eq. 42 and Eq. 46 respectively.

D.3 Covariance between production and depletion flux

From the above calculation, it is easy to show that to leading order in $M$ that the covariance between the production and depletion flux is zero. To see this, we use Eq. 45 and insert this into the definition of the covariance $\langle\delta\Phi_{\alpha}^{d}\delta\Phi_{\alpha}^{p}\rangle$ between the production and depletion flux. Since $\delta\Phi_{\alpha}^{p}$ has a single factor of the standard normal variable $y_{\alpha\beta}$ multiplying all terms, the only way to get a nonzero average is for $\delta\Phi_{\alpha}^{d}$ to contain factors of $y_{\alpha\beta}$ that could generate a $y_{\alpha\beta}^{2}$ . But the depletion flux does not have any dependence on $y_{\alpha\beta}$ , and so the correlation vanishes.

Appendix E Growth flux

In the previous section, we focused on the production and depletion fluxes for the resources. Here, we focus on the growth flux. Because we are using the cavity method, we consider a system where a species $i$ has been removed from the ecosystem. The steady-state abundances reached by the ecosystem in the absence of species $i$ are denoted by $N_{j/i}$ and $R_{\alpha/i}$ . We then ask about the growth flux that the species $i$ would have when it invades such an ecosystem:

\Phi_{i}^{g}=\sum_{\alpha}(1-l)gc_{i\alpha}R_{\alpha/i}+(1-l)g\mu_{c}{R_{E}\over M},

(50)

where once again we ignore fluctuations in the consumption of the externally supplied resource. As before, we write $c_{i\alpha}=\mu_{c}/M+{\sigma_{c}\over\sqrt{M}}x_{i\alpha}$ to get

	$\displaystyle\langle\Phi_{i}^{g}\rangle$	$\displaystyle=(1-l)\mu_{c}g\left(\langle R\rangle+{R_{E}\over M}\right)$
		$\displaystyle\approx(1-l)\mu_{c}g\left(\langle R\rangle+{R_{E}\over M}\right).$		(51)

Defining the fluctuating part of the growth flux

\delta\Phi_{i}^{g}\equiv\Phi_{i}^{g}-\langle\Phi_{i}^{g}\rangle=\sum_{\alpha}{\sigma_{c}\over\sqrt{M}}x_{i\alpha}R_{\alpha},

(52)

allows us to calculate the variance of the growth flux yielding:

\sigma_{g,i}^{2}\equiv\langle(\delta\Phi_{i}^{g})^{2}\rangle=(1-l)^{2}g^{2}\sigma_{c}^{2}\langle R^{2}\rangle.

(53)

With these expressions, within the replica-symmetry ansatz the growth flux can be though of as a normal random variable of the form

\Phi_{i}^{g}=(1-l)\mu_{c}g\langle R\rangle+\sigma_{g,i}z_{i}^{g}

(54)

Appendix F Cavity Corrections

The essence of the cavity calculations is to relate an ecosystem with $(S+1,M+1)$ species and resources to an ecosystem with $(S,M)$ system. We will use the language of adding an additional species and resource which by convention we will denote by $N_{0}$ and $R_{0}$ . The assumption is that when $S,M\gg 1$ , the addition of an additional species or resource is a small perturbation and we can relate the steady-state abundances in the presence and absence of the extra resource and species using perturbation theory.

F.1 Setting-up the perturbation theory

In the presence of the new resource $R_{0}$ and species $N_{0}$ , Eq. 12 is modified to as:

0={dR_{\alpha}\over dt}=\kappa_{\alpha}+\delta k_{\alpha}-\omega_{\alpha}R_{\alpha}-\Phi_{\alpha}^{d}+\Phi_{\alpha}^{p},

(55)

with

\delta k_{\alpha}=-c_{0\alpha}N_{0}R_{0}+\sum_{j}lD_{\alpha 0}c_{j0}R_{0}N_{j}+\sum_{\beta}D_{\alpha\beta}c_{0\beta}N_{0}R_{\beta}+D_{\alpha 0}R_{0}N_{0}

(56)

Similarly, Eq. 11 is modified to

{dN_{i}\over dt}=N_{i}\left(\Phi^{g}_{i}-m_{i}-\delta m_{i}\right)

(57)

with

\delta m_{i}=-(1-l)gc_{i0}R_{0}.

(58)

We can use perturbation theory to relate the species and resource abundances an ecosystem without the new resource and species (denoted by $N_{i/0}$ and $R_{\alpha/0}$ respectively) to the abundances in the presence of $N_{0}$ and $R_{0}$ (denoted by $N_{i}$ and $R_{\alpha}$ respectively) . To do so, we define the susceptibilities

\nu_{i}=-{\partial N_{i}\over\partial m_{i}}.

(59)

and

\chi_{\alpha}={\partial R_{\alpha}\over\partial\kappa_{\alpha}}.

(60)

We note that this latter susceptibility differs from the resource susceptibility in Cui et al. (2020) where we took the derivative with respect to $\omega_{\alpha}$ rather than $\kappa_{\alpha}$ .

Combining the expressions above yields

	$\displaystyle R_{\alpha}$	$\displaystyle=$	$\displaystyle R_{\alpha/0}+\delta R_{\alpha}=R_{\alpha/0}+\chi_{\alpha}\delta\kappa_{\alpha}$		(61)
		$\displaystyle=$	$\displaystyle R_{\alpha/0}+\chi_{\alpha}\left(-c_{0\alpha}N_{0}R_{\alpha}+l\sum_{j}D_{\alpha 0}c_{j0}N_{j}R_{0}+l\sum_{\beta}D_{\alpha\beta}c_{0\beta}R_{\beta}N_{0}\right)$		(61)

and

	$\displaystyle N_{i}$	$\displaystyle=$	$\displaystyle N_{i/0}+\delta N_{i}=N_{i/0}-\nu_{i}\delta m_{i}$		(62)
		$\displaystyle=$	$\displaystyle N_{i/0}+\nu_{i}(1-l)gc_{i0}R_{0}.$		(62)

In writing these expressions, we have ignored the “off-diagonal” susceptibilities whose contributions we have assumed scale lower order in $M$ .

F.2 Computing the TAP corrections

To proceed, we need to calculate the TAP corrections to the three fluxes $\Phi_{\alpha}^{p}$ , $\Phi_{\alpha}^{d}$ , and $\Phi_{i}^{g}$ due to the addition of $N_{0}$ and $R_{0}$ . Let us denote these cavity corrections to the fluxes by $\Delta\Phi$ . Notice through an abuse of notation we will drop the $/0$ notation so that $R_{\alpha/0}$ is represented by $R_{\alpha}$ . For the growth flux this amounts to calculating

$\displaystyle\Delta\Phi_{0}^{g}$	$\displaystyle=$	$\displaystyle\sum_{\alpha}(1-l_{\alpha})g_{0\alpha}c_{0\alpha}\delta R_{\alpha}$	(63)
	$\displaystyle=$	$\displaystyle\sum_{\alpha}(1-l)g\chi_{\alpha}\left(-\langle c_{0\alpha}c_{0\alpha}\rangle N_{0}R_{\alpha}+l\sum_{j}D_{\alpha 0}\langle c_{0\alpha}c_{j0}\rangle N_{j}R_{0}+l\sum_{\beta}D_{\alpha\beta}\langle c_{0\alpha}c_{0\beta}\rangle R_{\beta}N_{0}\right)$
	$\displaystyle=$	$\displaystyle(1-l)g\left[-{1\over M}\sum_{\alpha}\chi_{\alpha}\sigma_{c}^{2}N_{0}R_{\alpha}+{l\over M}\chi_{0}D_{00}\sigma_{c}^{2}N_{0}R_{0}+{1\over M}\sum_{\alpha}\sigma_{c}^{2}\chi_{\alpha}D_{\alpha\alpha}R_{\alpha}N_{0}\right]$
	$\displaystyle=$	$\displaystyle-(1-l)g\sigma_{c}^{2}N_{0}\langle\chi R\rangle$

where

\displaystyle\langle\chi R\rangle={1\over M}\sum_{\alpha\neq E}\chi_{\alpha}R_{\alpha}

(64)

and we have dropped the last two terms from the third line because they are higher order in $1/M$ . We also need

$\displaystyle\Delta\Phi_{0}^{d}$	$\displaystyle=$	$\displaystyle\sum_{j}c_{j0}\delta N_{j}R_{0}$	(65)
	$\displaystyle=$	$\displaystyle R_{0}^{2}\sum_{j}c_{j0}^{2}\nu_{j}(1-l)g_{j0}$
	$\displaystyle=$	$\displaystyle\gamma^{-1}R_{0}^{2}(1-l)g\sigma_{c}^{2}\nu$
	$\displaystyle\equiv$	$\displaystyle R_{0}^{2}\tilde{\nu}$

where we have introduced the notation

	$\displaystyle\nu$	$\displaystyle=$	$\displaystyle\frac{1}{S}\sum_{j}\nu_{j}$		(66)
	$\displaystyle\tilde{\nu}$	$\displaystyle=$	$\displaystyle\gamma^{-1}(1-l)g\sigma_{c}^{2}\nu.$		(67)

Finally,

\Delta\Phi_{0}^{p}=l\sum_{j,\beta}D_{0\beta}c_{j\beta}(R_{\beta}\delta N_{j}+\delta R_{\beta}N_{j})

(68)

One can show that to leading order $M$ that these terms disappear in the thermodynamic limit so that one has

\Delta\Phi_{0}^{p}\approx 0.

(69)

Appendix G Self-consistency equation

G.1 Equations for species

Let us start by defining an equation for species 0 at steady-state. We separate out the actual steady-state growth flux $\Phi_{0}^{g}$ for species 0 into three parts: the average growth flux $\langle\Phi^{g}\rangle$ over all species in the regional pool, a Gaussian random variable with variance $\sigma_{N}^{2}=Var[\Phi_{i}^{g}]$ to account for the variability in growth rates in the regional pool, and finally the TAP correction $\Delta\Phi_{0}^{g}$ that accounts for the specific feedback from species 0 on the rest of the community:

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle N_{0}(\langle\Phi_{0}^{g}\rangle+\sigma_{N}z_{0}^{g}+\Delta\Phi_{0}^{g}-m_{0})$
	$\displaystyle 0$	$\displaystyle=$	$\displaystyle N_{0}\left[(1-l)\mu_{c}g\langle R\rangle+\sigma_{N}z_{0}^{g}-(1-l)g\sigma_{c}^{2}N_{0}\langle\chi R\rangle-m_{0}\right]$		(70)

Requiring that the steady-states cannot be invaded (i.e are ecologically stable) translates into the statement that we choose $N_{0}$ to be

N_{0}=\mathrm{max}\,\left[0,{(1-l)\mu_{c}g\langle R\rangle-m+\sigma_{N}z_{0}^{g}\over(1-l)g\sigma_{c}^{2}\langle\chi R\rangle}\right],

(71)

with

\sigma_{N}^{2}=\sigma_{m}^{2}+\sigma_{g,0}^{2}=\sigma_{m}^{2}+(1-l)^{2}\sigma_{c}^{2}g^{2}\langle R^{2}\rangle,

(72)

where in writing this equation we have used the fact that $m_{0}$ is a normal random variable with variance $\delta m^{2}$ that is uncorrelated with $\Phi_{0}^{g}$ and Eq. 53.

Assuming replica symmetry, we can now write down the self-consistency equations for the moments of $N_{i}$ , using the special functions $w_{n}$ defined as

w_{n}(\Delta)=\int_{-\Delta}^{\infty}{1\over\sqrt{2\pi}}e^{-z^{2}\over 2}(z+\Delta)^{n}.

(73)

Using the usual arguments for replica-symmetric cavity solutions and matching the moments of Eq 71, one arrives at the expressions for the fraction of surviving species $\phi_{N}$ and the first and second moments of the species abundance distributions $\langle N\rangle$ and $\langle N^{2}\rangle$ Advani et al. (2018); Bunin (2017); Barbier et al. (2018):

$\displaystyle\phi_{N}$	$\displaystyle=w_{0}(\Delta_{N})$	(74)
$\displaystyle\langle N\rangle$	$\displaystyle=\frac{\sigma_{N}}{(1-l)g\sigma_{c}^{2}\langle\chi R\rangle}w_{1}(\Delta_{N})$	(75)
$\displaystyle\langle N^{2}\rangle$	$\displaystyle=\left(\frac{\sigma_{N}}{(1-l)g\sigma_{c}^{2}\langle\chi R\rangle}\right)^{2}w_{2}(\Delta_{N})$	(76)

where

\displaystyle\Delta_{N}=\frac{(1-l)\mu_{c}g\langle R\rangle-m}{\sigma_{N}}..

(77)

To calculate the average susceptibility we note that

\displaystyle\nu=-\left\langle\frac{\partial N_{0}}{\partial m_{0}}\right\rangle=\frac{\phi_{N}}{(1-l)g\sigma_{c}^{2}\langle\chi R\rangle},

(78)

and, from Eq 67 that:

\displaystyle\tilde{\nu}=\frac{\gamma^{-1}\phi_{N}}{\langle\chi R\rangle}.

(79)

Eq. 71 also allows us to compute the expectation values that occur in Eq. 42 and Eq. 46. To do so, we note that from

\displaystyle\sigma_{N}z_{i}^{g}=\delta m_{i}+(1-l)\frac{\sigma_{c}}{\sqrt{M}}\sum_{\beta}g_{i\beta}x_{i\beta}R_{\beta/i}

(80)

Using this expression and the definition of $x_{i\beta}$ , we obtain:

$\displaystyle\langle NxR\rangle$	$\displaystyle=\phi_{N}{\langle R^{2}\rangle\over\sqrt{M}\sigma_{c}\langle\chi R\rangle}$	(81)
$\displaystyle\langle NxR^{2}\rangle$	$\displaystyle=\phi_{N}{\langle R^{3}\rangle\over\sqrt{M}\sigma_{c}\langle\chi R\rangle}$	(82)
$\displaystyle\langle NxR^{2}xN\rangle$	$\displaystyle=\phi_{N}^{2}{\gamma\langle R^{2}\rangle^{2}+\langle R^{4}\rangle\over M\sigma_{c}^{2}\langle\chi R\rangle^{2}}$	(83)

The first three quantities have the correct scaling for the production flux and its variance to remain finite as $M\to\infty$ . In calculating these expressions, we have ignored correlation between $x_{i\alpha}$ and $R_{\alpha}$ which we have assumed are higher order in $1/M$ . We check that this assumption is reasonable numerically.

With these expressions, we can rewrite Eqs. 42 and 46 as

\langle\Phi_{\alpha}^{p}\rangle=\gamma^{-1}lD\left[\mu_{c}\langle N\rangle(\langle R\rangle+{R_{E}\over M})+\phi_{N}{\langle R^{2}\rangle\over\langle\chi R\rangle}\right]

(84)

and

\sigma_{p,\alpha}^{2}=\gamma^{-2}l^{2}\sigma_{D}^{2}\left[\mu_{c}^{2}\langle R^{2}\rangle\langle N\rangle^{2}+\phi_{N}^{2}{(\gamma+1)\langle R^{2}\rangle^{2}+\langle R^{4}\rangle\over\langle\chi R\rangle^{2}}+2\mu_{c}\langle N\rangle\phi_{N}{\langle R^{3}\rangle\over\langle\chi R\rangle}\right].

(85)

G.2 Equations for resources

Let us now think about the steady-state equation for resource 0. We will separate the fluxes out into the same three components – the mean over the regional pool, a random term with variance equal to the variance over the regional pool, and a feedback term. Since we already showed that the covariance between the production and depletion fluxes vanishes to leading order in $1/M$ , we can introduce independent random terms for each of these fluxes. We recall that $\kappa_{0}=0$ since we are considering a scenario with just a single externally supplied resource, but we will keep the term present for the sake of defining the susceptibilities. This gives:

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle\kappa_{0}-\omega_{0}R_{0}-\langle\Phi_{0}^{d}\rangle+\langle\Phi_{0}^{p}\rangle-\tilde{\sigma}_{d}R_{0}z_{0}^{d}+\sigma_{p,ex}z_{0}^{p}-\Delta\Phi_{0}^{d}+\Delta\Phi_{0}^{p}$
	$\displaystyle 0$	$\displaystyle\equiv$	$\displaystyle\kappa_{0}^{\mathrm{eff}}-\omega^{\mathrm{eff}}_{0}R_{0}-\tilde{\nu}R_{0}^{2},$		(86)

with

	$\displaystyle\kappa_{0}^{\mathrm{eff}}$	$\displaystyle=\kappa_{0}+\langle\Phi_{0}^{p}\rangle+\sigma_{p}z_{0}^{p}=\gamma^{-1}lD\left(\mu_{c}\langle N\rangle(\langle R\rangle+{R_{E}\over M})+\sqrt{M}\sigma_{c}\langle NxR\rangle\right)+\sigma_{p}z_{0}^{p}$		(87)
	$\displaystyle\omega^{\mathrm{eff}}_{0}$	$\displaystyle=\omega+\langle\Phi_{0}^{d}\rangle+\tilde{\sigma}_{d}z_{0}^{d}=\omega+\gamma^{-1}\mu_{c}\langle N\rangle+\tilde{\sigma}_{d}z_{0}^{d},$		(88)

where in first line we have used $\kappa_{0}=0$ and Eq. 42 and in the second equation we have used Eq. 35. This is identical to the case without cross-feeding considered in Cui et al. (2020) with $\kappa\rightarrow\kappa^{\mathrm{eff}}$ , $\omega\to\omega^{\rm eff}$

Solving for $R_{0}$ gives:

\displaystyle R_{0}

\displaystyle=\frac{-\omega^{\rm eff}_{0}+\sqrt{(\omega^{\rm eff}_{0})^{2}+4\kappa^{\rm eff}_{0}\tilde{\nu}}}{2\tilde{\nu}}.

(89)

This expression allows us to compute self-consistency equations for the first two moments of the distribution of $R_{\alpha}$ , which require performing some expansions in small parameters to solve, as in the case without crossfeeding. We can also compute

\displaystyle\chi_{0}

\displaystyle=\frac{1}{\sqrt{(\omega^{\rm eff}_{0})^{2}+4\kappa^{\rm eff}_{0}\tilde{\nu}}}

(90)

and find

\displaystyle\langle\chi R\rangle

\displaystyle=\frac{1}{2\tilde{\nu}}\left\langle 1-\frac{\omega^{\rm eff}_{\alpha}}{\sqrt{(\omega^{\rm eff}_{0})^{2}+4\kappa^{\rm eff}_{0}\tilde{\nu}}}\right\rangle.

(91)

G.3 Solving the equations for resources

Unfortunately, due to the non-linear nature of Eq. 89 it is not possible to solve for the moments of $R_{0}$ without making additional approximations. In what follows, we will perform an expansion assuming the fluctuations in the production flux $\Phi_{p}^{0}$ (Eqs. 42 and 46)and depletion flux $\Phi_{d}^{0}$ (Eqs. 35 and 37) are small. The first thing we need to do in order to solve these equations is to perform a Taylor expansion of the square root. To do so, we define three “dimensionless” parameters:

$\displaystyle\epsilon_{p}^{2}$	$\displaystyle=\frac{\langle(\delta\Phi_{0}^{p})^{2}\rangle}{\langle\Phi_{0}^{p}\rangle^{2}}=\frac{\sigma_{p}^{2}}{\langle\kappa^{\rm eff}\rangle^{2}}$	(92)
$\displaystyle\epsilon_{d}^{2}$	$\displaystyle=\frac{\langle(\delta\Phi_{0}^{d})^{2}\rangle}{\langle\Phi_{0}^{d}\rangle^{2}}\times{\langle R\rangle^{2}\over\langle R^{2}\rangle}=\frac{\tilde{\sigma_{d}}^{2}}{\langle\omega^{\rm eff}\rangle^{2}}$	(93)
$\displaystyle\epsilon_{\kappa}^{1}$	$\displaystyle=\frac{\langle\kappa^{\rm eff}\rangle\tilde{\nu}}{\langle\omega^{\rm eff}\rangle^{2}}$	(94)

From these definitions, it is clear that the first two parameters measure fluctuations in the production and depletion flux respectively. The parameter $\epsilon_{k}$ is a little harder to interpret. We will see below that the leading order contribution to $\epsilon_{\kappa}^{2}$ is the just the species-packing fraction and is always less than $1/2$ (see derivation of species packing bound below).

In terms of these parameters, we have:

\displaystyle\sqrt{(\omega_{\alpha}^{\rm eff})^{2}+4\kappa_{\alpha}^{\rm eff}\tilde{\nu}}=\langle\omega^{\rm eff}\rangle\sqrt{1+2\epsilon_{d}z_{\alpha}^{d}+\epsilon_{d}^{2}(z_{\alpha}^{d})^{2}+4\epsilon_{\kappa}^{2}(1+\epsilon_{p}z_{\alpha}^{p})}.

(95)

We will use the formula:

\displaystyle\sqrt{1+x}=1+\frac{1}{2}x-\frac{1}{8}x^{2}+\frac{1}{16}x^{3}-\frac{5}{128}x^{4}+\mathcal{O}(x^{5}).

(96)

Keeping terms up to order $\epsilon^{4}$ (and noting that several terms cancel which greatly simplifies the final expression), we have

	$\displaystyle\sqrt{(\omega_{\alpha}^{\rm eff})^{2}+4\kappa_{\alpha}^{\rm eff}\tilde{\nu}}$	$\displaystyle=1+\epsilon_{d}z_{\alpha}^{d}+2\epsilon_{\kappa}^{2}[1+\epsilon_{p}z_{\alpha}^{p}-\epsilon_{d}z_{\alpha}^{d}+\epsilon_{d}^{2}(z_{\alpha}^{d})^{2}$		(97)
		$\displaystyle-\epsilon_{d}\epsilon_{p}z_{\alpha}^{p}z_{\alpha}^{d}-\epsilon_{\kappa}^{2}]+\mathcal{O}(\epsilon^{5}).$		(98)

Using the expansion in eq. (98) in the expression for $R_{\alpha}$ found in eq. (89) above, we have:

\displaystyle R_{\alpha}

\displaystyle=\frac{\langle\omega^{\rm eff}\rangle}{2\tilde{\nu}}\left(2\epsilon_{\kappa}^{2}[1+\epsilon_{p}z_{\alpha}^{p}-\epsilon_{d}z_{\alpha}^{d}+\epsilon_{d}^{2}(z_{\alpha}^{d})^{2}-\epsilon_{d}\epsilon_{p}z_{\alpha}^{p}z_{\alpha}^{d}-\epsilon_{\kappa}^{2}]+\mathcal{O}(\epsilon^{5})\right).

(99)

Now we can use the definition of $\epsilon_{\kappa}$ to simplify the expression further:

\displaystyle R_{\alpha}

\displaystyle=\frac{\langle\kappa^{\rm eff}\rangle}{\langle\omega^{\rm eff}\rangle}[1+\epsilon_{p}z_{\alpha}^{p}-\epsilon_{d}z_{\alpha}^{d}+\epsilon_{d}^{2}(z_{\alpha}^{d})^{2}-\epsilon_{d}\epsilon_{p}z_{\alpha}^{p}z_{\alpha}^{d}-\epsilon_{\kappa}^{2}]+\mathcal{O}(\epsilon^{3}).

(100)

This expression makes sense, because to lowest order the resource abundance is just the ratio of supply to depletion.

We can immediately calculate the expectation of various moments. The mean resource abundance is given by

\displaystyle\langle R\rangle

\displaystyle=\frac{\langle\kappa^{\rm eff}\rangle}{\langle\omega^{\rm eff}\rangle}(1+\epsilon_{d}^{2}-\epsilon_{k}^{2})+\mathcal{O}(\epsilon^{3}).

(101)

Similarly, one can express higher order moments as

$\displaystyle\langle R^{2}\rangle$	$\displaystyle=\langle R\rangle^{2}(1+\epsilon_{d}^{2}+\epsilon_{p}^{2})+\mathcal{O}(\epsilon^{3})$
$\displaystyle\langle R^{3}\rangle$	$\displaystyle=\langle R\rangle^{3}(1+3\epsilon_{d}^{2}+3\epsilon_{p}^{2})+\mathcal{O}(\epsilon^{3})$
$\displaystyle\langle R^{4}\rangle$	$\displaystyle=\langle R\rangle^{4}(1+6\epsilon_{d}^{2}+6\epsilon_{p}^{2})+\mathcal{O}(\epsilon^{3})$	(102)

G.4 Average abundance of leaked resources

We now solve for the average resource abundance $\langle R\rangle$ . To proceed, we substitute the expression for $\tilde{\nu}$ in Eq. 79 into Eq. 91 to obtain the intermediate relation

\displaystyle\gamma^{-1}\phi_{N}

\displaystyle=\frac{1}{2}\left\langle 1-\frac{\omega^{\rm eff}_{0}}{\sqrt{(\omega^{\rm eff}_{0})^{2}+4\kappa^{\rm eff}_{0}\tilde{\nu}}}\right\rangle

(103)

Using eq. (98), this can be rewritten as

$\displaystyle\gamma^{-1}\phi_{N}$	$\displaystyle=\frac{1}{2}\left\langle 1-\frac{1}{1+\epsilon_{d}z_{\alpha}^{d}+2\epsilon_{\kappa}^{2}[1+\epsilon_{p}z_{\alpha}^{p}-\epsilon_{d}z_{\alpha}^{d}+\epsilon_{d}^{2}(z_{\alpha}^{d})^{2}-\epsilon_{d}\epsilon_{p}z_{\alpha}^{p}z_{\alpha}^{d}-\epsilon_{\kappa}^{2}]+\mathcal{O}(\epsilon^{5})}\right\rangle$	(104)
	$\displaystyle=\frac{1}{2}\left\langle\epsilon_{d}z_{\alpha}^{d}+2\epsilon_{\kappa}^{2}-\epsilon_{d}^{2}(z_{\alpha}^{d})^{2}\right\rangle+\mathcal{O}(\epsilon^{3})$	(105)
	$\displaystyle=\epsilon_{k}^{2}-{\epsilon_{d}^{2}\over 2}$	(106)

Explicitly this last equation can be written as

\gamma^{-1}\phi_{N}=\frac{\langle\kappa^{\rm eff}\rangle\tilde{\nu}}{\langle\omega^{\rm eff}\rangle^{2}}-\frac{\tilde{\sigma_{d}}^{2}}{2\langle\omega^{\rm eff}\rangle^{2}}+\mathcal{O}(\epsilon^{3}).

(107)

Rearranging this equation gives

\tilde{\nu}=\frac{\langle\omega^{\rm eff}\rangle^{2}\gamma^{-1}\phi_{N}+\frac{1}{2}\tilde{\sigma}_{d}^{2}}{\langle\kappa^{\rm eff}\rangle}.

(108)

Using the expression for $\langle R\rangle$ in Eq. 100, we can simplify this to:

\displaystyle\tilde{\nu}\langle R\rangle=\langle\omega^{\rm eff}\rangle\left(\gamma^{-1}\phi_{N}+\frac{1}{2}\epsilon_{d}^{2}\right)+\mathcal{O}(\epsilon^{3}).

(109)

Combining this with Eq. 79 also yields the relation:

\displaystyle\langle\chi R\rangle

\displaystyle=\frac{\gamma^{-1}\phi_{N}}{\tilde{\nu}}=\frac{\gamma^{-1}\phi_{N}\langle R\rangle}{\langle\omega^{\rm eff}\rangle\left(\gamma^{-1}\phi_{N}+\frac{1}{2}\epsilon_{d}^{2}\right)}.

(110)

Substituting our expression for $\tilde{\nu}\langle R\rangle$ , Eq. 83, and Eq. 101 into eq. (88) gives

$\displaystyle\langle\kappa_{0}^{\mathrm{eff}}\rangle$	$\displaystyle=\gamma^{-1}lD\left(\mu_{c}\langle N\rangle(\langle R\rangle+{R_{E}\over M})+\sqrt{M}\sigma_{c}\langle NxR\rangle\right)$
	$\displaystyle=\gamma^{-1}lD\left(\mu_{c}\langle N\rangle(\langle R\rangle+{R_{E}\over M})+{\phi_{N}\langle R^{2}\rangle\over\langle\chi R\rangle}\right)$
	$\displaystyle=\gamma^{-1}lD\left(\mu_{c}\langle N\rangle(\langle R\rangle+{R_{E}\over M}\right)+\gamma\langle\omega^{\rm eff}\rangle\left(\gamma^{-1}\phi_{N}+\frac{1}{2}\epsilon_{d}^{2}\right){\langle R^{2}\rangle\over\langle R\rangle}$
	$\displaystyle=\gamma^{-1}lD\left(\mu_{c}\langle N\rangle(\langle R\rangle+{R_{E}\over M}\right)+\gamma\langle\omega^{\rm eff}\rangle\left(\gamma^{-1}\phi_{N}+\frac{1}{2}\epsilon_{d}^{2}\right)\langle R\rangle(1+\epsilon_{d}^{2}+\epsilon_{p}^{2}).$	(111)

We can substitute this equation into Eq. 101 to get

	$\displaystyle\langle R\rangle$	$\displaystyle={R_{E}\over M}{q\over 1-q}\left(1+{\epsilon_{d}^{2}\over 1-q}-{(1+\gamma)\epsilon_{k}^{2}\over 1-q}\right)+\mathcal{O}(\epsilon^{3})$
	$\displaystyle q$	$\displaystyle=lD\left(1-{\omega\over\langle\omega_{eff}\rangle}\right)$		(112)

Substituting Eq. 30 into the expression above yields

\langle R\rangle={R_{E}\over M}{q\over 1-q}={\tilde{\kappa}_{E}\over\langle\omega_{eff}\rangle}{q\over 1-q}+\mathcal{O}(\epsilon^{2}).

(113)

One final approximation we can make is consider the limit where depletion is dominated by microbes and we expand the expression above dropping all terms $\mathcal{O}\left({\omega\over\langle\omega_{eff}\rangle}\right)$ . In this case, defining

q_{0}=lD

(114)

and using energy conservation relation Eq. 26 with $\mathrm{w}=1$ we have that

$\displaystyle\langle R\rangle$	$\displaystyle={\tilde{\kappa}_{E}\over\gamma^{-1}\mu_{c}\langle N\rangle}{q_{0}\over 1-q_{0}}+\mathcal{O}(\epsilon^{2})$
	$\displaystyle={\langle Nm\rangle\over g\mu_{c}\langle N\rangle}{q_{0}\over 1-q_{0}}+\mathcal{O}(\epsilon^{2})$
	$\displaystyle={\langle N\rangle\langle m\rangle\over g\mu_{c}\langle N\rangle}{q_{0}\over 1-q_{0}}+\mathcal{O}(\epsilon^{2})$
	$\displaystyle={m\over g\mu_{c}}{q_{0}\over 1-q_{0}}$	(115)

where we have defined the expectation

\langle Nm\rangle={1\over S}\sum_{i}N_{i}m_{i}

(116)

and in going from the second to third line we have used that $\langle Nm\rangle=\langle N\rangle\langle m\rangle$ at this order in $\epsilon$ .

G.5 Second and higher-order moments of resource abundance

In order to calculate higher order moments of resources (see Eq. 102) we need to calculate $\epsilon_{d}^{2}$ and $\epsilon_{p}^{2}$ . Substituting the expressions in Eq. 35 and Eq. 37 into Eq. 94 yields

\displaystyle\epsilon_{d}^{2}

\displaystyle=\frac{\langle(\delta\Phi_{0}^{d})^{2}\rangle}{\langle\Phi_{0}^{d}\rangle^{2}}\times{\langle R\rangle^{2}\over\langle R^{2}\rangle}=\gamma{\sigma_{c}^{2}\over\mu_{c}^{2}}{\langle N^{2}\rangle\over\langle N\rangle^{2}}.

(117)

This expression implies that, as expected, ${\sigma_{c}^{2}\over\mu_{c}^{2}}\sim\epsilon^{2}$ .

In order to calculate $\epsilon_{p}^{2}$ , we start by calculating the average of the production flux. Combining Eq. 42 with Eq. 91 and yields Eq. 106 yields

\langle\Phi_{\alpha}^{p}\rangle=\gamma^{-1}\mu_{c}\langle N\rangle lD(\langle R\rangle+{R_{E}\over M})+\gamma\langle\omega^{eff}\rangle{\langle R^{2}\rangle\over\langle R\rangle}\epsilon_{k}^{2}+\mathcal{O}(\epsilon^{3}).

(118)

To simplify this expression further, we note that $\gamma^{-1}\mu_{c}\langle N\rangle=\langle\omega^{eff}\rangle-\omega$ and then use Eq 112 to get

\langle\Phi_{\alpha}^{p}\rangle=\langle\omega^{eff}\rangle\langle R\rangle+\gamma\langle\omega^{eff}\rangle{\langle R^{2}\rangle\over\langle R\rangle}\epsilon_{k}^{2}+\mathcal{O}(\epsilon^{3}).

(119)

In order to calculate $\sigma_{p}^{2}$ in Eq. 42, we start by simplifying the expressions in Eq. 83 using Eq. 110, Eq. 102, and Eq. 106, yielding:

$\displaystyle\sqrt{M}\sigma_{c}\langle NxR\rangle$	$\displaystyle=\gamma\langle\omega^{eff}\rangle\langle R\rangle\epsilon_{k}^{2}+\mathcal{O}(\epsilon^{3})$	(120)
$\displaystyle\sqrt{M}\sigma_{c}\langle NxR^{2}\rangle$	$\displaystyle=\gamma\langle\omega^{eff}\rangle\langle R\rangle^{2}\epsilon_{k}^{2}+\mathcal{O}(\epsilon^{3})$	(121)
$\displaystyle M\sigma_{c}^{2}\langle NxR^{2}xN\rangle$	$\displaystyle=\gamma^{2}\langle\omega^{eff}\rangle^{2}\langle R\rangle^{2}(1+\gamma)\epsilon_{k}^{2}+\mathcal{O}(\epsilon^{3}).$	(122)

Plugging these into Eq. 42 gives

\sigma_{p}^{2}={\sigma_{D}^{2}\over D^{2}}\langle\omega^{eff}\rangle^{2}\langle R\rangle^{2}\left[q^{2}+q^{2}(\epsilon_{d}^{2}+\epsilon_{p}^{2})+(l^{2}D^{2}(\gamma+1)+2qlD\gamma)\epsilon_{k}^{2}+\mathcal{O}(\epsilon^{3})\right]

(123)

and Eq 46. This yields

\epsilon_{p}^{2}={\sigma_{D}^{2}\over D^{2}}q^{2}+\mathcal{O}(\epsilon^{3})

(124)

where we have noted this expression implies that as expected ${\sigma_{D}^{2}\over D^{2}}\sim\epsilon^{2}$ . Thus, to this order in $\epsilon_{p}^{2}$ we have

\displaystyle\epsilon_{p}^{2}

\displaystyle={\sigma_{D}^{2}\over D^{2}}q^{2}

(125)

Appendix H Summary of results

In this section, we summarize our results from the calculation. We start by defining the parameters

$\displaystyle q$	$\displaystyle=lD\left(1-{\omega\over\langle\omega_{eff}\rangle}\right)$
$\displaystyle\epsilon_{d}^{2}$	$\displaystyle=\gamma{\sigma_{c}^{2}\over\mu_{c}^{2}}{\langle N^{2}\rangle\over\langle N\rangle^{2}}$
$\displaystyle\epsilon_{p}^{2}$	$\displaystyle={\sigma_{D}^{2}\over D^{2}}q^{2}$
$\displaystyle\epsilon_{k}^{2}$	$\displaystyle=\gamma^{-1}\phi_{N}+{1\over 2}\epsilon_{d}^{2}$	(126)

The average resource abundance is given by

$\displaystyle\langle R\rangle$	$\displaystyle={R_{E}\over M}{q\over 1-q}\left(1+{\epsilon_{d}^{2}\over 1-q}-{(1+\gamma)\epsilon_{k}^{2}\over 1-q}\right)+\mathcal{O}(\epsilon^{3})$
	$\displaystyle={\tilde{k}_{E}\over\langle\omega^{eff}\rangle}{q\over 1-q}\left(1+{\epsilon_{d}^{2}\over 1-q}-{(1+\gamma)\epsilon_{k}^{2}\over 1-q}\right)+\mathcal{O}(\epsilon^{3})$
	$\displaystyle={m\over g\mu_{c}}{q\over 1-q}\left(1+{\epsilon_{d}^{2}\over 1-q}-{(1+\gamma)\epsilon_{k}^{2}\over 1-q}\right)+\mathcal{O}(\epsilon^{3}).$	(127)

The higher order resource moments are just

$\displaystyle\langle R^{2}\rangle$	$\displaystyle=\langle R\rangle^{2}(1+\epsilon_{d}^{2}+\epsilon_{p}^{2})+\mathcal{O}(\epsilon^{3})$
$\displaystyle\langle R^{3}\rangle$	$\displaystyle=\langle R\rangle^{3}(1+3\epsilon_{d}^{2}+3\epsilon_{p}^{2})+\mathcal{O}(\epsilon^{3})$
$\displaystyle\langle R^{4}\rangle$	$\displaystyle=\langle R\rangle^{4}(1+6\epsilon_{d}^{2}+6\epsilon_{p}^{2})+\mathcal{O}(\epsilon^{3})$
$\displaystyle q$	$\displaystyle=lD\left(1-{\omega\over\langle\omega_{eff}\rangle}\right)$
$\displaystyle\epsilon_{d}^{2}$	$\displaystyle=\gamma^{-1}{\sigma_{c}^{2}\over\mu_{c}^{2}}{\langle N^{2}\rangle\over\langle N\rangle^{2}}$
$\displaystyle\epsilon_{p}^{2}$	$\displaystyle={\sigma_{D}^{2}\over D^{2}}q^{2}$
$\displaystyle\epsilon_{k}^{2}$	$\displaystyle=\gamma^{-1}\phi_{N}+{1\over 2}\epsilon_{d}^{2}.$	(128)

The full resource distribution is given by

\displaystyle R_{\alpha}

\displaystyle=\frac{\langle R\rangle}{2\gamma^{-1}\phi_{N}+\epsilon_{d}^{2}}\left[\sqrt{(1+\epsilon_{d}z_{\alpha}^{d})^{2}+2(2\gamma^{-1}\phi_{N}+\epsilon_{d}^{2})(1+\epsilon_{p}z_{\alpha}^{p})}-1-\epsilon_{d}z_{\alpha}^{d}\right],

(129)

where are before $z_{\alpha}^{d}$ and $z_{\alpha}^{p}$ are standard random normal variables.

In addition, we have expressions for the species abundances.

$\displaystyle\phi_{N}$	$\displaystyle=w_{0}(\Delta_{N})$
$\displaystyle\langle N\rangle$	$\displaystyle=\frac{\sigma_{N}}{(1-l)g\sigma_{c}^{2}\langle\chi R\rangle}w_{1}(\Delta_{N})$
$\displaystyle\langle N^{2}\rangle$	$\displaystyle=\left(\frac{\sigma_{N}}{(1-l)g\sigma_{c}^{2}\langle\chi R\rangle}\right)^{2}w_{2}(\Delta_{N})$
$\displaystyle w_{n}(\Delta)$	$\displaystyle=\int_{-\Delta}^{\infty}{1\over\sqrt{2\pi}}e^{-z^{2}\over 2}(z+\Delta)^{n}$
$\displaystyle\Delta_{N}$	$\displaystyle=\frac{(1-l)\mu_{c}g\langle R\rangle-m}{\sigma_{N}}$
$\displaystyle\sigma_{N}^{2}$	$\displaystyle=\sigma_{m}^{2}+(1-l)^{2}\sigma_{c}^{2}g^{2}\langle R\rangle^{2}(1+\epsilon_{d}^{2}+\epsilon_{p}^{2}).$	(130)

Finally, we have

	$\displaystyle\tilde{\nu}\langle R\rangle$	$\displaystyle=\langle\omega^{\rm eff}\rangle\left(\gamma^{-1}\phi_{N}+\frac{1}{2}\epsilon_{d}^{2}\right)+\mathcal{O}(\epsilon^{3})$
	$\displaystyle\langle\chi R\rangle$	$\displaystyle=\frac{\gamma^{-1}\phi_{N}}{\tilde{\nu}}=\frac{\gamma^{-1}\phi_{N}\langle R\rangle}{\langle\omega^{\rm eff}\rangle\left(\gamma^{-1}\phi_{N}+\frac{1}{2}\epsilon_{d}^{2}\right)}.$		(131)

Appendix I Species Packing Bounds

The expressions above allow us to define a simple bound for the fraction of metabolic niches that an ecosystem can occupy Cui et al. (2020). If we denote the number of surviving species by $S^{*}$ and the number of metabolites by $M$ , then in the thermodynamic limit the species-packing fraction is just $S^{*}/M=\gamma^{-1}\phi_{N}$ . To derive the bound, we start with Eq. 91 which reads

\displaystyle\langle\chi R\rangle

\displaystyle=\frac{1}{2\tilde{\nu}}\left\langle 1-\frac{\omega^{\rm eff}_{0}}{\sqrt{(\omega^{\rm eff}_{0})^{2}+4\kappa^{\rm eff}_{0}\tilde{\nu}}}\right\rangle.

(132)

Substituting Eq. 79 for $\tilde{\nu}$ into the expression yields

\displaystyle\gamma^{-1}\phi_{N}

\displaystyle=\frac{1}{2}\left\langle 1-\frac{\omega^{\rm eff}_{0}}{\sqrt{(\omega^{\rm eff}_{0})^{2}+4\kappa^{\rm eff}_{0}\tilde{\nu}}}\right\rangle.

(133)

Since the second term is always positive we have that

\gamma^{-1}\phi_{N}\leq\frac{1}{2},

(134)

proving our bound. Note that in proving this bound, we have made no approximations beyond replica-symmetry and therefore expect this to be exact in the thermodynamic limit.

Appendix J Fraction of surviving species depends only on competition and not leakage rate

In this section, we show that to leading order in $\epsilon$ (see Eq. 94), the fraction of surviving species depends only on the competition between species ( $\mu_{c}$ and $\sigma_{c}$ ) but not the leakage rate $l$ when $\sigma_{m}^{2}$ , the fluctuations in the maintenance costs, can be ignored. Technically, this is a statement that $\Delta_{N}$ , and hence $\phi_{N}$ , depends on the variance over the mean squared of consumer preferences $\sigma_{c}^{2}/\mu_{c}^{2}$ but not the leakage parameter $l$ . There are many ways to show this and here we present one particularly simple derivation of this fact.

Our starting point is Eq. 130:

\langle N\rangle=\frac{\sigma_{N}}{(1-l)g\sigma_{c}^{2}\langle\chi R\rangle}w_{1}(\Delta_{N}).

(135)

We then use Eq. eq:chiR2 to rewrite this as

	$\displaystyle\langle N\rangle$	$\displaystyle={\sigma_{N}\over(1-l)g\sigma_{c}^{2}\langle\chi R\rangle}$
		$\displaystyle={\sigma_{N}\langle\omega^{\rm eff}\rangle\epsilon_{k}^{2}\over(1-l)g\sigma_{c}^{2}\gamma^{-1}\phi_{N}\langle R\rangle},$		(136)

where in the second line we have used Eq. 110 and Eq. 126 for $\epsilon_{k}^{2}$ . Now, using the fact that when $\sigma^{2}=0$ we have from Eq. 130 that

\sigma_{N}=(1-l)\sigma_{c}g\langle R\rangle(1+{1\over 2}\epsilon_{d}^{2}+{1\over 2}\epsilon_{p}^{2})+\mathcal{O}(\epsilon^{3}),

(137)

the above can be rewritten as

$\displaystyle\langle N\rangle$	$\displaystyle={\langle\omega^{\rm eff}\rangle\epsilon_{k}^{2}\over\sigma_{c}\gamma^{-1}\phi_{N}}+\mathcal{O}(\epsilon^{3})$
$\displaystyle\langle N\rangle$	$\displaystyle=\gamma^{-1}\langle N\rangle{\mu_{c}\over\sigma_{c}}\left(1+{\epsilon_{d}^{2}\over 2\gamma^{-1}\phi_{N}}\right)+\mathcal{O}(\epsilon^{3})$
$\displaystyle\langle N\rangle$	$\displaystyle=\gamma^{-1}\langle N\rangle{\mu_{c}\over\sigma_{c}}\left(1+{\sigma_{c}^{2}\over\mu^{2}}{\langle N^{2}\rangle\over 2\phi_{N}\langle N\rangle^{2}}\right)+\mathcal{O}(\epsilon^{3})$
$\displaystyle\langle N\rangle$	$\displaystyle=\gamma^{-1}\langle N\rangle{\mu_{c}\over\sigma_{c}}\left(1+{\sigma_{c}^{2}\over\mu^{2}}{w_{2}(\Delta_{N})\over 2w_{0}(\Delta_{N})w_{1}(\Delta_{N})^{2}}\right)+\mathcal{O}(\epsilon^{3})$	(138)

Canceling the $\langle N\rangle$ from both sides yields

1=\gamma^{-1}{\mu_{c}\over\sigma_{c}}\left(1+{\sigma_{c}^{2}\over\mu^{2}}{w_{2}(\Delta_{N})\over 2w_{0}(\Delta_{N})w_{1}(\Delta_{N})^{2}}\right)+\mathcal{O}(\epsilon^{3}),

(139)

showing that $\Delta_{N}$ and hence $\phi_{N}$ only depends on ${\sigma_{c}^{2}\over\mu_{c}}$ at this order.

Appendix K Details of numerical simulations

K.1 Numerical Simulations

We now give an overview of how we performed numerical simulations. All code can be found on the corresponding Github page. To numerically simulate the communities, we used the Community Simulator package to find steady-states with indicated parameters Marsland et al. (2020a). We used the parameters: $q=0$ (generalist consumers), supply = ’external’, regulation= ’independent’, response=’type I’ , and sparsity $s=0.01$ . These choices ensure that the dynamics follow Eq. 12. We set the number of species $S=600$ , the number of metabolites equal $M=300$ , the average maintenance cost $m=1$ with standard deviation $\sigma_{m}=0.1$ . A single resource was supplied externally at a rate $\tilde{\kappa}_{E}=\kappa_{E}/M=100$ and $\omega=1$ . The leakage rate was chosen to be $l=0.8$ if fixed, or varied over $10$ equal intervals from $l=0.7-0.95$ as indicated in figures.

The consumption coefficients were chose from binary distribution with $c_{i\alpha}\in\{0,c_{1}\}$ with probability having non-zero coefficient $c_{1}$ given by $p=\mu_{c}/(c_{1}M)$ . This was done to ensure that the average consumption rate $\langle c_{i\alpha}\rangle=\mu_{c}/M$ , with $\mu_{c}=1$ . Under these assumptions $\sigma_{c}=\sqrt{Mc_{1}p(1-p)}$ . In order to vary competition, $c_{1}$ was varied over $10$ equal interval from $c_{1}=1/M-10/M$ with $M=300$ . Default choice for $c_{1}=9/M=0.03$ . See section “Sampling parameters” in Marsland et al. (2020a) for details and Github for code.

Averages $\phi_{N}$ , $\langle N\rangle$ , $\langle N^{2}\rangle$ , $\langle R\rangle$ , and $\langle R^{2}\rangle$ reflect averages over 28 independent realizations. In order to ensure numerical stability, we treat species as abundances below $10^{-3}$ as going extinct.

K.2 Numerical solution of cavity equations

To numerically check the solutions to our self-consistent cavity equations, we used the expressions summarized in section H. In particular, given the parameters we solved for $\Delta_{N}$ numerically. To do so, we used the second line of Eq. 138, Eq. 128, and Eq. 130 to get:

$\displaystyle\phi_{N}{\sigma_{c}\over\mu_{c}}$	$\displaystyle=\left(\gamma^{-1}\phi_{N}+{\epsilon_{d}^{2}\over 2}\right)$
$\displaystyle\phi_{N}{\sigma_{c}\over\mu_{c}}$	$\displaystyle=\left(\gamma^{-1}\phi_{N}+{\sigma_{c}^{2}\over\mu_{c}^{2}}{\langle N\rangle^{2}\over\langle N\rangle^{2}}\right)$
$\displaystyle w_{0}(\Delta_{N}){\sigma_{c}\over\mu_{c}}$	$\displaystyle=\left(\gamma^{-1}w_{0}(\Delta_{N})+{\sigma_{c}^{2}\over\mu_{c}^{2}}{w_{2}(\Delta_{N})\over w_{1}(\Delta_{N})^{2}}\right)$	(140)

This allowed us define the cost function

\displaystyle\mathcal{C}(\Delta_{N})=\left(w_{0}(\Delta_{N}){\sigma_{c}\over\mu_{c}}-\left(\gamma^{-1}w_{0}(\Delta_{N})+{\sigma_{c}^{2}\over\mu_{c}^{2}}{w_{2}(\Delta_{N})\over w_{1}(\Delta_{N})^{2}}\right)\right)^{2},

(141)

which we can minimize to solve for $\Delta_{N}$ using the “minimize scalar” function from the scipy.opitmize package.

Once we solve for $\Delta_{N}$ , we can use Eq. 126 along with Eq. 129 to numerically calculate the predicted resource distribution:

\displaystyle R_{\alpha}

\displaystyle=\frac{\langle R\rangle}{2\gamma^{-1}\phi_{N}+\epsilon_{d}^{2}}\left[\sqrt{(1+\epsilon_{d}z_{\alpha}^{d})^{2}+2(2\gamma^{-1}\phi_{N}+\epsilon_{d}^{2})(1+\epsilon_{p}z_{\alpha}^{p})}-1-\epsilon_{d}z_{\alpha}^{d}\right].

(142)

To do so, we randomly sample from this distribution $10^{6}$ times and then plot the resulting pdf.

We can also use Eq. 71 to generate the distribution of the species abundances. In particular, this equation states that the distribution of surviving species follows a truncated Gaussian distribution with a fraction $\phi_{N}$ of the species surviving. To proceed, note that from Eq. 130 that we can define the quantity

\Sigma={\langle N\rangle\over w_{1}(\Delta_{N})}=\frac{\sigma_{N}}{(1-l)g\sigma_{c}^{2}\langle\chi R\rangle}

(143)

In terms of $\Sigma$ , we can rearrange Eq. 71 to derive that surviving species will be described by the distribution

P_{N}^{sur}(N)={1\over\sqrt{2\pi}\phi_{N}\Sigma}e^{-0.5({N\over\Sigma}-\Delta_{N})^{2}}.

(144)

$\Delta_{N}$ is obtained by minimizing $\mathcal{C}(\Delta_{N})$ as discussed above. To obtain $\langle N\rangle$ to leading order in $\epsilon$ , we use energy conservation Eq. Eq:EnergyCons:

$\displaystyle\kappa_{E}$	$\displaystyle\approx\sum_{i}\frac{N_{i}m_{i}}{\mathrm{w}g}$
$\displaystyle{\kappa_{E}\over S}$	$\displaystyle\approx{1\over S}\sum_{i}\frac{N_{i}m_{i}}{\mathrm{w}g}$
$\displaystyle{\kappa_{E}\over S}$	$\displaystyle\approx\frac{\langle Nm\rangle}{\mathrm{w}g}$
$\displaystyle{\kappa_{E}\over S}$	$\displaystyle\approx\frac{\langle N\rangle\langle m\rangle}{\mathrm{w}g},$	(145)

which can be rearrange to yield

\langle N\rangle\approx{\mathrm{w}g\kappa_{E}\over mS},

(146)

allowing us to calculate $\Sigma$ .

As a further check the prediction of the cavity method and our approximations we use the value of $\Delta_{N}$ obtained from the procedure above along with Eq. 130 to compute

{\langle N^{2}\rangle\over\langle N\rangle^{2}}={w_{2}(\Delta_{N})\over w_{1}(\Delta_{N})}.

(147)

To check that the average resource abundance $\langle R\rangle$ only depend on the leakage rate in the limit of fast dilution (i.e. when $q\approx q_{0}=l$ ), we plotted the leading order contribution to Eq. 127:

\langle R\rangle\approx{m\over g\mu_{c}}{l\over 1-l}

(148)

versus the average obtained from direct numerical simulation.

Cross-feeding shapes both competition and cooperation in microbial ecosystems

Abstract

References

Appendix A Basic Setup

Appendix B Energy Conservation

Appendix C Externally Supplied Resource

Appendix D Means, variances, and covariances for fluxes of self-generated resources

D.1 Depletion flux

D.2 Production flux

D.3 Covariance between production and depletion flux

Appendix E Growth flux

Appendix F Cavity Corrections

F.1 Setting-up the perturbation theory

F.2 Computing the TAP corrections

Appendix G Self-consistency equation

G.1 Equations for species

G.2 Equations for resources

G.3 Solving the equations for resources

G.4 Average abundance of leaked resources

G.5 Second and higher-order moments of resource abundance

Appendix H Summary of results

Appendix I Species Packing Bounds

Appendix J Fraction of surviving species depends only on competition and not leakage rate

Appendix K Details of numerical simulations

K.1 Numerical Simulations

K.2 Numerical solution of cavity equations

Appendix L Additional Figures and Numerical Checks