#  Phase transition to chaos in complex ecosystems with non-reciprocal
species-resource interactions

Emmy Blumenthal<EMAIL_ADDRESS>Jason W. Rocks<EMAIL_ADDRESS>Pankaj Mehta
<EMAIL_ADDRESS>Department of Physics, Boston University, Boston, MA 02215, USA
Faculty of Computing and Data Science, Boston University, Boston, MA 02215,
USA

###### Abstract

Non-reciprocal interactions between microscopic constituents can profoundly
shape the large-scale properties of complex systems. Here, we investigate the
effects of non-reciprocity in the context of theoretical ecology by analyzing
a generalization of MacArthur’s consumer-resource model with asymmetric
interactions between species and resources. Using a mixture of analytic cavity
calculations and numerical simulations, we show that such ecosystems
generically undergo a phase transition to chaotic dynamics as the amount of
non-reciprocity is increased. We analytically construct the phase diagram for
this model and show that the emergence of chaos is controlled by a single
quantity: the ratio of surviving species to surviving resources. We also
numerically calculate the Lyapunov exponents in the chaotic phase and
carefully analyze finite-size effects. Our findings show how non-reciprocal
interactions can give rise to complex and unpredictable dynamical behaviors
even in the simplest ecological consumer-resource models.

††preprint: APS/123-QED

Many complex systems operate out of equilibrium where components generically
interact non-reciprocally. Significant current research aims to untangle the
implications of non-reciprocal interactions for self-organization and pattern
formation. While much progress has been made towards understanding non-
reciprocity in systems composed of a few types of species or fields, the
consequences of non-reciprocity in more complex systems composed of many
interacting components are less clear [1, 2, 3, 4].

Large, diverse ecosystems with many types of species and resources provide a
natural setting for exploring this open problem. Over the last decade,
researchers have adapted methods from the statistical physics of disordered
systems (e.g., replicas, the cavity method, Random Matrix Theory) to analyze
such ecosystems [5, 6, 7, 8, 9, 10, 11, 12]. Much of this work has focused on
systems with reciprocal interactions in which dynamics are often implicitly
governed by an optimization function and reach a fixed point [13].

One notable exception are recent studies of the random Generalized
Lotka–Volterra model in which species interact non-reciprocally [14, 15, 16,
17, 18, 19]. These systems can exhibit novel behaviors such as dynamic
fluctuations and chaos, including unpredictable “boom-and-bust” dynamics where
low-abundance species suddenly bloom to high abundance [20]. These
observations suggest that non-reciprocal interactions can qualitatively change
ecological dynamics in species-only models. However, the generalization of
these observations to more complex ecosystems with multiple trophic layers or
environmentally-mediated interactions remains unexplored.

Here, we introduce a generalization of the classic MacArthur Consumer Resource
Model (MCRM) that includes non-reciprocal interactions between species and
resources. Consumer-resource models, first introduced by MacArthur and Levins
[21, 22, 23], have played a foundational role in modern theoretical ecology
and undergird many powerful theoretical frameworks for understanding
ecological competition, including contemporary niche theory and Tilman’s R*
principle [24, 25].

Theoretical Setup. We consider an ecosystem with ${i=1,\ldots,S}$ species
which may consume ${\alpha=1,\ldots,M}$ distinct self-replenishing resources
with dynamics governed by the equations,

$\displaystyle\derivative{N_{i}}{t}$
$\displaystyle=N_{i}\quantity(\sum_{\alpha=1}^{M}c_{i\alpha}R_{\alpha}-m_{i}),$
(1) $\displaystyle\derivative{R_{\alpha}}{t}$
$\displaystyle=R_{\alpha}\quantity(K_{\alpha}-R_{\alpha})-\sum_{i=1}^{S}N_{i}e_{i\alpha}R_{\alpha},$
(2)

where $N_{i}$ is the population size of species $i$, $R_{\alpha}$ is the
abundance of resource $\alpha$, $c_{i\alpha}$ is the relative consumption
preference of species $i$ for resource $\alpha$, $e_{i\alpha}$ describes the
impact of species $i$ on resource $\alpha$, $m_{i}$ is the natural mortality
rate of species $i$, and $K_{\alpha}$ is the carrying capacity of resource
$\alpha$ in the absence of consumption. We call this model the asymmetric
MacArthur Consumer Resource Model (aMCRM) with a schematic provided in Fig. 1.
When $e_{i\alpha}=c_{i\alpha}$ the species-resource interactions become
reciprocal, or symmetric, and the aMCRM reduces to the classical MacArthur
Consumer Resource Model (MCRM).

To develop intuition for the role of non-reciprocity in the aMCRM, we consider
the limit where the resource dynamics are fast and the resource abundances
become entrained to species dynamics. In this case, we take the RHS of Eq. (2)
to be zero and solve to find
${R_{\alpha}=\max\\{0,K_{\alpha}-\sum_{i}N_{i}e_{i\alpha}R_{\alpha}\\}}$.
Substituting this result into the equation for species dynamics yields an
effective Lotka–Volterra equation,

$\displaystyle\derivative{N_{i}}{t}=$ $\displaystyle
N_{i}\quantity(\kappa_{i}-\sum_{j=1}^{S}A_{ij}N_{j}),$ (3)
$\displaystyle\kappa_{i}=$
$\displaystyle\sum_{\alpha=1}^{M}c_{i\alpha}K_{\alpha}-m_{i},$ $\displaystyle
A_{ij}=$
$\displaystyle\sum_{\alpha=1}^{M}c_{i\alpha}e_{j\alpha}\Theta(R_{\alpha}),$

where $\kappa_{i}$ is the effective carrying capacity for species $i$ and
$A_{ij}$ is the effective species-species interaction matrix, encoding how
species $j$ impacts species $i$ ($\Theta$ is the Heaviside function). Although
typically not quantitatively accurate, this approximation provides useful
qualitative insight into the nature of the non-reciprocal interactions.

In MacArthur’s original consumer-resource model, impacts and benefits are
identical, ${e_{i\alpha}=c_{i\alpha}}$. In this case, $A_{ij}$ is symmetric,
all interactions are reciprocal, the ecosystem has a unique fixed point, and
the resulting steady state can be derived using an optimization principle
[13]. Such behavior is expected because choosing ${c_{i\alpha}=e_{i\alpha}}$
implicitly assumes that each species consumes resources proportional to the
marginal utility conferred to that species (in the context of game theory and
microeconomics, this is a “rational strategy”). When the resource-species
interactions are non-reciprocal, ${e_{i\alpha}\neq c_{i\alpha}}$, $A_{ij}$ is
no longer symmetric, the resulting dynamics can no longer be described using
an optimization principle, and there is no guarantee that the dynamics will
reach a stable fixed point.

Figure 1:  Schematic of the asymmetric MacArthur Consumer Resource Model
(aMCRM). Species $i$ benefits with relative weight $c_{i\alpha}$ from
consuming resource $\alpha$ and impacts the abundance of the resource with
relative weight $e_{i\alpha}$.

Thermodynamic Limit. To investigate the aMCRM, we work in the thermodynamic
limit where the numbers of species $S$ and resources $M$ become very large
while their ratio $M/S$ is held fixed. We assume that parameters are drawn
randomly from a fixed distribution analogous to quenched disorder. To ensure a
proper thermodynamic limit, parameters are drawn as follows:

$\begin{gathered}K_{\alpha}=K+\sigma_{K}\delta K_{\alpha},\quad
m_{i}=m+\sigma_{m}\delta m_{i},\\\
c_{i\alpha}=\frac{\mu_{c}}{M}+\frac{\sigma_{c}}{\sqrt{M}}d_{i\alpha},\\\
e_{i\alpha}=\frac{\mu_{e}}{M}+\frac{\sigma_{e}}{\sqrt{M}}\quantity(\rho
d_{i\alpha}+\sqrt{1-\rho^{2}}x_{i\alpha})\end{gathered}$ (4)

where $\delta K_{\alpha},\delta m_{i},d_{i\alpha},x_{i\alpha}$ are independent
standard random variables (i.e., zero mean and unit variance) and $|\rho|\leq
1$ is the interaction reciprocity parameter. For simplicity, we take
$\mu_{c}=\mu_{e}\equiv\mu$ and $\sigma_{c}=\sigma_{e}\equiv\sigma$ in all
figures and simulations. The central limit theorem ensures that, in the
thermodynamic limit, our results are agnostic to the exact form of the
underlying distributions and depend only on first and second moments.
Therefore, we sample all parameters from normal distributions unless otherwise
specified.

With this parameterization, $\rho$ controls the level of reciprocity of
species-resource interactions through the correlation of consumption benefits
and impacts:

$\mathrm{Cor}(c_{i\alpha},e_{j\beta})=\rho\;\delta_{ij}\delta_{\alpha\beta}.$
(5)

When $\rho=1$, the aMCRM reduces to the fully symmetric MCRM; when $\rho=0$,
the aMCRM models completely non-reciprocal species-resource interactions. By
tuning $\rho$, we can systematically explore the effects of non-reciprocity.

Cavity Method. Just as in the original MCRM, we can analytically calculate the
thermodynamic-limit behavior using the cavity method [11, 12, 26, 27]. Unlike
replicas, the cavity method does not require the existence of an energy
function and therefore can be extended to the aMCRM. We assume dynamics are
self-averaging and described by a replica-symmetric ansatz. Using this ansatz,
we derive self-consistent mean-field equations for the fraction of surviving
species, the fraction of non-depleted resources, the first and second moments
of the steady-state species and resource abundances, and the trace of two
relevant susceptibility matrices (see Appendix A for detailed derivations). As
seen in Figs. A6 and A8, numerical simulations and analytical predictions
agree remarkably well for moderate non-reciprocity.

Figure 2:  Example dynamics of the aMCRM in a community of $S=M=256$ species
and resources. Left: dynamics in the stable phase; species-resource
interactions are nearly reciprocal. Right: dynamics in the dynamic phase;
species-resource interactions are less reciprocal. The parameter values for
the stable-phase and dynamic-phase simulations are respectively marked with a
circle and star in Fig. 3(a).  Figure 3:  Phase diagram of the aMCRM and
diverging susceptibility. (a) Heatmap of the fraction of simulations which
reached steady state in finite simulation time for various values of $\rho$,
the level of reciprocity of species-resource interactions, and $\sigma$, the
magnitude of fluctuations in species-resource interactions. Overlain is the
cavity method-calculated phase boundary. (b) Variances of susceptibilities of
mean-field species and resources as a function of $\rho$, with $\sigma$ fixed
at the value indicated by the dashed line in (a).

Transition to Dynamic Phase. Without reciprocal interactions, the aMCRM has no
guarantee of reaching a steady state. In fact, we find that when the
interaction reciprocity $\rho$ is less than a critical $\rho^{\star}$, the
aMCRM exhibits a phase transition from a unique self-averaging steady state to
a chaotic dynamic phase. Fig. 2 shows numerical simulations of typical
resource and species dynamics observed in each phase (see Appendix D for
simulation details [28, 29, 30, 31, 32, 33]).

Using the cavity method, we can analytically compute the phase boundary
between the stable and dynamic phases [26]. We perturb the non-zero steady-
state species and resource abundances, ${N_{i}\to
N_{i}+\varepsilon\eta_{i}^{(N)}}$ and $R_{\alpha}\to
R_{\alpha}+\varepsilon\eta_{\alpha}^{(R)}$, where $\varepsilon$ is a small
parameter and $\eta_{i}^{(N)},\eta_{\alpha}^{(R)}$ are independent standard
random variables, and calculate the susceptibilities
$\derivative*{N_{i}}{\varepsilon}$, $\derivative*{R_{\alpha}}{\varepsilon}$.
Because of the disordered nature of the perturbation, the expectations of the
first moments of the susceptibilities are zero, but the second moments,
$\expectationvalue*{(\derivative*{N_{i}}{\varepsilon})^{2}}$,
$\expectationvalue*{(\derivative*{R_{\alpha}}{\varepsilon})^{2}}$, are non-
zero (see Appendix B for details).

The phase transition to the dynamic phase is signaled by the divergence of the
these susceptibilities’ second moments (see Fig. 3). Surprisingly, we find
that $\rho^{\star}$, the critical value marking the phase transition to chaos,
depends on model parameters only through the species packing fraction, the
ratio of surviving species to non-depleted resources, via the expression (see
Appendix B):

$\left(\rho^{\star}\right)^{2}=\frac{(\text{\\# of surviving
species})}{(\text{\\# of non-depleted resources})}.$ (6)

When $\rho<\rho^{\star}$ the ecosystem undergoes a phase transition to chaos.
As the number of surviving species and non-depleted resources are fixed by
model parameters, the above equation defines a co-dimension-one phase boundary
in the parameter space. Beyond this boundary in the dynamic phase, the second
moments of the susceptibilities become negative, indicating that the replica-
symmetric ansatz no longer holds, and its results are unstable to any
perturbation.

Fig. 3(a) shows a phase diagram overlain on a heatmap of the fraction of
simulations that reach steady state within a chosen finite runtime. We
highlight the locations of the simulations in the stable and dynamic phases in
Fig. 2 with a circle and a star, respectively. In Fig. 3(b), we plot the
second moments of the susceptibilities as a function of $\rho$ with fixed
$\sigma$ along the slice of phase space indicated by the dashed line in Fig.
3(a). The susceptibilities’ variances diverge at the phase transition and
become invalidly negative in the dynamic phase. As the phase transition is
approached, the fraction of simulations that reach steady state in a finite
simulation time sharply decreases. An alternative phase diagram with
parameters drawn from uniform distributions is shown in Fig. A10.

Finally, we note that for certain choices of parameters, the replica-symmetric
self-consistent equations do not have a solution. This transition to
infeasibility has an interesting interpretation but is not physically realized
because it occurs within the dynamic phase where the replica symmetric
solution is unstable (see Appendix A.5).

Chaos. In order to better understand the transition to chaos, we numerically
computed the maximal Lyapunov exponent $\lambda_{1}$ of the aMCRM in the
dynamic and stable phases using the “H2” method of Geist [34, 35, 36, 37]. The
maximal Lyapunov exponent characterizes how quickly trajectories from nearby
initial conditions diverge (positive exponent) or converge (negative
exponent). As seen in Fig. 4(a), in the dynamic phase, $\lambda_{1}>0$, while
in the stable phase, $\lambda_{1}<0$. For the parameters used in Fig. 2,
$|\lambda_{1}|\approx 5\times 10^{-3}$, indicating that the divergence or
convergence of nearby trajectories occurs on a timescale of
$\lambda_{1}^{-1}=2\times 10^{2}$ time units. We further confirmed the
existence of chaos by analyzing the generalized alignment index (GALI) which
measures how a volume element formed by tangent vectors to a trajectory
changes over time [38, 37, 36] (see Fig. C14). Further details are given in
Appendix C.

Figure 4:  Chaos in the dynamic phase of the aMCRM. (a) Dot plot of
$\lambda_{1}$s, maximal Lyapunov exponents, for simulations classified by
whether they reach a steady state for various values of $\rho$. (b) Two
trajectories (red and blue) with slightly different initial conditions in the
dynamic phase of the aMCRM. A species and a resource are highlighted to
emphasize the chaotic dynamics; all other species and resources are shown at
low opacity for clarity. The units of time are given by the inverse of the
maximal Lyapunov exponent, $\lambda_{1}^{-1}=190$.

A direct signature of chaotic dynamics is high sensitivity to initial
conditions as observed in Fig. 4(b). The red and blue lines show the simulated
trajectory of a single species (top) and resource (bottom) started from
initial conditions with slight differences. Initially, the trajectories are
almost identical before diverging from each other significantly after a few
Lyapunov times.

Finite-Size Effects. Like most phase transitions, the transition between the
stable and dynamic phases is a thermodynamic-limit phenomenon. In small
ecosystems, the aMCRM may approach steady state even when in the dynamic phase
due to finite-size effects. As a result, it is not clear in Fig. 3 what the
true probability of steady state is in the thermodynamic limit. In Appendix E,
we quantify these affects by performing a numerical analysis to extrapolate
the steady-state probabilities to infinite system size for each of the two
points highlighted in Fig. 3. For both sets of parameters, we measure the
distribution of steady-state times for many simulations for a variety of
system sizes. Using a custom method based on maximum-likelihood estimation, we
then perform a finite-size scaling collapse on these distributions, allowing
us to approximately determine the steady-state probabilities as a function of
system size. Our scaling collapses provide strong evidence that the
probability of reaching steady state in the thermodynamic limit approaches
exactly zero in the dynamic phase and one in the stable phase.

Discussion. In this letter, we analyzed the effects of non-reciprocal species-
resource interactions on the stability of ecosystems. We introduced the
asymmetric MacArthur Consumer Resource Model (aMCRM), a generalization of the
MacArthur Consumer Resource Model (MCRM). Using the cavity method, we
identified a phase transition between a stable phase in which a unique,
uninvadable, self-averaging steady state exists and a dynamic phase with
chaotic fluctuations. Remarkably, the phase boundary depends on model
parameters only through the species-packing ratio—the ratio of surviving
species to non-depleted resources.

We found that the chaotic regime is generic and occurs robustly, in contrast
with some recent works on Generalized Lotka–Volterra models with completely
antisymmetric interactions where chaos is often hard to nucleate [14, 15]. In
addition, the chaotic dynamics in consumer-resource models generically occurs
when the systems are well below the competitive exclusion bound, while the
dynamics in Generalized Lotka–Volterra systems can violate the competitive
exclusion principle.

Collectively, these works suggest that non-reciprocal interactions can lead to
complex, chaotic dynamics in systems with many different types of
species/fields. In particular, like Generalized Lotka–Volterra models, we also
find that species and resources often jump rapidly between low and high
abundances. In the future, it will be interesting to see if the methods
developed in Ref. 20 in the context of Lotka–Volterra systems generalize to
explain such boom-and-bust dynamics in consumer-resource models. Finally, it
will also be interesting to understand these phenomena in the context of
ecological processes such as immigration, alternative resource dynamics [39],
the addition of network and metabolic structure into interactions [40, 41,
42], the inclusion of additional trophic structure [43], and spatial and
temporal structure [44].

Acknowledgements. We would like to thank Zhijie (Sarah) Feng, Claudio Chamon,
and Chris Laumann for useful discussions. Additionally, we thank the Boston
University Research Computing Services for managing computational resources.
This work was funded by NIH NIGMS R35GM119461 to P.M. and the Boston
University Undergraduate Research Opportunities Program to E.B.

## References

  * Fruchart  _et al._ [2021] M. Fruchart, R. Hanai, P. B. Littlewood, and V. Vitelli, Non-reciprocal phase transitions, Nature 592, 363 (2021).
  * Ivlev  _et al._ [2015] A. V. Ivlev, J. Bartnick, M. Heinen, C.-R. Du, V. Nosenko, and H. Löwen, Statistical mechanics where newton’s third law is broken, Phys. Rev. X 5, 011035 (2015).
  * Sompolinsky and Kanter [1986] H. Sompolinsky and I. Kanter, Temporal association in asymmetric neural networks, Phys. Rev. Lett. 57, 2861 (1986).
  * Gómez-Nava  _et al._ [2022] L. Gómez-Nava, R. Bon, and F. Peruani, Intermittent collective motion in sheep results from alternating the role of leader and follower, Nature Physics 18, 1494 (2022).
  * MAY [1972] R. M. MAY, Will a large complex system be stable?, Nature 238, 413 (1972).
  * Yoshino  _et al._ [2007] Y. Yoshino, T. Galla, and K. Tokita, Statistical mechanics and stability of a model eco-system, Journal of Statistical Mechanics: Theory and Experiment 2007, P09003 (2007).
  * Diederich and Opper [1989] S. Diederich and M. Opper, Replicators with random interactions: A solvable model, Phys. Rev. A 39, 4333 (1989).
  * Fisher and Mehta [2014] C. K. Fisher and P. Mehta, The transition between the niche and neutral regimes in ecology, Proceedings of the National Academy of Sciences 111, 13111 (2014), https://www.pnas.org/doi/pdf/10.1073/pnas.1405637111 .
  * Dickens  _et al._ [2016] B. Dickens, C. K. Fisher, and P. Mehta, Analytically tractable model for community ecology with many species, Phys. Rev. E 94, 022423 (2016).
  * Marsland  _et al._ [2020a] R. Marsland, W. Cui, and P. Mehta, A minimal model for microbial biodiversity can reproduce experimentally observed ecological patterns, Scientific Reports 10, 3308 (2020a).
  * Advani  _et al._ [2018] M. Advani, G. Bunin, and P. Mehta, Statistical physics of community ecology: a cavity solution to macarthur’s consumer resource model, Journal of Statistical Mechanics: Theory and Experiment , 033406 (2018).
  * Cui  _et al._ [2021] W. Cui, R. Marsland III, and P. Mehta, Diverse communities behave like typical random ecosystems, Physical Review E 104, 034416 (2021).
  * Marsland  _et al._ [2020b] R. Marsland, W. Cui, and P. Mehta, The minimum environmental perturbation principle: A new perspective on niche theory, The American Naturalist 196, 291 (2020b).
  * Pearce  _et al._ [2020] M. T. Pearce, A. Agarwala, and D. S. Fisher, Stabilization of extensive fine-scale diversity by ecologically driven spatiotemporal chaos, Proceedings of the National Academy of Sciences 117, 14572 (2020), https://www.pnas.org/doi/pdf/10.1073/pnas.1915313117 .
  * Mahadevan  _et al._ [2023] A. Mahadevan, M. T. Pearce, and D. S. Fisher, Spatiotemporal ecological chaos enables gradual evolutionary diversification without niches or tradeoffs, eLife 12, e82734 (2023).
  * Altieri  _et al._ [2021] A. Altieri, F. Roy, C. Cammarota, and G. Biroli, Properties of equilibria and glassy phases of the random lotka-volterra model with demographic noise, Phys. Rev. Lett. 126, 258301 (2021).
  * Roy  _et al._ [2020] F. Roy, M. Barbier, G. Biroli, and G. Bunin, Complex interactions can create persistent fluctuations in high-diversity ecosystems, PLOS Computational Biology 16, 1 (2020).
  * Hu  _et al._ [2022] J. Hu, D. R. Amor, M. Barbier, G. Bunin, and J. Gore, Emergent phases of ecological diversity and dynamics mapped in microcosms, Science 378, 85 (2022).
  * Arnoulx de Pirey and Bunin [2023] T. Arnoulx de Pirey and G. Bunin, Aging by near-extinctions in many-variable interacting populations, Phys. Rev. Lett. 130, 098401 (2023).
  * de Pirey and Bunin [2023] T. A. de Pirey and G. Bunin, Many-species ecological fluctuations as a jump process from the brink of extinction (2023), arXiv:2306.13634 [q-bio.PE] .
  * MacArthur and Levins [1967] R. MacArthur and R. Levins, The limiting similarity, convergence, and divergence of coexisting species, The american naturalist 101, 377 (1967).
  * Macarthur [1970] R. Macarthur, Species packing and competitive equilibrium for many species*, POPULATION BIOLOGY 1, 11 (1970).
  * Chesson [1990] P. Chesson, Macarthur’s consumer-resource model, Theoretical Population Biology 37, 26 (1990).
  * Tilman [1982] D. Tilman,  _Resource competition and community structure_ (Princeton university press, 1982).
  * Chase and Leibold [2009] J. M. Chase and M. A. Leibold,  _Ecological niches: linking classical and contemporary approaches_ (University of Chicago Press, 2009).
  * Bunin [2017] G. Bunin, Ecological communities with lotka-volterra dynamics, Phys. Rev. E 95, 042414 (2017).
  * Nishimori [2001] H. Nishimori,  _Statistical Physics of Spin Glasses and Information Processing: An Introduction_, International series of monographs on physics (Oxford University Press, 2001).
  * Rackauckas and Nie [2017] C. Rackauckas and Q. Nie, DifferentialEquations.jl–a performant and feature-rich ecosystem for solving differential equations in Julia, Journal of Open Research Software 5 (2017).
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A fresh approach to numerical computing (2017).
  * Tsitouras [2011] C. Tsitouras, Runge–kutta pairs of order 5 (4) satisfying only the first column simplifying assumption, Computers & Mathematics with Applications 62, 770–775 (2011).
  * Ernst Hairer [1993] S. P. N. Ernst Hairer, Gerhard Wanner,  _Solving Ordinary Differential Equations I_ (Springer Berlin Heidelberg, 1993).
  * SciML [2022] J. SciML, Blackboxoptim.jl, GitHub (2022).
  * Dixit and Rackauckas [2023] V. K. Dixit and C. Rackauckas, Optimization.jl: A unified optimization package (2023).
  * Benettin  _et al._ [1980] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. part 1: Theory, Meccanica 15, 9 (1980).
  * Geist  _et al._ [1990] K. Geist, U. Parlitz, and W. Lauterborn, Comparison of Different Methods for Computing Lyapunov Exponents, Progress of Theoretical Physics 83, 875 (1990), https://academic.oup.com/ptp/article-pdf/83/5/875/5302061/83-5-875.pdf .
  * Datseris and Parlitz [2022] G. Datseris and U. Parlitz,  _Nonlinear Dynamics_ (Springer International Publishing, 2022).
  * Datseris [2018] G. Datseris, Dynamicalsystems.jl: A julia software library for chaos and nonlinear dynamics, Journal of Open Source Software 3, 598 (2018).
  * Skokos  _et al._ [2007] C. Skokos, T. Bountis, and C. Antonopoulos, Geometrical properties of local dynamics in hamiltonian systems: The generalized alignment index (gali) method, Physica D: Nonlinear Phenomena 231, 30 (2007).
  * Cui  _et al._ [2020] W. Cui, R. Marsland III, and P. Mehta, Effect of resource dynamics on species packing in diverse ecosystems, Physical Review Letters 125, 048101 (2020).
  * Marsland  _et al._ [2020c] R. Marsland, W. Cui, J. Goldford, and P. Mehta, The community simulator: A python package for microbial ecology, Plos one 15, e0230430 (2020c).
  * Dal Bello  _et al._ [2021] M. Dal Bello, H. Lee, A. Goyal, and J. Gore, Resource–diversity relationships in bacterial communities reflect the network structure of microbial metabolism, Nature Ecology & Evolution 5, 1424 (2021).
  * Marcus  _et al._ [2023] S. Marcus, A. M. Turner, and G. Bunin, Local and extensive fluctuations in sparsely-interacting ecological communities, arXiv preprint arXiv:2308.01828 (2023).
  * Feng  _et al._ [2023] Z. Feng, R. M. I. au2, J. W. Rocks, and P. Mehta, Emergent competition shapes the ecological properties of multi-trophic ecosystems (2023), arXiv:2303.02983 [q-bio.PE] .
  * Ebrahimi  _et al._ [2022] A. Ebrahimi, A. Goyal, and O. X. Cordero, Particle foraging strategies promote microbial diversity in marine environments, Elife 11, e73948 (2022).

###### Contents

  1. A Cavity calculation
    1. A.1 Setup
    2. A.2 Cavity solution
      1. A.2.1 Self-consistency equations for species populations
      2. A.2.2 Self-consistency equations for resource abundances
    3. A.3 ReLU function-transformed normal distributions
    4. A.4 Final self-consistency equations
    5. A.5 Infeasibility of the cavity solution
    6. A.6 Comparison of cavity and simulation results
  2. B Stability phase transition in the thermodynamic limit
  3. C Chaotic nature of dynamics
    1. C.1 Analysis of the Lyapunov exponents
    2. C.2 Analysis of the generalized alignment index (GALI)
  4. D Simulation and numerical methods
    1. D.1 Numerical integration
    2. D.2 Surviving/depletion cutoffs, numerical instability, and small immigration
    3. D.3 Numerically solving the cavity self-consistency equations and calculating phase boundaries
  5. E Finite size scaling analysis
    1. E.1 Modeling simulation outcomes
    2. E.2 Likelihood function
    3. E.3 Fitting the model
    4. E.4 Scaling of parameters with system size and discussion
  6. F Parameters in figures

##  Appendix A Cavity calculation

The objective of the cavity calculation is to find the steady-state behavior
of the aMCRM. In particular, we will find the distribution of the steady-state
abundances of the species and resources. The cavity method takes advantage of
the system’s self-averaging behavior. A system is said to have self-averaging
behavior if the distribution of properties of constituents are independent of
the exact realization of quenched disorder. This means that constituents’
properties can be treated as random variates drawn from a distribution that is
common to all systems with the same distribution of quenched disorder. A
consequence of this is that an average of some observable taken over
constituents given a fixed realization of quenched disorder is equal to the
average of the observable for one constituent taken over all realizations of
quenched disorder. The emergence of this self-averaging behavior is a
consequence of the central limit theorem and quenched disorder, as we will
see. In Fig. A5, we compare the distribution of steady-state resource and
species abundance for a single realization of quenched disorder to the
distribution of abundances predicted by the cavity method.

Figure A5:  Comparison of the empirical cumulative distribution functions
(eCDFs) of the steady-state abundances of species ($S=512$) and resources
($M=512$) for a single realization of quenched disorder and the distribution
predicted by the cavity calculation; simulations are performed in the stable
phase. See appendix F for simulation parameters.

###  A.1 Setup

We begin by introducing the constituent averages,

$\expectationvalue{R}\equiv\frac{1}{M}\sum_{\alpha=1}^{M}{\overline{R}_{\alpha}},\qquad\expectationvalue{N}\equiv\frac{1}{S}\sum_{i=1}^{S}{\overline{N}_{i}},$
(7)

where we have introduced the notation $\overline{X}$ to denote the steady-
state value of a quantity $X$. With these quantities defined, we can write the
steady-state aMCRM as,

$\displaystyle 0=\derivative{N_{i}}{t}$
$\displaystyle=\overline{N}_{i}\quantity[g+\frac{\sigma_{c}}{\sqrt{M}}\sum_{\alpha=1}^{M}d_{i\alpha}\overline{R}_{\alpha}-\sigma_{m}\delta
m_{i}],$ (8) $\displaystyle 0=\derivative{R_{\alpha}}{t}$
$\displaystyle=\overline{R}_{\alpha}\quantity[\kappa-\overline{R}_{\alpha}-\frac{\sigma_{e}}{\sqrt{M}}\sum_{i=1}^{S}N_{i}\left(\rho
d_{i\alpha}+\sqrt{1-\rho^{2}}x_{i\alpha}\right)+\sigma_{K}\delta K_{\alpha}],$
(9)

where

$\displaystyle g\equiv\mu_{c}\expectationvalue{R}-m,\qquad\kappa\equiv
K-\mu_{e}\gamma^{-1}\expectationvalue{N},\qquad\gamma\equiv\frac{M}{S}.$ (10)

###  A.2 Cavity solution

The cavity method begins by introducing a new species $i=0$ and a new resource
$\alpha=0$. We will later treat the new terms introduced by these variables as
small perturbations. With the new species and resource, we can write the
steady-state aMCRM for the existing species $i=1,\dots,S$ and resources
$\alpha=1,\dots,M$ as,

$\displaystyle 0=\derivative{N_{i}}{t}$
$\displaystyle=\overline{N}_{i}\quantity[\mu_{c}\expectationvalue{R}-\quantity(m-\frac{\sigma_{c}}{\sqrt{M}}d_{i0}R_{0})+\frac{\sigma_{c}}{\sqrt{M}}\sum_{\alpha=1}^{M}d_{i\alpha}\overline{R}_{\alpha}-\sigma_{m}\delta
m_{i}],$ (11) $\displaystyle 0=\derivative{R_{\alpha}}{t}$
$\displaystyle=\overline{R}_{\alpha}\left[\quantity(K-\frac{\sigma_{e}}{\sqrt{M}}N_{0}\quantity(\rho
d_{0\alpha}+\sqrt{1-\rho^{2}}x_{0\alpha}))-\mu_{e}\gamma^{-1}\expectationvalue{N}-\overline{R}_{\alpha}\right.$
$\displaystyle\left.\qquad-\frac{\sigma_{e}}{\sqrt{M}}\sum_{i=1}^{S}\overline{N}_{i}\quantity(\rho
d_{i\alpha}+\sqrt{1-\rho^{2}}x_{i\alpha})+\sigma_{K}\delta K_{\alpha}\right].$
(12)

For the new species and resource, the steady-state aMCRM says,

$\displaystyle 0=\derivative{N_{0}}{t}$
$\displaystyle=\overline{N}_{0}\quantity[g+\frac{\sigma_{c}}{\sqrt{M}}\sum_{\alpha=1}^{M}d_{0\alpha}\overline{R}_{\alpha}+M^{-1/2}\sigma_{c}d_{00}\overline{R}_{0}-\sigma_{m}\delta
m_{0}],$ (13) $\displaystyle 0=\derivative{R_{0}}{t}$
$\displaystyle=\overline{R}_{0}\left[\kappa-
R_{0}-\frac{\sigma_{e}}{\sqrt{M}}\sum_{i=1}^{S}\overline{N}_{j}\left(\rho
d_{i0}+\sqrt{1-\rho^{2}}x_{i0}\right)\right.$
$\displaystyle\qquad\qquad\left.-\frac{\sigma_{e}}{\sqrt{M}}\overline{N}_{0}\left(\rho
d_{00}+\sqrt{1-\rho^{2}}x_{00}\right)+\sigma_{K}\delta K_{\alpha}\right[.$
(14)

Next, we will analyze the perturbed system relative to the unperturbed system.
A quantity with $\setminus 0$ represents the value before adding the new
variables to the system. Looking to Eqs. 11 and 12, we see that the presence
of the new species and resource effectively perturbs the model parameters as,

$\displaystyle m_{i}\to
m_{i}-\frac{\sigma_{c}}{\sqrt{M}}d_{i0}\overline{R}_{0},\qquad K_{\alpha}\to
K_{\alpha}-\frac{\sigma_{e}}{\sqrt{M}}\overline{N}_{0}\left(\rho
d_{0\alpha}+\sqrt{1-\rho^{2}}x_{0\alpha}\right).$ (15)

In the thermodynamic limit where $M$ and $S$ are large, we model the
perturbation using linear response:

$\displaystyle\overline{N}_{i}$ $\displaystyle=\overline{N}_{i\setminus
0}-\frac{\sigma_{e}}{\sqrt{M}}\sum_{\beta=1}^{M}\chi_{i\beta}^{(N)}\left(\rho
d_{0\beta}+\sqrt{1-\rho^{2}}x_{0\beta}\right)\overline{N}_{0}-\frac{\sigma_{c}}{\sqrt{M}}\sum_{j=1}^{S}\nu_{ij}^{(N)}d_{j0}\overline{R}_{0},$
(16) $\displaystyle\overline{R}_{\alpha}$
$\displaystyle=\overline{R}_{\alpha\setminus
0}-\frac{\sigma_{e}}{\sqrt{M}}\sum_{\beta=1}^{M}\chi_{\alpha\beta}^{(R)}\left(\rho
d_{0\beta}+\sqrt{1-\rho^{2}}x_{0\beta}\right)\overline{N}_{0}-\frac{\sigma_{c}}{\sqrt{M}}\sum_{j=1}^{S}\nu_{\alpha
j}^{(R)}d_{j0}\overline{R}_{0},$ (17)

where we have defined the susceptibility matrices:

$\displaystyle\chi_{i\beta}^{(N)}\equiv\frac{\partial\overline{N}_{i}}{\partial
K_{\beta}},$
$\displaystyle\chi_{\alpha\beta}^{(R)}\equiv\frac{\partial\overline{R}_{\alpha}}{\partial
K_{\beta}},$ (18)
$\displaystyle\nu_{ij}^{(N)}\equiv\frac{\partial\overline{N}_{i}}{\partial
m_{j}},$ $\displaystyle\nu_{\alpha
j}^{(R)}\equiv\frac{\partial\overline{R}_{\alpha}}{\partial m_{j}}.$ (19)

####  A.2.1 Self-consistency equations for species populations

With the perturbation introduced and susceptibilities defined, we will exploit
the linear response approximation to derive the self-consistency equations for
the species populations. Substituting Eqs. 16, 17 into the aMCRM equation for
the additional species [Eq. (13)], we obtain,

$\displaystyle 0$
$\displaystyle=\overline{N}_{0}\left[g+\frac{\sigma_{c}}{\sqrt{M}}\sum_{\alpha=1}^{M}d_{0\alpha}\overline{R}_{\alpha\setminus
0}-\frac{\sigma_{c}\sigma_{e}}{M}\sum_{\alpha,\beta=1}^{M}\chi_{\alpha\beta}^{(R)}d_{0\alpha}\left(\rho
d_{0\beta}+\sqrt{1-\rho^{2}}x_{0\beta}\right)\overline{N}_{0}\right.$
$\displaystyle\qquad\qquad\left.-\frac{\sigma_{c}^{2}}{M}\sum_{\alpha=1}^{M}\sum_{j=1}^{S}\nu_{\alpha
j}^{(R)}d_{0\alpha}d_{j0}\overline{N}_{0}+\frac{\sigma_{c}}{\sqrt{M}}d_{00}\overline{R}_{0}-\sigma_{m}\delta
m_{0}\right].$ (20)

Taking the mean of the third term with respect to the new matrix elements
$d_{0\alpha}$ and $x_{0\alpha}$, we obtain

$\displaystyle\expectationvalue{\frac{\sigma_{c}\sigma_{e}}{M}\sum_{\alpha,\beta=1}^{M}\chi_{\alpha\beta}^{(R)}d_{0\alpha}\left(\rho
d_{0\beta}+\sqrt{1-\rho^{2}}x_{0\beta}\right)\overline{N}_{0}}$
$\displaystyle\qquad=\overline{N}_{0}\frac{\sigma_{c}\sigma_{e}}{M}\sum_{\alpha,\beta=1}^{M}{\chi_{\alpha\beta}^{(R)}}\left(\rho\expectationvalue{d_{0\alpha}d_{0\beta}}+\expectationvalue{d_{0\alpha}}\expectationvalue{x_{0\beta}}\sqrt{1-\rho^{2}}\right)$
$\displaystyle\qquad=\overline{N}_{0}\frac{\sigma_{c}\sigma_{e}}{M}\sum_{\alpha,\beta=1}^{M}{\chi_{\alpha\beta}^{(R)}}\left(\rho\delta_{\alpha\beta}+0\times
0\sqrt{1-\rho^{2}}\right)=\overline{N}_{0}\sigma_{c}\sigma_{e}\rho\chi,$ (21)

where we have used the fact that $d_{0\alpha}$ and $x_{0\beta}$ are
independent and have defined

$\displaystyle\chi\equiv\frac{1}{M}\sum_{\alpha=1}^{M}{\chi_{\alpha\alpha}^{(R)}},$
(22)

to be the trace of the $\overline{R}_{\alpha}\leftarrow K_{\beta}$
susceptibility matrix divided by the number of resources. The variance of this
term is $O(M^{-1})$, which can be verified by expanding its second moment. The
mean of the fourth term is zero because $d_{0\alpha}$ and $d_{j0}$ are
uncorrelated when $\alpha=1,\dots,M,\;j=1,\dots,S$ and its variance is
$O(M^{-1})$. Discarding terms of order $O(M^{-1/2})$ and higher, we obtain,

$\displaystyle
0=\overline{N}_{0}\left(g-\sigma_{c}\sigma_{e}\rho\chi\overline{N}_{0}+\frac{\sigma_{c}}{\sqrt{M}}\sum_{\alpha=1}^{M}d_{0\alpha}\overline{R}_{\alpha\setminus
0}-\sigma_{m}\delta m_{0}\right)+O(M^{-1/2}).$ (23)

The last two terms above are a sum of many independent random variables, so,
by the central limit theorem, we can model these terms as a sum of normal
random variables. The mean of these terms is,

$\displaystyle\expectationvalue{\frac{\sigma_{c}}{\sqrt{M}}\sum_{\alpha=1}^{M}d_{0\alpha}\overline{R}_{\alpha\setminus
0}-\sigma_{m}\delta m_{0}}$
$\displaystyle=\frac{\sigma_{c}}{\sqrt{M}}\sum_{\alpha=1}^{M}\expectationvalue{d_{0\alpha}}\overline{R}_{\alpha\setminus
0}-\sigma_{m}\expectationvalue{\delta m_{0}}$
$\displaystyle=\frac{\sigma_{c}}{\sqrt{M}}\sum_{\alpha=1}^{M}0\times{\overline{R}_{\alpha\setminus
0}}-\sigma_{m}\times 0=0.$ (24)

The variance of these terms is,

$\displaystyle\sigma_{g}^{2}$
$\displaystyle\equiv\mathrm{Var}\left[\frac{\sigma_{c}}{\sqrt{M}}\sum_{\alpha=1}^{M}d_{0\alpha}\overline{R}_{\alpha\setminus
0}-\sigma_{m}\delta
m_{0}\right]=\frac{\sigma_{c}^{2}}{M}\sum_{\alpha=1}^{M}{\overline{R}^{2}_{\alpha\setminus
0}}\mathrm{Var}\left[d_{0\alpha}\right]+\sigma_{m}^{2}\mathrm{Var}\left[\delta
m_{0}\right]$
$\displaystyle=\frac{\sigma_{c}^{2}}{M}\sum_{\alpha=1}^{M}{\overline{R}_{\alpha\setminus
0}^{2}}+\sigma_{m}^{2}=\sigma_{c}^{2}q_{R}+\sigma_{m}^{2},$ (25)

where we have defined $q_{R}\equiv
M^{-1}\sum_{\alpha=1}^{M}{\overline{R}_{\alpha\setminus 0}^{2}}.$ Due to the
self-averaging nature of the system and that the size of the perturbation is
of order $O(M^{-1/2})$, we have
$q_{R}\approx\expectationvalue*{\overline{R}_{0}^{2}}$. Let
$Z_{N}\sim\mathcal{N}(0,1)$ be a unit normal random variable. The large-$M$
limit approximate steady-state condition for the perturbing species becomes,

$\displaystyle
0=\overline{N}_{0}\left(g-\sigma_{c}\sigma_{e}\rho\chi\overline{N}_{0}+\sigma_{g}Z_{N}\right).$
(26)

Solving for $\overline{N}_{0}$ and discarding invadable solutions, we obtain,

$\displaystyle\boxed{\overline{N}_{0}=\max\left\\{0,\frac{g+\sigma_{g}Z_{N}}{\sigma_{c}\sigma_{e}\rho\chi}\right\\}}.$
(27)

Additionally, observe that if $\chi\to 0$, the cavity solution diverges; this
will be discussed in further detail in Appendix A.5.

####  A.2.2 Self-consistency equations for resource abundances

Now, we repeat this process to find a self-consistency equation for the
resources. We substitute the linear response approximation for species into
the aMCRM steady-state equation for the additional resource:

$\displaystyle 0$
$\displaystyle=\overline{R}_{0}\left[\kappa-\overline{R}_{0}-\frac{\sigma_{e}}{\sqrt{M}}\sum_{i=1}^{S}\overline{N}_{i\setminus
0}\left(\rho d_{i0}+\sqrt{1-\rho^{2}}x_{i0}\right)\right.$ (28)
$\displaystyle\qquad\left.+\frac{\sigma_{e}^{2}}{M}\sum_{i=1}^{S}\sum_{\beta=1}^{M}\chi_{i\beta}^{(N)}\left(\rho
d_{i0}+\sqrt{1-\rho^{2}}x_{i0}\right)\left(\rho
d_{0\beta}+\sqrt{1-\rho^{2}}x_{0\beta}\right)\overline{N}_{0}\right.$
$\displaystyle\qquad\left.+\frac{\sigma_{e}\sigma_{c}}{M}\sum_{i,j=1}^{S}\nu_{ij}^{(N)}\left(\rho
d_{i0}+\sqrt{1-\rho^{2}}x_{i0}\right)d_{j0}\overline{R}_{0}-\frac{\sigma_{e}}{\sqrt{M}}\overline{N}_{0}\left(\rho
d_{00}+\sqrt{1-\rho^{2}}x_{00}\right)+\sigma_{K}\delta K_{\alpha}\right].$

Observe that the fourth term (involving $\chi_{i\beta}^{(N)}$) has zero mean
and variance of order $O(M^{-1})$; this can be seen by recalling that
$d_{i0},d_{0\beta},x_{i0},x_{0\beta}$ are all independent for $i,\beta\geq 1$.
Similarly, we see that the variance of the fifth term (involving
$\nu_{ij}^{(N)}$) is of order $O(M^{-1})$. We will ignore fluctuations of
order $O(M^{-{1/2}})$. The mean of the fifth term is,

$\displaystyle{\frac{\sigma_{e}\sigma_{c}}{M}\sum_{i,j=1}^{S}\nu_{ij}^{(N)}\left(\rho
d_{i0}+\sqrt{1-\rho^{2}}x_{i0}\right)d_{j0}\overline{R}_{0}}=\overline{R}_{0}\frac{\sigma_{e}\sigma_{c}}{M}\sum_{i,j=1}^{S}{\nu_{ij}^{(N)}}\left(\rho\expectationvalue{d_{i0}d_{j0}}+\sqrt{1-\rho^{2}}\expectationvalue{x_{i0}}\expectationvalue{d_{j0}}\right)$
$\displaystyle\qquad=\overline{R}_{0}\sigma_{e}\sigma_{c}\frac{S}{M}\frac{1}{S}\sum_{i,j=1}^{S}{\nu_{ij}^{(N)}}\left(\rho\delta_{ij}+\sqrt{1-\rho^{2}}\times
0\times 0\right)=\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu\overline{R}_{0},$ (29)

where we have defined,

$\displaystyle\nu\equiv\frac{1}{S}\sum_{i=1}^{S}{\nu_{ii}^{(N)}},$ (30)

which is the trace of the $\overline{N}_{i}\leftarrow m_{j}$ susceptibility
matrix divided by $S$. Discarding terms of order $O(M^{-1/2})$, we obtain,

$\displaystyle
0=\overline{R}_{0}\left(\kappa-\overline{R}_{0}-\frac{\sigma_{e}}{\sqrt{M}}\sum_{i=1}^{S}\overline{N}_{i\setminus
0}\left(\rho
d_{i0}+\sqrt{1-\rho^{2}}x_{i0}\right)+\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu\overline{R}_{0}+\sigma_{K}\delta
K_{0}\right)+O(M^{-1/2}).$ (31)

Now, observe the third and last terms are a sum of many independent random
variables, meaning we can apply the central limit theorem and model the sum of
these terms as a normal random variable. The mean of these terms is,

$\displaystyle\expectationvalue{\sigma_{K}\delta
K_{0}-\frac{\sigma_{e}}{\sqrt{M}}\sum_{i=1}^{S}\overline{N}_{i\setminus
0}\left(\rho
d_{i0}+\sqrt{1-\rho^{2}}x_{i0}\right)}=\sigma_{K}\expectationvalue{\delta
K_{0}}-\frac{\sigma_{e}}{\sqrt{M}}\sum_{i=1}^{S}{\overline{N}_{i\setminus
0}}\left(\rho\expectationvalue{d_{i0}}+\sqrt{1-\rho^{2}}\expectationvalue{x_{i0}}\right)$
$\displaystyle\qquad=\sigma_{K}\times
0-\frac{\sigma_{e}}{\sqrt{M}}\sum_{i=1}^{S}{\overline{N}_{i\setminus
0}}\left(\rho\times 0+\sqrt{1-\rho^{2}}\times 0\right)=0.$ (32)

The variance is,

$\displaystyle\sigma_{\kappa}^{2}$
$\displaystyle\equiv\mathrm{Var}\left[\sigma_{K}\delta
K_{0}-\frac{\sigma_{e}}{\sqrt{M}}\sum_{i=1}^{S}\overline{N}_{i\setminus
0}\left(\rho
d_{i0}+\sqrt{1-\rho^{2}}x_{i0}\right)\right]=\sigma_{K}^{2}\mathrm{Var}\left[\delta
K_{0}\right]+\frac{\sigma_{e}^{2}}{M}\sum_{i=1}^{S}\mathrm{Var}\left[\overline{N}_{i\setminus
0}\left(\rho d_{i0}+\sqrt{1-\rho^{2}}x_{i0}\right)\right]$
$\displaystyle=\sigma_{K}^{2}+\frac{\sigma_{e}^{2}}{M}\sum_{i=1}^{S}{\overline{N}_{i\setminus
0}^{2}}\left(\rho^{2}\mathrm{Var}\left[d_{i0}\right]+(1-\rho^{2})\mathrm{Var}\left[x_{i0}\right]\right)=\sigma_{K}^{2}+{\sigma_{e}^{2}}\frac{S}{M}\frac{1}{S}\sum_{i=1}^{S}{\overline{N}_{i\setminus
0}^{2}}=\sigma_{K}^{2}+\gamma^{-1}\sigma_{e}^{2}q_{N},$ (33)

where, $q_{N}\equiv\frac{1}{S}\sum_{i=1}^{S}{N_{i\setminus 0}^{2}}$ which is
approximately equal to $\expectationvalue*{\overline{N}_{0}^{2}}$. The
approximate steady-state condition for the added resource then becomes,

$\displaystyle
0=\overline{R}_{0}\left(\kappa-\overline{R}_{0}+\sigma_{\kappa}Z_{R}+\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu\overline{R}_{0}\right),$
(34)

where $Z_{R}\sim\mathcal{N}(0,1)$ is a standard normal random variable.
Solving for $\overline{R}_{0}$ and discarding invadable solutions gives,

$\displaystyle\boxed{\overline{R}_{0}=\max\left\\{0,\frac{\kappa+\sigma_{\kappa}Z_{R}}{1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu}\right\\}}.$
(35)

###  A.3 ReLU function-transformed normal distributions

In these computations, we regularly work with normal distributions that are
transformed by the ‘ReLU’ function:
$\mathrm{ReLU}(x)=\max\\{0,x\\}=x\Theta(x)$. If $Z$ is a standard normal
random variable, the PDF of $\mathrm{ReLU}(\sigma Z+\mu)$ is,

$\displaystyle p_{\mathrm{ReLU}(\sigma
Z+\mu)}(z)=\delta(z)\Phi(-\mu/\sigma)+\frac{1}{\sqrt{2\pi}\sigma}e^{-(z-\mu)^{2}/2\sigma^{2}}\Theta(z),$
(36)

where,

$\displaystyle\Phi(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}\mathrm{d}z\,e^{-z^{2}/2}=\frac{1}{2}\left(1+\mathrm{erf}(x/\sqrt{2})\right),$
(37)

is the standard normal CDF. The $j$th ($j\geq 1$) moment is then,

$\displaystyle W_{j}(\mu,\sigma)$ $\displaystyle=\langle\mathrm{ReLU}(\sigma
Z+\mu)^{j}\rangle=0+\frac{1}{\sqrt{2\pi}\sigma}\int_{0}^{\infty}\mathrm{d}z\,z^{j}e^{-(z-\mu)^{2}/2\sigma^{2}}$
(38)
$\displaystyle=\sigma^{j}\int_{-\mu/\sigma}^{\infty}\frac{\mathrm{d}z}{\sqrt{2\pi}}e^{-z^{2}/2}(z+\mu/\sigma)^{j},$
$\displaystyle=\frac{2^{-3/2}}{\sqrt{\pi}}(\sqrt{2}\sigma)^{j}\left[j\frac{\mu}{\sigma}\Gamma\left(\frac{j}{2}\right)\,_{1}F_{1}\left(\frac{1-j}{2};\frac{3}{2};-\frac{\mu^{2}}{2\sigma^{2}}\right)+\sqrt{2}\Gamma\left(\frac{j+1}{2}\right)\,_{1}F_{1}\left(-\frac{j}{2};\frac{1}{2};-\frac{\mu^{2}}{2\sigma^{2}}\right)\right],$

where $\,{}_{1}F_{1}$ is the confluent hypergeometric function of the first
kind, and $\Gamma$ is the gamma function. Observe that
$W_{j}(\mu/\alpha,\sigma/\alpha)=\alpha^{-j}W_{j}(\mu,\sigma)$. Additionally,

$\displaystyle W_{0}(x,1)$ $\displaystyle=1,$ (39) $\displaystyle W_{1}(x,1)$
$\displaystyle=\frac{1}{\sqrt{2\pi}}e^{-x^{2}/2}+x\Phi(x),$ (40)
$\displaystyle W_{2}(x,1)$
$\displaystyle=\frac{1}{\sqrt{2\pi}}xe^{-x^{2}/2}+(1+x^{2})\Phi(x).$ (41)

It follows from integration by parts,

$\displaystyle W_{2}(x,1)=\Phi(x)+xW_{1}(x,1).$ (42)

Additionally, for a random variable $\Theta(\sigma Z+\mu)$, the PDF is,

$\displaystyle p_{\Theta(\sigma
Z+\mu)}(z)=\frac{1}{2}\quantity[1+\mathrm{erf}\left(\frac{\mu}{\sigma\sqrt{2}}\right)]\delta(z-1)+\frac{1}{2}\mathrm{erfc}\left(\frac{\mu}{\sigma\sqrt{2}}\right)\delta(z),$
(43)

so the $j$th moment ($j\geq 1$) is,

$\displaystyle\langle\Theta(\sigma
Z+\mu)^{j}\rangle=0+\frac{1}{2}\quantity[1+\mathrm{erf}\left(\frac{\mu}{\sigma\sqrt{2}}\right)]1^{j}=\frac{1}{2}\quantity[1+\mathrm{erf}\left(\frac{\mu}{\sigma\sqrt{2}}\right)]=\Phi(\mu/\sigma).$
(44)

###  A.4 Final self-consistency equations

Some essential quantities of interest are the expected fraction of surviving
species $\phi_{N}$ and fraction of non-depleted resources $\phi_{R}$. These
quantities are computed using the moments calculated in section A.3 and Eqs.
26, 34:

$\displaystyle\boxed{\phi_{N}=\expectationvalue{\Theta(\overline{N}_{0})}=\Phi(\Delta_{g})},$
(45)
$\displaystyle\boxed{\phi_{R}=\expectationvalue{\Theta(\overline{R}_{0})}=\Phi(\Delta_{\kappa})},$
(46)

where $\Delta_{g}=g/\sigma_{g}$ and $\Delta_{\kappa}=\kappa/\sigma_{\kappa}$
and $\Theta$ is the Heaviside step function with the convention $\Theta(0)=0$.
Next, we can differentiate our expressions for $\overline{N}_{0}$ and
$\overline{R}_{0}$ to find,

$\displaystyle\frac{\partial\overline{N}_{0}}{\partial m}$
$\displaystyle=\frac{\partial}{\partial
m}\frac{g+\sigma_{g}Z_{N}}{\sigma_{c}\sigma_{e}\rho\chi}\Theta\left(\overline{N}_{0}\right)=-\frac{1}{\sigma_{c}\sigma_{e}\rho\chi}\Theta(\overline{N}_{0})+\left[\text{$\overline{N}_{0}\delta(\overline{N}_{0})$-term}\right]$
$\displaystyle\implies\expectationvalue{\frac{\partial\overline{N}_{0}}{\partial
m}}$
$\displaystyle=\boxed{\nu=-\frac{\phi_{N}}{\sigma_{c}\sigma_{e}\rho\chi}}$
(47) $\displaystyle\frac{\partial\overline{R}_{0}}{\partial K}$
$\displaystyle=\frac{\partial}{\partial
K}\frac{\kappa+\sigma_{\kappa}Z_{R}}{1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu}\Theta(\overline{R}_{0})=\frac{1}{1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu}\Theta(\overline{R}_{0})+\left[\text{$\overline{R}_{0}\delta(\overline{R}_{0})$-term}\right]$
$\displaystyle\implies\expectationvalue{\frac{\partial\overline{R}_{0}}{\partial
K}}$
$\displaystyle=\boxed{\chi=\frac{\phi_{R}}{1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu}}$
(48)

We can solve these two equations for $\chi,\nu$ to obtain the relations,

$\displaystyle\nu=\frac{1}{\rho\sigma_{c}\sigma_{e}\gamma^{-1}(1-\gamma\phi_{R}/\phi_{N})},\qquad\chi=\phi_{R}-\gamma^{-1}\phi_{N}.$
(49)

Next, we use Eqs. 26 and 34 and invoke our assumption that the system self-
averages to find,

$\displaystyle\boxed{\expectationvalue{N}=\expectationvalue{\overline{N}_{0}}=\frac{\sigma_{g}}{\sigma_{c}\sigma_{e}\rho\chi}W_{1}(\Delta_{g},1)=\frac{\sigma_{g}}{\sigma_{c}\sigma_{e}\rho\chi}\left(\frac{e^{-\Delta_{g}^{2}/2}}{\sqrt{2\pi}}+\Delta_{g}\Phi(\Delta_{g})\right)},$
(50)
$\displaystyle\boxed{\expectationvalue{R}=\expectationvalue{\overline{R}_{0}}=\frac{\sigma_{\kappa}}{1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu}W_{1}(\Delta_{\kappa},1)=\frac{\sigma_{\kappa}}{1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu}\left(\frac{e^{-\Delta_{\kappa}^{2}/2}}{\sqrt{2\pi}}+\Delta_{\kappa}\Phi(\Delta_{\kappa})\right)},$
(51)
$\displaystyle\boxed{q_{N}=\expectationvalue{\overline{N}_{0}^{2}}=\left(\frac{\sigma_{g}}{\sigma_{c}\sigma_{e}\rho\chi}\right)^{2}W_{2}(\Delta_{g},1)=\left(\frac{\sigma_{g}}{\sigma_{c}\sigma_{e}\rho\chi}\right)^{2}\left(\frac{\Delta_{g}e^{-\Delta_{g}^{2}/2}}{\sqrt{2\pi}}+(1+\Delta_{g}^{2})\Phi(\Delta_{g})\right)},$
(52)
$\displaystyle\boxed{q_{R}=\expectationvalue{\overline{R}_{0}^{2}}=\left(\frac{\sigma_{\kappa}}{1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu}\right)^{2}W_{2}(\Delta_{\kappa},1)=\left(\frac{\sigma_{\kappa}}{1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu}\right)^{2}\left(\frac{\Delta_{\kappa}e^{-\Delta_{\kappa}^{2}/2}}{\sqrt{2\pi}}+(1+\Delta_{\kappa}^{2})\Phi(\Delta_{\kappa})\right)}.$
(53)

The equations 45, 46, 47, 48, 50, 51, 52, 53 constitute the cavity self-
consistency equations for the aMCRM model. They are eight independent
nonlinear equations to solve for eight variables:
$\phi_{N},\phi_{R},\chi,\nu,\expectationvalue{N},\expectationvalue{R},q_{N},q_{R}$.
These equations are solved numerically using nonlinear least squares as
discussed in Appendix D.3.

Figure A6:  Comparison of cavity results to numerical simulations for various
$\rho$ and fixed $\sigma_{c}=\sigma_{e}=3.5$, corresponding to the gray dashed
line in Fig. 3. (a) Fractions of surviving species, $\phi_{N}$, and non-
depleted resources, $\phi_{R}$. (b) Average populations of species, $\langle
N\rangle$, and abundances of resources, $\langle R\rangle$. (c) Standard
deviations of species populations, $\sqrt{\langle N^{2}\rangle-\langle
N\rangle^{2}}$, and resource abundances, $\sqrt{\langle R^{2}\rangle-\langle
R\rangle^{2}}$. Points are means of numerical simulations, lines are cavity
results. Error bars are standard deviations of numerical simulations; some
error bars are smaller than the points. A dashed line is shown at
$\rho^{\star}=0.72$, the critical value of $\rho$ at which the system
transitions between the stable and dynamic phases. For each value of $\rho$,
$64$ simulations are performed with $M,S=512$ resources and species.  Figure
A7:  Cavity susceptibilities for the same parameters as in Fig. A6. The
susceptibilities $\chi$, $\nu$ as defined in lines 22, 30 are shown on the
left and right, respectively.

Figure A8:  Heatmap comparison of cavity predictions and simulation results.
Results from 32 simulations for various values of $\sigma$ and $\rho$ are
displayed in the left column of each pane and the cavity predictions are
displayed in the right column. The instability boundary is shown as a solid
black line while the infeasibility boundary is shown as a dashed black line.
For the cavity prediction plots, the parameter regions for which no cavity
solution exists have a diamond pattern overlaid (c.f., Fig. D16 and Appendix
D.3). The parameters used are those used in Fig. 3(a).  Figure A9:  Heatmap of
cavity susceptibilities $\chi$, $\nu$ as defined in lines 22, 30 for various
values of $\sigma$ and $\rho$, as in Fig. A8.

###  A.5 Infeasibility of the cavity solution

Next, we will investigate when there may exist a solution to the cavity self-
consistency equations. Observe that Eq. (47) implies that as $\chi\to 0$,
$\nu$ diverges; Eqs. 27, 35 indicate that $\overline{N}_{0}$ and
$\overline{R}_{0}$ become singular. This is when, numerically, there fails to
exist a solution to the cavity self-consistency equations. From Eq. (48), we
see that $\chi\to 0$ when,

$\displaystyle\phi_{R}^{\text{CSI}}-\gamma^{-1}\phi_{N}^{\text{CSI}}=0\Longleftrightarrow\frac{\phi_{R}^{\text{CSI}}}{\phi_{N}^{\text{CSI}}}\frac{M}{S}=1$
(54) $\displaystyle\Downarrow$ $\displaystyle\boxed{\text{cavity solution
infeasibility boundary:}\qquad\text{\\# of non-depleted resources}=\text{\\#
of surviving species}.}$ (55)

The cavity solution no longer has a solution when the ecosystem approaches the
bound set by the principle of competitive exclusion. According to this
calculation, $\nu$ diverging here means that if a species becomes marginally
more fit when all niches are packed, it will disrupt the ecosystem. This
relation defines a ninth equation for eight variables, so it determines a co-
dimension-one boundary in the space of parameters. This boundary can be
calculated numerically using nonlinear least squares as discussed in Appendix
D.3. When solving the cavity self-consistency equations numerically, we will
see that there is a region of parameter space beyond this boundary where the
least squares objectives are large compared to machine error (see Fig. D16 in
Appendix D.3).

###  A.6 Comparison of cavity and simulation results

Figure A10:  Heatmap of the fraction of simulations at steady state when the
sampling distributions are uniform. Parameters are the same as in Fig. 3(a).
Because the calculations are performed in the thermodynamic limit, the phase
boundary is agnostic to the choice of sampling distribution. Further details
about sampling distribution and parameters are in appendix F.

For higher reciprocity levels $\rho$ when the system is in the stable phase,
the analytic solution produced from the cavity method matches the simulation
results remarkably well. In Fig. A6, we compare the cavity results to
numerical simulations for various $0\leq\rho\leq 1$ and fixed
$\sigma_{c}=\sigma_{e}=3.5$, corresponding to the gray dashed line in Fig.
3(a). For all predicted quantities, the chi-squared statistic is near the
number of degrees of freedom, indicating that the cavity solution is a good
fit to the simulation results. The cavity susceptibilities $\chi,\nu$ for this
slice are shown in Fig. A7. In this figure, we see that the infeasibility
boundary occurs when $\chi\to 0$ and $\nu\to-\infty$.

In Fig. A8, we compare the cavity results to numerical simulations for the
various values $\rho$ and $\sigma_{c}=\sigma_{e}\equiv\sigma$ shown on the
grid in Fig. 3(a). Again, the analytic solutions produced from the cavity
method match the simulation results remarkably well in the stable phase.
Cavity susceptibilities $\chi,\nu$ for this grid are shown in Fig. A9.

##  Appendix B Stability phase transition in the thermodynamic limit

In nonlinear dynamics, a common approach to assess the appearance of
instability is to begin with an assumption of stability and see what breaks
down when the system is perturbed. In the case of the aMCRM, when there is
stability, the results from the cavity solution will be valid, so we will use
these results to determine the critical threshold for instability.

We consider perturbing all surviving species and resources,
$\overline{N}_{i}^{+}$, $\overline{R}_{i}^{+}$, by small amounts,
$\varepsilon\eta_{i}^{(N)}$, $\varepsilon\eta_{\alpha}^{(R)}$, respectively,
where $\eta_{i}^{(N)}$ and $\eta_{\alpha}^{(R)}$ are independent random
variables all with mean zero and variance one [26]. By considering
perturbations on surviving species and resources and using cavity results
(c.f., Eqs. 35, 27), we are assuming that we are working with an uninvadable
steady state. From Eqs. 23, 31, for surviving species and non-depleted
resources,

$\displaystyle\overline{N}_{0}^{+}=\frac{1}{\sigma_{c}\sigma_{e}\rho\chi}\quantity[g+\frac{\sigma_{c}}{\sqrt{M}}\sum_{\alpha:\overline{R}_{\alpha}>0}d_{0\alpha}\overline{R}_{\alpha\setminus
0}^{+}-\sigma_{m}\delta m_{0}],$ (56)
$\displaystyle\overline{R}_{0}^{+}=\frac{1}{1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu}\quantity[\kappa-\frac{\sigma_{e}}{\sqrt{M}}\sum_{i:\overline{N}_{i}>0}\overline{N}_{i\setminus
0}^{+}\left(\rho d_{i0}+\sqrt{1-\rho^{2}}x_{i0}\right)+\sigma_{K}\delta
K_{0}].$ (57)

Applying the perturbation gives,

$\displaystyle\overline{N}_{0}^{+}$
$\displaystyle=\frac{1}{\sigma_{c}\sigma_{e}\rho\chi}\quantity[g+\frac{\sigma_{c}}{\sqrt{M}}\sum_{\alpha:\overline{R}_{\alpha}>0}d_{0\alpha}\left(\overline{R}_{\alpha\setminus
0}^{+}+\varepsilon\eta_{\alpha}^{(R)}\right)-\sigma_{m}\delta m_{0}],$ (58)
$\displaystyle\overline{R}_{0}^{+}$
$\displaystyle=\frac{1}{1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu}\quantity[\kappa-\frac{\sigma_{e}}{\sqrt{M}}\sum_{i:\overline{N}_{i}>0}\left(\overline{N}_{i\setminus
0}^{+}+\varepsilon\eta_{i}^{(N)}\right)\left(\rho
d_{i0}+\sqrt{1-\rho^{2}}x_{i0}\right)+\sigma_{K}\delta K_{0}].$ (59)

Differentiating with respect to $\varepsilon$ yields,

$\displaystyle\derivative{\overline{N}_{0}^{+}}{\varepsilon}$
$\displaystyle=\frac{1}{\sigma_{e}\rho\chi\sqrt{M}}\sum_{\alpha:\overline{R}_{\alpha}>0}d_{0\alpha}\left(\derivative{\overline{R}_{\alpha\setminus
0}^{+}}{\varepsilon}+\eta_{\alpha}^{(R)}\right),$ (60)
$\displaystyle\derivative{\overline{R}_{0}^{+}}{\varepsilon}$
$\displaystyle=-\frac{{\sigma_{e}}/{\sqrt{M}}}{1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu}\sum_{i:\overline{N}_{i}>0}\left(\derivative{\overline{N}_{i\setminus
0}^{+}}{\varepsilon}+\eta_{i}^{(N)}\right)\left(\rho
d_{i0}+\sqrt{1-\rho^{2}}x_{i0}\right).$ (61)

Because $\overline{R}_{\alpha\setminus 0}^{+}$, $d_{0\alpha}$,
$\eta_{\alpha}^{(R)}$ are all independent, $\langle
d\overline{R}_{0}^{+}/\mathrm{d}\varepsilon\rangle=0$; similarly, $\langle
d\overline{N}_{0}^{+}/\mathrm{d}\varepsilon\rangle=0$. The first moment of the
response to the perturbation does not give useful information about the
stability; therefore, we turn to the second moment. First, we square the above
quantities:

$\displaystyle\left[\derivative{\overline{N}_{0}^{+}}{\varepsilon}\right]^{2}$
$\displaystyle=\frac{1/M}{(\sigma_{e}\rho\chi)^{2}}\sum_{\alpha,\beta:\overline{R}_{\alpha}>0,\overline{R}_{\beta}>0}d_{0\alpha}d_{0\beta}\left(\derivative{\overline{R}_{\alpha\setminus
0}^{+}}{\varepsilon}+\eta_{\alpha}^{(R)}\right)\left(\derivative{\overline{R}_{\beta\setminus
0}^{+}}{\varepsilon}+\eta_{\beta}^{(R)}\right),$ (62)
$\displaystyle\left[\derivative{\overline{R}_{0}^{+}}{\varepsilon}\right]^{2}$
$\displaystyle=\frac{\sigma_{e}^{2}/M}{(1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu)^{2}}\sum_{i,j:\overline{N}_{i}>0,\overline{N}_{j}>0}\left(\derivative{\overline{N}_{i\setminus
0}^{+}}{\varepsilon}+\eta_{i}^{(N)}\right)\left(\derivative{\overline{N}_{j\setminus
0}^{+}}{\varepsilon}+\eta_{j}^{(N)}\right)$ (63) $\displaystyle\hskip
142.26378pt\times\left(\rho d_{i0}+\sqrt{1-\rho^{2}}x_{i0}\right)\left(\rho
d_{j0}+\sqrt{1-\rho^{2}}x_{j0}\right).$

Averaging over all sources of randomness and using
$\expectationvalue{\left[\derivative{N_{0}}{\varepsilon}\right]^{2}}=S^{-1}\sum_{i=1}^{S}\left[\derivative{N_{i\setminus
0}}{\varepsilon}\right]^{2}$ and
$\expectationvalue{\left[\derivative{R_{0}}{\varepsilon}\right]^{2}}=M^{-1}\sum_{\alpha=1}^{M}\left[\derivative{R_{\alpha\setminus
0}}{\varepsilon}\right]^{2}$, which follows from the assumption that the
system self averages, we obtain:

$\displaystyle\expectationvalue{\left[\derivative{\overline{N}_{0}^{+}}{\varepsilon}\right]^{2}}$
$\displaystyle=\frac{1/M}{(\sigma_{e}\rho\chi)^{2}}\sum_{\alpha,\beta:\overline{R}_{\alpha}>0,\overline{R}_{\beta}>0}\expectationvalue{d_{0\alpha}d_{0\beta}}\left(\derivative{\overline{R}_{\alpha\setminus
0}^{+}}{\varepsilon}\derivative{\overline{R}_{\beta\setminus
0}^{+}}{\varepsilon}+\derivative{\overline{R}_{\alpha\setminus
0}^{+}}{\varepsilon}\expectationvalue{\eta_{\beta}^{(R)}}\right.$ (64)
$\displaystyle\hskip
170.71652pt\left.+\expectationvalue{\eta_{\alpha}^{(R)}}\derivative{\overline{R}_{\beta\setminus
0}^{+}}{\varepsilon}+\expectationvalue{\eta_{\alpha}^{(R)}\eta_{\beta}^{(R)}}\right)$
$\displaystyle=\frac{1/M}{(\sigma_{e}\rho\chi)^{2}}\sum_{\alpha,\beta:\overline{R}_{\alpha},\overline{R}_{\beta}>0}\delta_{\alpha\beta}\left(\derivative{\overline{R}_{\alpha}\setminus}{\varepsilon}\derivative{\overline{R}_{\beta\setminus
0}}{\varepsilon}+\delta_{\alpha\beta}\right)$ (65)
$\displaystyle=\frac{\phi_{R}}{(\sigma_{e}\rho\chi)^{2}}\left(\expectationvalue{\left[\derivative{\overline{R}_{0}^{+}}{\varepsilon}\right]^{2}}+1\right),$
(66)
$\displaystyle\expectationvalue{\left[\derivative{\overline{R}_{0}^{+}}{\varepsilon}\right]^{2}}$
$\displaystyle=\frac{\sigma_{e}^{2}/M}{(1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu)^{2}}\sum_{i,j:\overline{N}_{i}>0,\overline{N}_{j}>0}\left(\derivative{\overline{N}_{i\setminus
0}^{+}}{\varepsilon}\derivative{\overline{N}_{j\setminus
0}^{+}}{\varepsilon}+\derivative{\overline{N}_{i\setminus
0}^{+}}{\varepsilon}\expectationvalue{\eta_{j}^{(N)}}\right.$ (67)
$\displaystyle\hskip
11.38092pt\left.+\expectationvalue{\eta_{i}^{(N)}}\derivative{\overline{N}_{j\setminus
0}^{+}}{\varepsilon}+\expectationvalue{\eta_{i}^{(N)}\eta_{j}^{(N)}}\right)$
$\displaystyle\qquad\times\left(\rho^{2}\expectationvalue{d_{i0}d_{j0}}+\rho\sqrt{1-\rho^{2}}\left(\expectationvalue{d_{i0}}\expectationvalue{x_{j0}}+\expectationvalue{x_{i0}}\expectationvalue{d_{j0}}\right)+(1-\rho^{2})x_{i0}x_{j0}\right)$
$\displaystyle=\frac{\sigma_{e}^{2}(S/M)/S}{(1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu)^{2}}\sum_{i,j:\overline{N}_{i}>0,\overline{N}_{j}>0}\left(\derivative{\overline{N}_{i\setminus
0}^{+}}{\varepsilon}\derivative{\overline{N}_{j\setminus
0}^{+}}{\varepsilon}+\delta_{ij}\right)\left(\rho^{2}\delta_{ij}+(1-\rho^{2})\delta_{ij}\right)$
$\displaystyle=\frac{\sigma_{e}^{2}\gamma^{-1}\phi_{N}}{(1-\rho\sigma_{e}\sigma_{c}\gamma^{-1}\nu)^{2}}\left(\expectationvalue{\left[\derivative{\overline{N}_{0}^{+}}{\varepsilon}\right]^{2}}+1\right).$
(68)

Solving this system of equations yields:

$\displaystyle\boxed{\expectationvalue{\left[\derivative{\overline{N}_{0}^{+}}{\varepsilon}\right]^{2}}=\frac{\phi_{R}\left((1-\nu\rho\sigma_{c}\sigma_{e}\gamma^{-1})^{2}\sigma_{e}^{-2}+\gamma^{-1}\phi_{N}\right)}{\left[\rho\chi\left(1-\nu\rho\sigma_{c}\sigma_{e}\gamma^{-1}\right)\right]^{2}-\gamma^{-1}\phi_{N}\phi_{R}}},$
(69)
$\displaystyle\boxed{\expectationvalue{\left[\derivative{\overline{R}_{0}^{+}}{\varepsilon}\right]^{2}}=\frac{\gamma^{-1}\phi_{N}\left(\phi_{R}+(\rho\sigma_{e}\chi)^{2}\right)}{\left[\rho\chi\left(1-\nu\rho\sigma_{c}\sigma_{e}\gamma^{-1}\right)\right]^{2}-\gamma^{-1}\phi_{N}\phi_{R}}}.$
(70)

Figure B11:  Heatmaps of the variances of the response of a surviving species
and resources to a random perturbation (Eqs. 69, 70) as a function of the
parameters $\sigma$ and $\rho$ with other parameters matching those in Fig.
3(a). The instability boundary is shown as a gray solid line, and the
infeasibility boundary is shown as a gray dashed line.  Figure B12:  Log-log
plot of the absolute value of the variances of the susceptibilities
$\langle\left(dN_{0}^{+}/\mathrm{d}\varepsilon\right)^{2}\rangle$,
$\langle\left(dR_{0}^{+}/\mathrm{d}\varepsilon\right)^{2}\rangle$ plotted as a
function of the distance from the instability boundary, $|\rho-\rho^{\star}|$.
The data is exactly that which is used in Fig. 3(b). A fit $\propto
1/|\rho-\rho^{\star}|$ is shown for each susceptibility with dashed lines.

These susceptibilities are the variance of the response of a surviving species
or resource when the system is subject to a random perturbation. The
divergence of these variances represent the breakdown of the mean-field
approximation. Further, when these susceptibilities diverge, the uninvadable
fixed point predicted by the mean-field approximation becomes dynamically
unstable. Therefore, when these susceptibilities are divergent, whenever
invasion is attempted, the system will exhibit dynamical fluctuations. We call
the boundary in parameter space at which these susceptibilities diverge the
instability boundary. These susceptibilities diverge when
$0=\left[\rho\chi\left(1-\nu\rho\sigma_{c}\sigma_{e}\gamma^{-1}\right)\right]^{2}-\gamma^{-1}\phi_{N}\phi_{R}.$
This condition along with the cavity self-consistency equations (Eqs. 45, 46,
47, 48, 50, 51, 52, 53) determines a co-dimension-one boundary in the space of
model parameters. By using the relations in line 49, we obtain the following
relation for the boundary at which the system becomes unstable:

$\displaystyle\phi_{R}^{\star}\left((\rho^{\star})^{2}\phi_{R}^{\star}-\gamma^{-1}\phi_{N}^{\star}\right)=0\implies\left(\rho^{\star}\right)^{2}=\gamma^{-1}\frac{\phi_{N}^{\star}}{\phi_{R}^{\star}}$
(71) $\displaystyle\Downarrow$ $\displaystyle\boxed{\text{instability boundary
occurs when:}\quad\left(\rho^{\star}\right)^{2}=\frac{(\text{\\# of surviving
species})^{\star}}{(\text{\\# of non-depleted resources})^{\star}},}$ (72)

where $X^{\star}$ denotes the value of a quantity $X$ at the instability
boundary. We are able to determine the critical level of asymmetry,
$\rho^{\star}$, at which the ecosystem becomes unstable by solving a nonlinear
least squares problem (see Appendix D.3). The variances of the
susceptibilities
$\langle\left(\mathrm{d}N_{0}^{+}/\mathrm{d}\varepsilon\right)^{2}\rangle$,
$\langle\left(\mathrm{d}R_{0}^{+}/\mathrm{d}\varepsilon\right)^{2}\rangle$ are
shown for various values of $\sigma=\sigma_{c}=\sigma_{e}$ and $\rho$ in Fig.
B11 with the instability and infeasibility boundaries overlaid. These
variances of susceptibilities additionally have an asymptotic power law
dependence on the distance from the instability boundary,
$|\rho-\rho^{\star}|$, as shown in Fig. B12.

##  Appendix C Chaotic nature of dynamics

In the dynamic phase of the aMCRM, when an ecosystem is sufficiently large,
the dynamics are chaotic. There are classical definitions, such as the
Lyapunov exponent, which can be used to determine whether a system is chaotic.
However, chaos is a subtle concept, and the Lyapunov exponent and dependence
on initial conditions does not always give the whole picture. For this reason,
there are other tests for chaos and predictability of a dynamical system
which, in conjunction with more traditional tests, can be used to determine
the nature of the chaos. As shown in Fig. 4(b) in the main text, the
trajectories of the abundances of the surviving species and resources diverge
from one another, indicating that the dynamics are chaotic in the classical
sense. The aMCRM in the dynamic phase exhibits other characteristics of chaos,
providing further evidence that the dynamics are chaotic and unpredictable.

For example, in Fig. C13, we show that the dynamics of the aMCRM in the
dynamic phase are unpredictable in the sense that a trajectory does not
display any clear patterns in the projection onto the first three principal
components of the correlation matrix of the time series of the abundances of
the surviving species and resources. Such a projection is necessary to
visualize the dynamics of the system, as they are very high-dimensional.

Figure C13:  Projection of $M+S=256+256=512$-dimensional dynamics onto the 3
highest-ranked eigenvectors of the correlation matrix (first 3 principal
components) of the time series of the abundances of the surviving species and
resources. The lack of clustering of the trajectories in the projection
indicates that the dynamics are chaotic and unpredictable. The parameters
correspond to the gray star in Fig. 3(a). Before plotting, the simulation was
run for $10^{4}$ time units to eliminate any potential transients.

###  C.1 Analysis of the Lyapunov exponents

In Fig. 4(a), we show the maximal Lyapunov exponent for simulations classified
by whether they reach a steady state for various values of $\rho$. To
determine the Lyapunov spectrum, we use the lyapunovspectrum function in the
ChaosTools.jl Julia package which employs the ‘H2’ method of Geist, originally
stated in Benettin et al. [37, 36, 34, 35]. This algorithm is described in
detail in Appendix A of [36], and its applications are described in Chapter
3.2. Conceptually, in order to compute the $k$th largest Lyapunov exponent,
the algorithm evolves $n\geq k$ deviation vectors in the tangent space of the
system and computes how the shape of the $k$-dimensional parallelepiped
spanned by the deviation vectors evolve; the eigenvalues of the matrix
describing the evolution of the parallelepiped are asymptotically related to
the Lyapunov exponents. The Jacobian of the system when the dynamical variable
is considered as $(N_{1},\dots,N_{S},R_{1},\dots,R_{M})$ is:

$\displaystyle\mathbf{J}(N_{1},\dots,N_{S},R_{1},\dots,R_{M})=\left[\begin{array}[]{c|c}\frac{\partial\dot{N}_{i}}{\partial
N_{j}}&\frac{\partial\dot{N}_{i}}{\partial R_{\beta}}\\\
\hline\cr\frac{\partial\dot{R}_{\alpha}}{\partial
N_{j}}&\frac{\partial\dot{R}_{\alpha}}{\partial
R_{\beta}}\end{array}\right]=\left[\begin{array}[]{c|c}\delta_{ij}(\sum_{i=1}^{S}c_{i\alpha}R_{\alpha}-m_{i})&N_{i}c_{i\beta}\\\
\hline\cr-
e_{j\alpha}R_{\alpha}&\delta_{\alpha\beta}(K_{\alpha}-2R_{\alpha}-\sum_{i=1}^{S}c_{i\alpha}N_{i})\end{array}\right].$
(77)

The Lyapunov spectrum is computed with $M=S=512$ species and resources, and
the initial deviation vectors are chosen to be $8$ unit vectors in the
direction of randomly chosen species and resources. The parameters passed to
the lyapunovspectrum function are the number of steps, N = 1000, the time
step, $\Delta$t = 5, and the time to wait before starting to record the
Lyapunov spectrum, Ttr = 1000. The initial conditions are $N_{i}(0)=1/S$ and
$R_{\alpha}(0)=1/M$.

As seen in Fig. 4(a), the maximal Lyapunov exponent for a simulation that does
not reach a steady state is positive while the maximal Lyapunov exponent for a
simulation that does reach a steady state is negative. A simulation is
categorized as reaching a steady state if the mean absolute value of all
derivatives of the abundances at the end of the simulation is less than
$10^{-6}$. These results indicates that the dynamics are chaotic in the
classical sense for the simulations which have persistent fluctuations. The
impact of system size on the probability of reaching a steady state is
discussed in Appendix E.

###  C.2 Analysis of the generalized alignment index (GALI)

Figure C14:  Generalized alignment indices (GALI) for simulations in the
dynamic and stable phases with $M=S=512$ species and resources. The
simulations all have the same sampled parameters; only $\rho$ and order of the
GALI are varied. For the simulations in the dynamic phase (solid lines), the
GALI asymptotically decays exponentially with time, indicating chaos, while
for the simulations in the stable phase (dashed lines), the GALI
asymptotically approaches a non-zero value, indicating the dynamics achieve a
steady state. Simulation parameters for the GALI in the stable and dynamic
phases correspond to the gray dot and star in Fig. 3(a), respectively.

An alternative method to determine whether a system is chaotic is to use the
generalized alignment index (GALI) [38, 37, 36]. The GALI is a measure of how
vectors in the tangent space of a trajectory align with each other. Let
$\hat{w}_{1}(0),\dots,\hat{w}_{k}(0)\in\mathbb{R}^{M+S}$ be linearly-
independent unit deviation vectors that evolve according to,

$\displaystyle\derivative{\mathrm{d}t}\hat{w}_{i}(t)=\mathbf{J}(x(t))\hat{w}_{i}(t),$
(78)

where $x(t)=(N_{1}(t),\dots,N_{S}(t),R_{1}(t),\dots,R_{M}(t))$ and
$\mathbf{J}(x(t))$ is the Jacobian matrix of the system in state $x(t)$. These
deviation vectors are normalized to have unit length at regular (small) time
intervals. The order-$k$ GALI is defined as,

$\displaystyle\mathrm{GALI}_{k}(t)=\left\|\hat{w}_{1}(t)\wedge\hat{w}_{2}(t)\wedge\cdots\wedge\hat{w}_{k}(t)\right\|=|\det(\hat{w}_{1}(t),\hat{w}_{2}(t),\dots,\hat{w}_{k}(t))|,$
(79)

which is the volume of the $k$-dimensional parallelepiped spanned by the time-
evolved deviation vectors. If the system is chaotic, the GALI will decay
exponentially with time, indicating that the deviation vectors are becoming
more aligned with each other due to the exponential divergence of nearby
trajectories. If the system is not chaotic and reaches a steady state, the
GALI will asymptotically approach a non-zero value, indicating that the
deviation vectors form a non-zero volume in the tangent space because no
single direction dominates the dynamics. Alternatively, if the system is
asymptotically periodic or quasi-periodic and does not reach a steady state,
the GALI will decay as a power law with time.

In Fig. C14, we show the GALI for simulations in the dynamic and stable phases
for $k=2,\dots,7$ with $M=S=512$ species and resources. We use the gali
function from the ChaosTools.jl Julia package [37]. The simulations all have
the same sampled parameters; only $\rho$ and the order of the GALI are varied.
For the simulations in the dynamic phase (solid lines), the GALI
asymptotically decays exponentially with time, indicating chaos, while for the
simulations in the stable phase (dashed lines), the GALI asymptotically
approaches a non-zero value, indicating the dynamics achieve a steady state.
The initial deviation vectors are chosen to be $k$ unit vectors in the
direction of randomly chosen species and resources. The arguments passed to
the gali function are the run time, T = 2e4, the time step, $\Delta$t = 5.,
the threshold at which to stop the simulation, threshold = 1e-22, and the
order of the GALI, k.

##  Appendix D Simulation and numerical methods

###  D.1 Numerical integration

Numerical integration of the differential equations is performed using the
Tsit5 solver from the DifferentialEquations.jl package [28] in the Julia
programming language [29]. Tsit5 is the Tsitouras 5/4 Runge-Kutta method, a
5th order Runge-Kutta method with an embedded 4th order method for error
estimation and step size control [30]. The solver is configured to use a
relative tolerance of $10^{-14}$ and an absolute tolerance of $10^{-14}$. The
aMCRM differential equations are also effectively integrated using the VCABM
solver which an adaptive order adaptive time Adam-Moulton method [31, 28].
When using VCABM, the solver is configured to use a relative tolerance of
$10^{-11}$ and an absolute tolerance of $10^{-11}$. The VCABM is often faster
than Tsit5 for the aMCRM differential equations, but is less stable for some
parameter values, particularly in the dynamic phase with low $\rho$ and high
$\sigma$. Therefore, in nearly all cases, especially those where the parameter
space is explored, Tsit5 is used.

###  D.2 Surviving/depletion cutoffs, numerical instability, and small
immigration

Simulating ecological dynamics often involves including very small immigration
rates, $\lambda_{i}$, $\lambda_{\alpha}$ for species and resources,
respectively, to numerically regularize dynamics and ensure an uninvadable
steady state is approached in simulations. Furthermore, the inclusion of small
immigration ensures that chaotic fluctuations in the dynamics do not drive
species or resources to extinction, killing chaos due to finite size effects.
The equations used for numerical simulations are,

$\displaystyle\derivative{N_{i}}{t}$
$\displaystyle=N_{i}\left(\sum_{\alpha=1}^{M}c_{i\alpha}R_{\alpha}-m_{i}\right)+\lambda_{i},$
(80) $\displaystyle\derivative{R_{\alpha}}{\mathrm{d}t}$
$\displaystyle=R_{\alpha}\left(K_{\alpha}-R_{\alpha}\right)-\sum_{i=1}^{S}N_{i}e_{i\alpha}R_{\alpha}+\lambda_{\alpha}.$
(81)

In all our simulations, we take $\lambda_{i},\lambda_{\alpha}=10^{-10}$.
Including this small level of immigration regularizes the dynamics in that the
simulated steady-state abundances of species and resources now have a clear
gap between the abundances of those that are surviving and those that are
extinct, as shown in Fig. D15. This gap is used to define a cutoff to identify
fractions of surviving species and resources. Including this small level of
immigration additionally stabilizes the numerical integrator, leading to fewer
simulations that report errors due to numerical instability. The effects of
immigration on dynamics in the random Generalized Lotka–Volterra model is
discussed in detail in [20]; we expect that many of the intuitions presented
there apply to the aMCRM as well.

Figure D15:  Log histogram of the species and resource abundances at steady
state for a system with $M,S=256$ and $\rho=0.9$. There is a clear gap in
abundance between species/resources that are surviving and those that are
extinct. The immigration rate used is $\lambda=10^{-10}$.

###  D.3 Numerically solving the cavity self-consistency equations and
calculating phase boundaries

Figure D16:  Least-squares objectives for solving the cavity self-consistency
equations. The solid curve is the instability boundary, and the dashed curve
is the infeasibility boundary. The region in which the objective function is
very small is the region in which the cavity self-consistency equations are
solved; the region in which the objective function is large represents invalid
solutions which are not included in plots. The parameters used are those used
in Fig. 3(a).

In order to solve the self-consistency equations for the cavity method, we the
use nonlinear least squares method. We define the objective function to be the
sum of the squares of the differences between the left and right-hand sides of
the self-consistency equations (45, 46, 50, 51, 52, 53, 48, 47):

$\displaystyle L(\phi_{R},\phi_{N},\langle N\rangle,\langle
R\rangle,q_{N},q_{R},\chi,\nu)=(\text{LHS $-$ RHS of
Eq.~{}\eqref{eq:phiNselfconsist}})^{2}+\cdots+(\text{LHS $-$ RHS of
Eq.~{}\eqref{eq:nuselfconsist}})^{2}.$ (82)

When the self-consistency equations are solved, the objective function is
zero. The objective is minimized using adaptive differential evolution solver,
implemented as
BBO_adaptive_de_rand_1_bin_radiuslimited in the Julia package BlackBoxOptim.jl
[32].

There are choices of parameters
$\mu_{c},\mu_{e},\sigma_{c},\sigma_{e},K,\sigma_{K},m,\sigma_{m},\rho$ for
which the objective function does not reach zero. In Fig. D16, we plot the
objective function values for various choices of $\rho$ and
$\sigma=\sigma_{e}=\sigma_{c}$; the other parameters for this figure are those
used in Fig. 3. We can see that the objective function values are very small
for $\rho$ near one and suddenly jumps to a much larger value for $\rho$ near
zero. The cause of this increase in the objective function is the term (LHS
$-$ RHS of Eq. (47)) in the objective function. Observe that as $\chi\to 0$,
$-\nu\to\infty$, which leads to the objective function value blowing up. In
order to find the boundary at which the objective function blows up, we simply
add a term to the objective function that is zero when $\chi=0$ and allow one
model parameter (here, $\rho$) to vary:

$\displaystyle L^{\text{CSI}}(\phi_{R},\phi_{N},\langle N\rangle,\langle
R\rangle,q_{N},q_{R},\chi,\nu,\rho^{\text{CSI}})=$ $\displaystyle
L(\phi_{R},\phi_{N},\langle N\rangle,\langle
R\rangle,q_{N},q_{R},\chi,\nu,\rho^{\text{CSI}})$ (83)
$\displaystyle\qquad\qquad\qquad+(\chi-0)^{2}.$

Similarly, to find the instability boundary, we add a term to the objective
function corresponding to the instability condition (Eq. (72)):

$\displaystyle L^{\star}(\phi_{R},\phi_{N},\langle N\rangle,\langle
R\rangle,q_{N},q_{R},\chi,\nu,\rho^{\star})=$ $\displaystyle
L(\phi_{R},\phi_{N},\langle N\rangle,\langle
R\rangle,q_{N},q_{R},\chi,\nu,\rho^{\star})$ (84)
$\displaystyle\qquad\quad+(\text{LHS $-$ \text{RHS} of
Eq.~{}\eqref{eq:instability-boundary}})^{2}.$

Minimizing these objective functions while allowing one additional model
parameter to vary gives the infeasibility and instability boundaries shown in
all figures in this paper.

##  Appendix E Finite size scaling analysis

In order to understand whether the aMCRM achieves a steady state in the
thermodynamic limit, we must understand how finite system size and a finite
simulation time impact simulation results. To assess the existence and
dynamical nature of a steady state in a simulation, we numerically calculate a
time to steady state (TtSS). After running a simulation, we iterate through
the time series and find the time at which the mean absolute value of the
derivatives last dips below a threshold, $10^{-10}$. The tolerance of the
solver (Tsit5 [30, 28]) is $10^{-14}$, so this threshold is well above the
numerical error of the simulation. Because we simulate to a finite time
$T_{\text{max}}=2^{16}+10^{4}=75536$, it is not possible to distinguish
between a system that will reach a steady state at some time beyond
$T_{\text{max}}$ and a system that will never reach a steady state.
Furthermore, when performing these simulations, especially with large system
sizes, errors occasionally occur in the solver, and the simulation is aborted;
while we do not present finite-size scaling results in regions where the
simulations are frequently aborted, a robust analysis should account for this.
Here, we develop a custom maximum likelihood estimation (MLE) technique to fit
the distributions TtSS’s at different system sizes, taking into account these
two situations. We then examine the scaling of the resulting fit parameters to
determine their behavior in the thermodynamic limit.

###  E.1 Modeling simulation outcomes

To construct a MLE model, we first must model the statistical process which
determines each possible outcome of a simulation. We observe that simulation
outcomes fall into three categories:

  1. (i)

The simulation reaches steady state within the maximum simulation time
$T_{\max}$.

  2. (ii)

The simulation does not reach steady state within $T_{\max}$.

  3. (iii)

The solver encounters an error and the simulation aborts.

First, we represent the probability of a simulation encountering an error as
$p_{\text{err}}$. Next, we consider whether a simulation could possibly reach
a steady state if given enough simulation time. We represent the probability
of a simulation reach steady state in the infinite time limit as
$\phi_{\text{SS}}$. Finally, we define the probability density function
$p_{\text{SS}}(t)$ of TtSS’s for those simulations that can reach steady state
when given enough time.

Using these definitions, we can now construct the probabilities for each of
the cases defined above. Let $T$ be a random variable representing the
observed outcome of a simulation with possible values
${T\in[0,T_{\max}]\cup\\{\texttt{noSS},\texttt{error}\\}}$, where $T$ takes a
finite value in the interval $[0,T_{\max}]$ for case (i), or the special
values noSS or error in cases (ii) and (iii), respectively.

To model case (i), we note that it results in an outcome represented by a
continuous numerical value, while cases (ii) and (iii) represent categorical
outcomes. For ease of derivation, we also convert $T$ to a categorical outcome
in the prior case by artificially breaking our simulation interval into small
discrete time steps of size $\Delta t$. Later, we will take the limit $\Delta
t\rightarrow 0$, removing the dependency of our model on this step size. The
probability of a simulation completing successfully and reaching a steady
state in the time interval $t<T<t+\Delta t$ with $T<T_{\max}$ is approximately

${\mathrm{P}}(t<T\leq t+\Delta t\mid
T<T_{\max})\approx(1-p_{\text{err}})\phi_{\text{SS}}\Delta
t\,p_{\text{SS}}(t).$ (85)

For case (ii), simulations do not encounter an error, but do not reach steady
state within $T_{\max}$. We must consider two separate possibilities: a
simulation may never reach steady state, even if $T_{\max}\rightarrow\infty$,
or would reach steady-state after further simulation, $T>T_{\max}$. In the
first case, the probability is simply
$(1-p_{\text{err}})(1-\phi_{\text{SS}})$. If a simulation would eventually
have reached steady-state if we had continued simulating, we must use our
probability density function $p_{\text{SS}}(t)$. However, we do not know
exactly when the simulation would have finished, so we must sum across all
discrete time intervals greater than $T_{\max}$. The total probability for
case (ii) is then

$\displaystyle{\mathrm{P}}(\texttt{noSS})$
$\displaystyle\approx(1-p_{\text{err}})(1-\phi_{\text{SS}})+(1-p_{\text{err}})\phi_{\text{SS}}\sum_{n=0}^{\infty}\Delta
t\,p_{\text{SS}}(T_{\max}+n\Delta t)$ (86)
$\displaystyle\approx(1-p_{\text{err}})\quantity(1-\phi_{\text{SS}}F_{\text{SS}}(T_{\max})),$

where we have taken the limit $\Delta t\rightarrow 0$ and simplified, defining
$F_{\text{SS}}(T)=\int_{0}^{T}dtp_{\text{SS}}(t)$ as the cumulative
distribution function of $p_{\text{SS}}(t)$.

Finally, for case (iii), the probability is simply

${\mathrm{P}}(\texttt{error})=p_{\text{err}}.$ (87)

###  E.2 Likelihood function

We now use the outcome probabilities in the previous section to construct a
likelihood function to fit the data. We are given $N$ (independent)
simulations resulting in observed TtSSs $\\{T_{1},\cdots,T_{N}\\}$ which may
take on values ${T_{i}\in[0,T_{\max}]\cup\\{\texttt{noSS},\texttt{error}\\}}$.
Our fit parameters $\theta$ include $p_{\text{err}}$, $\phi_{\text{SS}}$,
along with the parameters that define the probability density function
$p_{\text{SS}}(t)$. Defining the likelihood of the data given the parameters
$\theta$ as ${\mathrm{P}}(\\{T_{i}\\}_{i=1}^{N}|\theta)$, the MLE estimate of
the parameters is found by minimizing the negative log-likelihood,

$\theta^{\star}=\operatorname*{argmax}\limits_{\theta}{\mathrm{P}}(\\{T_{i}\\}_{i=1}^{N}|\theta)=-\operatorname*{argmin}\limits_{\theta}\frac{1}{N}\log{\mathrm{P}}(\\{T_{i}\\}_{i=1}^{N}|\theta).$
(88)

Using our outcome probabilities, the negative log-likelihood is then

$\displaystyle\mathcal{L}(\\{T_{i}\\}_{i=1}^{N};\theta)=$
$\displaystyle-\frac{1}{N}\log{\mathrm{P}}(\\{T_{i}\\}_{i=1}^{N}|\theta)=-\frac{1}{N}\log\prod_{i=1}^{N}{\mathrm{P}}(T_{i}|\theta)=-\frac{1}{N}\sum_{i=1}^{N}\log{\mathrm{P}}(T_{i}|\theta)$
(89) $\displaystyle=$
$\displaystyle-\frac{1}{N}\sum_{i=1:T_{i}<T_{\max}}^{N}\log\quantity[(1-p_{\text{err}})\phi_{\text{SS}}\Delta
t\,p_{\text{SS}}(T_{i})]$
$\displaystyle-\frac{1}{N}\sum_{i=1:T_{i}=\texttt{noSS}}^{N}\log\quantity[(1-p_{\text{err}})\quantity(1-\phi_{\text{SS}}F_{\text{SS}}(T_{\max}))]-\frac{1}{N}\sum_{i=1:T_{i}=\texttt{noSS}}^{N}\log
p_{\text{err}}.$

Next, we define $f_{\text{err}}$ as the fraction of simulations that encounter
errors, and the $f_{\text{SS}}$ as the fraction of simulations that did not
encounter errors and reached a steady state within time $T_{\max}$. Using
these definitions, we get

$\displaystyle\mathcal{L}(\\{T_{i}\\}_{i=1}^{N};\theta)=$
$\displaystyle-\frac{1}{N}\sum_{i=1:T_{i}<T_{\max}}^{N}\log
p_{\text{SS}}(T_{i})-f_{\text{SS}}\log\quantity[(1-p_{\text{err}})\phi_{\text{SS}}\Delta
t]$ (90)
$\displaystyle-(1-f_{\text{err}}-f_{\text{SS}})\log\quantity[(1-p_{\text{err}})\quantity(1-\phi_{\text{SS}}F_{\text{SS}}(T_{\max}))]-f_{\text{err}}\log
p_{\text{err}}$ $\displaystyle=$
$\displaystyle-\frac{1}{N}\sum_{i=1:T_{i}<T_{\max}}^{N}\log
p_{\text{SS}}(T_{i})-(1-f_{\text{err}}-f_{\text{SS}})\log\quantity(1-\phi_{\text{SS}}F_{\text{SS}}(T_{\max}))$
$\displaystyle-
f_{\text{SS}}\log\phi_{\text{SS}}-(1-f_{\text{err}})\log(1-p_{\text{err}})-f_{\text{err}}\log
p_{\text{err}}-f_{\text{SS}}\log\Delta t.$

We note that the last term is a constant, so we may choose to ignore it and
take the limit $\Delta t\rightarrow 0$ as mentioned previously, giving us the
final form for the likelihood,

$\displaystyle\mathcal{L}(\\{T_{i}\\}_{i=1}^{N};\theta)=$
$\displaystyle-\frac{1}{N}\sum_{i=1:T_{i}<T_{\max}}^{N}\log
p_{\text{SS}}(T_{i})-(1-f_{\text{err}}-f_{\text{SS}})\log\quantity(1-\phi_{\text{SS}}F_{\text{SS}}(T_{\max}))$
(91) $\displaystyle-
f_{\text{SS}}\log\phi_{\text{SS}}-(1-f_{\text{err}})\log(1-p_{\text{err}})-f_{\text{err}}\log
p_{\text{err}}$

###  E.3 Fitting the model

To find the maximum likelihood estimates of the parameters, we minimize Eq.
(91) using the BFGS algorithm implemented in the Julia package Optimization.jl
[33]. We found empirically that the Fréchet distribution is a good choice for
$p_{\text{SS}}(t)$. The Fréchet distribution has the following forms for the
cumulative distribution function and probability density function:

$\displaystyle F_{\text{SS}}(t)=\begin{cases}e^{-(t/\tau)^{-\alpha}},&t>0,\\\
0,&t\leq 0,\end{cases}\qquad
p_{\text{SS}}(t)=\begin{cases}\frac{\alpha}{\tau}\left(\frac{t}{\tau}\right)^{-\alpha-1}e^{-(t/\tau)^{-\alpha}},&t>0,\\\
0,&t\leq 0,\end{cases}$ (92)

where $\tau$ is the timescale parameter and $\alpha$ is the shape parameter.
The errors in the fitted parameters
$\phi_{\text{SS}},p_{\text{err}},\tau,\alpha$ are computed by bootstrapping
the data. That is, we randomly sample $N$ data points from the simulation data
set with replacement and fit the model to the sampled data repeatedly,
providing us with an empirical distribution for each fit parameter.

###  E.4 Scaling of parameters with system size and discussion

The simulated data for different system sizes in the stable and dynamic phases
can be visualized by analyzing the empirical cumulative distribution functions
(CDFs) of the observed times to steady state called $T_{\text{SS}}$ above. In
Fig. E17, we show the CDFs for the stable and dynamic phases for $4\leq
M,S\leq 1024$ along with the CDFs for the fitted distributions. A clear
difference is apparent between the stable and dynamic phases, and the fitted
distributions match the simulated data well. A scaling collapse of the CDF
based on the fitted distributions is shown in Fig. E18.

The fit parameters for the stable and dynamic phases are shown in Fig. E19. In
the stable phase, the estimated probability of reaching steady state,
$\phi_{\text{SS}}$ remains approximately equal to $1$ for all system sizes,
and the estimated timescale $\tau$ asymptotically increases with system size
as a power law. In the dynamic phase, $\phi_{\text{SS}}$ decreases with system
size towards $0$, and $\tau$ asymptotically increases with system size also as
a power law, but with a larger exponent. In both phases, $\alpha$ approaches
$1$ asymptotically from above.

We provide fits of $\phi_{\text{SS}}$ as a function of system size $M$ using
curves of the form
$\phi_{\text{SS}}(M)=\phi_{\text{max}}-\Delta(1-\exp(-(M/\xi)^{\kappa}))$. For
the dynamic phase $\phi_{\text{max}}=1.0$, $\Delta=1.0$, $\xi=278.5$, and
$\kappa=1.5$; for the stable phase, $\phi_{\text{max}}=1.0,\Delta=0.0$,
$\kappa$ remains approximately equal to the value at which it is initialized
in the optimization algorithm, and $\xi$ approaches the maximum value set in
the optimization algorithm.

We also provide fits of the timescale $\tau$ as a function of system size $M$
using power laws of the form $\tau(M)=bM^{a}$. Curves are fit for data
$M,S\geq 2^{9.5}$. In the stable phase, $a=0.4$ and $b=580$; in the dynamic
phase, $a=0.9$ and $b=70$.

This finite size scaling is run in the range $512\leq M,S\leq 1024$ but is not
shown in Fig. E19 because solver errors are occasionally present. More
significantly, as very few simulations reach steady state in the dynamic phase
when the system size is sufficiently large, the MLE procedure cannot clearly
distinguish between cases where $\tau$ is very large or $\phi_{\text{SS}}$ is
very small; this ambiguity can be seen more clearly in Fig. E17. With the
prior knowledge of the finite size scaling as shown in Fig. E19 where the MLE
procedure succeeds and there are no solver errors, we conclude that at these
larger system sizes, it is both the case that $\tau$ grows large (in
polynomial order of system size) and $\phi_{\text{SS}}$ approaches either zero
or one. That fact both the timescale $\tau$ diverges and $\phi_{\text{SS}}$
approaches either zero at large system size leads us to conclude that in the
thermodynamic limit, all systems in the stable phase eventually reach steady
state for long enough times, while in the dynamic phase, no systems reach
steady state.

Figure E17:  Empirical cumulative distribution functions of $T$, the time to
steady state, for simulations of different system sizes in the stable phase
(left) and the dynamic phase (right). The maximum likelihood estimate
distributions as described in Appendix E for the time to steady state for the
various system sizes are shown with dashed lines. The maximum value shown on
the horizontal axis is the maximum time simulated,
$T_{\text{max}}=2^{16}+10^{4}=75536$.  Figure E18:  Scaling collapse of TtSS
CDFs. CDFs of observed TtSSs are scaled by the MLE-fitted value of
$\phi_{\text{SS}}$ where the horizontal axis is scaled by the MLE-fitted
values of $\tau$, $\alpha$. Note that the curves collapse onto each other in
both phases for system sizes that are moderately large $M,S\gtrapprox
2^{4.5}$. This indicates that we have successfully identified the asymptotic
scaling of the form of the TtSS distribution with system size.  Figure E19:
Impact of finite size effects on the dynamics of the aMCRM. Rows correspond to
different parameters describing dynamics: the estimated probability a
simulation of a given system size will reach steady state (top) and the
timescale of the time to steady state of simulations of a given system size
that do reach steady state (bottom). Columns correspond to simulation data
from different phases: the stable phase (left) and the dynamic phase (right);
the parameters corresponding these simulations are marked in the phase diagram
in Fig. 3 with a circle and star, respectively. The probability of reaching
steady state in the stable phase remains at $1$ for all system sizes, while
the probability of reaching steady state in the unstable phase decreases to
zero with increasing system size, indicating that stable dynamics are not
observed within the dynamic phase in the thermodynamic limit. The timescale of
the time to steady state increases with system size in both the stable and
dynamic phases.

##  Appendix F Parameters in figures

#### Figure 1 (schematic) details

Species drawings are traces of images generated using Adobe Firefly Generative
AI. Prompts: “simple 2D flat unshaded vector graphics of microbes,” “simple 2D
flat digital art sketches of bacteria and microbes in icon-style with no
shading.”

#### Figure 2 (example dynamics) simulation parameters

Both simulations have random variables
$d_{i\alpha},x_{i\alpha},K_{\alpha},m_{i}$ drawn from standard normal
distributions. The parameters are,

$\displaystyle K=1,\quad\sigma_{K}=0.1,\quad
m=1,\quad\sigma_{m}=0.1,\quad\mu_{c}=200,\quad\sigma_{c}=3.5,\quad\mu_{e}=200,\quad\sigma_{e}=3.5,\quad\lambda_{i}=\lambda_{\alpha}=10^{-10},$
(93)

with $M=S=256$ resources and species. For dynamics in the stable phase,
$\rho=0.9$, and for the dynamics in the unstable phase, $\rho=0.5$. Between
the two simulations, the model parameters ($\delta K_{\alpha},\delta
m_{i},d_{i\alpha},x_{i\alpha}$) are sampled once and only the level of
nonreciprocity ($\rho$) is changed. For these choices of parameters,
$\rho^{\star}=0.72$. The initial conditions are $N_{i}(0)=1/S$ and
$R_{\alpha}(0)=1/M$. Numerical integration is performed using a ‘fourth-order,
five-stage explicit Runge-Kutta method with embedded error estimator of
Tsitouras’ implemented as Tsit5 in the Julia DifferentialEquations.jl package
[28, 30]

#### Figure 3 (phase diagram) simulation parameters

Each point in the heatmap represents an average over 32 simulations where the
parameters were independently sampled from normal distributions:
$d_{i\alpha},x_{i\alpha},K_{\alpha},m_{i}\sim\mathcal{N}(0,1)$. The
distribution parameters are,

$\displaystyle K=1,\quad\sigma_{K}=0.1,\quad
m=1,\quad\sigma_{m}=0.1,\quad\mu_{c}=200,\quad\mu_{e}=200,\quad\lambda_{i}=\lambda_{\alpha}=10^{-10},$
(94)

with $M=S=256$ resources and species. The parameters $\sigma_{c}$ and
$\sigma_{e}$ are taken to be the same: $\sigma\equiv\sigma_{c}=\sigma_{e}$.
The dashed gray line represents a slice through the parameter space with the
above parameters fixed and $\sigma=3.5$ where $\rho$ varies from $0$ to $1$.
The gray star represents the parameters $\sigma=3.5$ and $\rho=0.5$; the gray
circle represents the parameters $\sigma=3.5$ and $\rho=0.9$. The values of
parameters $\mu_{c},\mu_{e},\sigma_{c},\sigma_{e}$ are chosen so that when
$M=S=256$, $\langle c_{i\alpha}\rangle,\langle e_{i\alpha}\rangle=O(1)$ and
$\mathrm{Var}\left[c_{i\alpha}\right],\mathrm{Var}\left[e_{i\alpha}\right]=O(1)^{2}$.
The initial conditions are $N_{i}(0)=1/S$ and $R_{\alpha}(0)=1/M$. Numerical
integration is performed using a ‘fourth-order, five-stage explicit Runge-
Kutta method with embedded error estimator of Tsitouras’ implemented as Tsit5
in the Julia DifferentialEquations.jl package [28, 30]. The simulations are
run until $5\times 10^{4}$ time units, and a simulation is considered to have
reached steady state if the average absolute value of all derivatives is less
than $10^{-7}$ within the last $0.05\times 10^{4}$ time units.
The instability phase boundary (solid black) is found by solving the cavity
self-consistency equations (see Appendix A and Eq. (72) simultaneously using
nonlinear least squares, as described in Appendix D.3).

#### Figure 4 (Lyapunov dot plot and diverging trajectories) simulations
details

(a) Details on the calculation and analysis of Lyapunov exponents are
discussed in appendix C.1. Parameters correspond to a slice of the phase
diagram shown in Fig. 3(a) along the dashed gray line.
(b) Simulation parameters correspond to the gray star ($\sigma=3.5$,
$\rho=0.5$) in Fig. 3(a). The initial conditions for the trajectory plotted in
(solid) red is the state given after evolving the system for $5\times 10^{3}$
time units from initial conditions, $N_{i}(0)=1/S$, $R_{\alpha}(0)=1/M$. The
initial conditions for the trajectory plotted in (dotted) blue are generated
by selecting 8 species with abundance greater than $10^{-4}$ at random and
perturbing their abundances by a random amount drawn from a uniform
distribution $U([-10^{-4},10^{-4}])$; the highlighted species is not one of
these 8 species. The Lyapunov exponent is computed using the ‘H2’ method of
Geist [34, 35, 36, 37] for the trajectory plotted in red; see Appendix C.1 for
details.

#### Figure A5 (cavity distribution comparison) simulation parameters

$\displaystyle K=1,\quad\sigma_{K}=0.1,\quad
m=1,\quad\sigma_{m}=0.1,\quad\mu_{c}=8,\quad\sigma_{c}=1,\quad\mu_{e}=8,\quad\sigma_{e}=1,\quad\rho=0.9$
(95)
$\displaystyle\phi_{N}=0.318,\;\phi_{R}=0.729,\;\nu=-0.858,\;\chi=0.412,\;\langle
N\rangle=0.105,\;\langle R\rangle=0.114,\;\langle
N^{2}\rangle=0.0575,\;\langle R^{2}\rangle=0.0258.$ (96)

The uninvadable steady state is found using an iterative constrained
optimization method described in [13]. Cavity parameters are found by solving
the self-consistency equations using non-linear least squares as described in
appendix D.3.

#### Figure A10 (phase diagram with uniform sampling) simulation parameters

All parameters are identical to those in Fig. 3(a) except that the random
matrices and vectors are sampled as:

$\displaystyle d_{i\alpha},x_{i\alpha},\delta K_{\alpha},\delta
m_{i}\sim\text{Uniform}([-\sqrt{3},\sqrt{3}]).$ (97)

This choice of distributions ensures that $d_{i\alpha},x_{i\alpha},\delta
K_{\alpha},\delta m_{i}$ are still standard random variables and appropriate
thermodynamic scaling are used to obtain cavity results that are identical to
those in Fig. 3(a).