{'fraction_non_alphanumeric': 0.04437117489871563, 'fraction_numerical': 0.013231617963968623, 'mean_word_length': 4.310104130907426, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Innovation rather than improvement: a solvable high-dimensional model highlights the limitations of scalar fitness

Much of our understanding of ecological and evolutionary mechanisms derives from analysis of low-dimensional models: with few interacting species, or few axes defining"fitness". It is not always clear to what extent the intuition derived from low-dimensional models applies to the complex, high-dimensional reality. For instance, most naturally occurring microbial communities are strikingly diverse, harboring a large number of coexisting species, each of which contributes to shaping the environment of others. Understanding the eco-evolutionary interplay in these systems is an important challenge, and an exciting new domain for statistical physics. Recent work identified a promising new platform for investigating highly diverse ecosystems, based on the classic resource competition model of MacArthur. Here, we describe how the same analytical framework can be used to study evolutionary questions. Our analysis illustrates how, at high dimension, the intuition promoted by a one-dimensional (scalar) notion of fitness becomes structurally incorrect. Specifically, while the low-dimensional picture emphasizes organism cost or efficiency, we exhibit a regime where cost becomes irrelevant for survival, and link this observation to generic properties of high-dimensional geometry.

The image of evolution proceeding via "improvement" is such a convenient metaphor that, although clearly known to be wrong, it still influences our intuition. Mutations make organisms not better, but different; whether such differences are beneficial depends on the context, which itself is subject to change. The issue becomes especially relevant in an ecological context, where the environment of any one organism is defined, among other factors, by the presence and activity of other species, and can therefore change as fast or faster than the timescale of evolutionary processes.
The lack of timescale separation is particularly pronounced in microbial ecology, and the subject of ecoevolutionary interplay is increasingly in the spotlight today in the context of microbiome research [1][2][3][4][5]. The issue is exacerbated by the observation that most naturally occurring microbial communities are highly diverse, harboring a large number of coexisting species [6,7], each of which contributes to shaping the environment of others [8][9][10][11]. The environment is therefore an intrinsically high-dimensional object, and it is well established that large dimensionality of a problem can lead to qualitatively novel effects [12].
Developing a theoretical understanding of ecology and evolution in the high-diversity regime is therefore an important challenge. So far, however the barrier of dimensionality proved difficult to cross. Although the improvement metaphor is very clearly understood to be misleading [13][14][15], and many others], most quantitative results in evolutionary theory have been derived in a simplistic onedimensional picture (the "fitness landscape" of classic population genetics), which inevitably reinforces certain expectations. It is not always clear to what extent the intuition derived from low-dimensional models applies to the complex, high-dimensional reality [8,10,11,[16][17][18][19][20][21].
Recent work identified a promising platform for investigating highly diverse ecosystems using statistical physics of disordered systems [22], based on the classic resource competition model of MacArthur [23][24][25]. Resource depletion is the simplest form of feedback of organisms onto their environment, and the high-diversity limit of this model was recently shown to be analytically tractable [10,11,26]. As the number of resources becomes large, the community was shown to acquire increasing control over the immediate environment experienced by its members. This should clearly have important implications for how evolution would act in such a community. Uncovering such implications is the focus of this work. Using the high-diversity MacArthur model as our platform, we identify some important deviations from the low-dimensional intuition. We then argue that our conclusions are in fact more general that this particular model, namely that the breakdown of the improvement metaphor at high diversity stems directly from the properties of high-dimensional geometry.
We begin by briefly defining the model. We largely follow the notations of [10], but will point out how the model used there for a purely ecological discussion can be applied to study evolutionary questions. Consider a multi-species community in a well-mixed habitat where a single limiting element X exists in N forms ("resources" i ∈ {1 . . . N }). For example, this could be carbon-limited growth of bacteria in a medium supplied with N sugars. Let n µ denote the population size of species µ ∈ {1 . . . S}. Briefly, MacArthur's model of resource competition can be described as a feedback loop: The availability h i of each resource i determines the dynamics of n µ . The changes in species abundance translate into changes in the total demand for resources, denoted T i . This total demand, in turn, depletes resource availability h i .
Assuming for simplicity that resource dynamics are faster than changes in species abundance, we can assume that resource availability h i at any moment quickly equilibrates to reflect the instantaneous demand T i at that same moment. Specifically, we will posit that all organisms are sharing a fixed total influx of resource R i , so that h i (T i ) = R i /T i , with resource supply R i = 1 + δRi √ N . In this small-fluctuation regime, the specific functional form of the relation h i (T i ) is not important as long as its linearization around equilibrium has a decreasing slope [10], which is a natural condition for stability (higher demand should result in lower supply). The variance of δR i over i is a control parameter describing the heterogeneity of resource supply, and is denoted δR 2 . It is worth noting that even the large-fluctuation regime of the MacArthur model can be rendered analytically tractable [11]; however, for our purposes the above model is sufficient.
In these equations, a species µ is characterized by its metabolic strategy {σ µi } and its requirement χ µ for the limiting element X . The population growth rate of species µ is determined by its resource surplus ∆ µ . In the expression for ∆ µ above, the first term is the total harvest of X from all sources, and the second is the requirement an individual must meet to survive; R i is the total influx of resource i. At equilibrium, each species is either absent (n µ = 0), or its resource intake and expenditure are balanced (∆ µ = 0). A stable (non-invadeable) equilibrium is characterized by an extra condition that all the absent species, if introduced, would be driven back to extinction: if n µ = 0, then ∆ µ < 0. The dynamics defined above always has a stable equilibrium, uniquely defined [10,19,23] by the set of competing species.
We have defined the ecological dynamics for a given set of species. To specify an evolutionary process, we now need to describe how this set is constructed and evolves. In our model, we posit that evolution "discovers" random new species one by one, and these are added to the pool of competitors. Each new species has a random strategy vector σ * , which we take to be binary for simplicity (each entry σ * i is 1 with probability p and 0 otherwise), and a random cost χ * = i σ * i + x * with a normally distributed x * . This cost model ensures a level playing field for all species; see Refs. [10,19]. We assume that new species are generated sufficiently slowly that the ecosystem has time to equilibrate before a new species is introduced. The "evolutionary sequence" we consider is the sequence of non-invadeable ecological equilibria resulting from the competition of all species discovered up to that moment. In other words, we consider the increasing number of species in the pool S as a proxy for evolutionary time. It is convenient to normalize S by N , defining a parameter α = S/N . We will say that a species is "alive (dead) at time α" if it is present (absent) at the respective equilibrium.
To summarize, the control parameters of the model include N (the number of resources), α (the size of the species pool), (the width of the cost distribution), p (the sparsity of the metabolic strategy vectors), and the heterogeneity of resource supply δR 2 . All simulations below use = 10 4 , p = 0.5, and are performed as described in Ref. [10]. Parameters N and δR 2 are specified in the respective legends, and α is taken as the variable against which all results are plotted. A Matlab script reproducing all figures is provided as Supplementary file 1.
It is worth stressing the simplifying assumptions made above. First, we take each new species to be fully random, rather than a small modification of an existing one. Second, the evolutionary process we consider proceeds through a sequence of equilibria. These assumptions make the model analytically tractable while preserving the main feature of interest, namely the eco-evolutionary feedback in a high-dimensional environment, and so provide a reasonable starting point for investigation.
Our model has a convenient geometric interpretation, first introduced by Tilman [24] and illustrated in Fig. 1a. Following Ref. [10], we can think of the metabolic strategies {σ µi } as S = αN vectors in the N -dimensional space of resource availability. Below each hyperplane h · σ µ = χ µ lies the half-space where available resources are insufficient to support species µ. The intersection of such regions over all competing strategies { σ µ , χ µ } defines the "unsustainable region", which we denote Ω: When resource availability h is inside Ω, no species can harvest enough resources to sustain its population. Outside Ω, at least one species can increase its abundance. Therefore, the equilibrium state can only be located at the boundary of Ω, and can in fact be found by solving an optimization problem over this region [10].  1. (a) The geometric interpretation of MacArthur's model. In the space of resource availability, a given set of species defines a convex region Ω where no species can survive. The dynamical equilibrium is always located at the boundary of Ω, and successful invaders slice off regions of Ω. In this cartoon, N = 2 and a community of two specialists (1, 0) and (0, 1) is invaded first by one generalist (0.7, 0.3), and then by another (0.4, 0.6), with progressively lower costs; arrows connect successive equilibria. (For the purposes of illustration, here we allow arbitrary metabolic strategies, without enforcing the constraint σµi ∈ {0, 1}, as in the text.) This geometric intuition suggests that evolution proceeds towards stronger depletion of resources, and consists of organisms with progressively smaller costs, which can then be interpreted as higher efficiency. One might expect such evolutionary process to be qualitatively similar to the simple "best N species" model, shown in panel (b). The "best N species" model ignores any metabolic considerations; instead, each newly generated species is simply assigned a "fitness" value, and the community consists of the N species with highest fitness. The panel shows species presence (in black; each species is a row) for 50 invasion events in one random realization of such a model with N = 10. Each species is present for a period of time, until outcompeted by a more efficient species. This geometric picture proved highly influential, and was used to study competition for N = 1 and N = 2 resources in great detail [24, and others]. It also provides a clear intuition for the evolutionary sequence in our model: as new species are generated, some will have a cost low enough to allow them to invade the community, slicing another piece off the unsustainable region (Fig. 1a). As a result, resource depletion becomes progressively stronger, and community composition shifts towards progressively more efficient species, where efficiency is measured as the species' cost per pathway. In fact, the dynamics (1) has a Lyapunov function, and each successive equilibrium corresponds to a lower value of i R i log h i , precluding any rock-paper-scissors scenarios, and reinforcing the expectation of a "linear" progression [19].
This intuitive picture emphasizes the role played by the individual species' efficiency, and suggests that the evolutionary sequence should be qualitatively similar to the "best N species" model illustrated in Fig. 1b. This is a model where every species is described by a fitness value, and at each point in time, the community consists of the best N species "discovered" to date. In other words, this model is a natural generalization of "survival of the fittest" to the case where N resources allow coexistence of up to N species. In the "best N " model, each species enjoys a period of existence, until it is outcompeted by someone better.
However, the expectations set up by a low-dimensional picture often prove incorrect at high dimension. Fig. 1c shows an evolutionary sequence for one simulation of the MacArthur model at N = 15. The striking qualitative difference between Fig. 1b and Fig. 1c is that the supposedly "outcompeted" species keep returning to again become part of the equilibrium state.
In our model, species are never removed from the pool of competitors. Although a particular species may be driven to zero abundance at a certain equilibrium state, it never goes completely extinct: biologically speaking, we assume that each previously discovered species is preserved in some other spatial patch, or in a dormant form (a spore). Thus, the model is set up to allow a previously outcompeted species to keep trying its luck again at subsequent times. Nevertheless, such "returns from the dead" appear surprising.
A clear qualitative phenomenon that contradicts the naive intuition is a good window into how evolution might be acting differently at high diversity. We therefore start by characterizing this effect within our model: for a species that was alive at time α 1 but dead at α 2 , what is the probability of its return at time α 3 ?

II. THE REPLICA-THEORETIC CALCULATION
MacArthur model is special in that it is a global optimization problem [23], and this makes it analytically tractable using methods of statistical physics of disor-dered systems [10,11,22]. The order parameters are: Both m and q can be explicitly computed as functions of the control parameters of the model.
The surviving species are those whose resource surplus ∆ is exactly zero. For any one α, the probability distribution p α (∆) was also computed in Ref. [10] and is shown in Fig. 2a. Its Gaussian part is centered at −pm with variance p(1 − p)q + 2 , and the delta-peak corresponds to the survivors.
To understand evolutionary sequences in our model, we need to know what happens if new species are added to the system. We know the new number of survivors, but how many of these were already present and stayed on, and how many are newcomers? This information is encoded in the joint probability p α1...α k (∆ (1) . . . ∆ (k) ), computed for a typical evolutionary sequence. In particular, the return probability of a previously outcompeted species (alive at α 1 , dead at α 2 , alive again at α 3 ) is encoded is the three-point joint probability p α1α2α3 .
Before tackling the 3-point problem, let us begin by computing the 2-point joint probability p αα (∆, ∆ ). To do so, consider two copies of a system, where the first has αN species, while the second has those same species plus (α − α)N extra ones. The calculation is similar to that of Ref. [10], except in addition to m, q and m , q characterizing each of the two copies (and satisfying the same equations as in Ref. [10] for respectively α and α ), there is now an additional order parameter describing the coupling: The replica calculation is relatively straightforward and is detailed in Appendix A. We find that the 2-point distribution p αα is a double Gaussian restricted to the negative quadrant (Fig. 2b), with excess probability accumulating at the quadrant boundaries (corresponding to species alive only at α or only at α ) and at the origin (corresponding to species that are alive at both α and α ). Just as the variance of p α (∆) was given by p(1 − p)q + 2 , we find that the correlation matrix of the Gaussian in p αα (∆, ∆ ) is given by: The truncated double-Gaussian form of p αα (∆, ∆ ) was expected, since the two marginals of this distribution p α (∆) and p α (∆ ) must take the form shown in Fig. 2a.
The key parameter to determine from the replica calculation is r. For a homogeneous resource supply (δR 2 = 0, i.e. all R i are identical), the equation we find takes an especially simple form: Here ω, ω are auxiliary noise variables, and D(ω, ω ) is a 2-dimensional Gaussian measure with the correlation matrix C above: We use the notation C(r) to stress that the unknown r enters on both sides of the equation (3). For a given set of parameters, this equation can be solved numerically.
The general equation for a heterogeneous resource supply (δR 2 = 0) and its derivation is provided in Appendix A.
We now return to the problem of computing the complete multi-point joint distribution. Thanks to the mean-field nature of our model, and similarly to the 2dimensional case, the 3-point joint distribution p α1α2α3 can only take the form of a 3-dimensional Gaussian restricted to the negative octant ∆ (1),(2),(3) < 0 (Fig. 2c). Now, however, this distribution is entirely determined by its projections, with no free parameters. Our 2-point calculation above is therefore sufficient to write the expression for the 3-point (or indeed any n-point joint distribution). Specifically, inside the octant the probability distribution is a 3-dimensional Gaussian with the correlation matrix: and excess probability accumulates at the octant boundary. In this expression, the quantities q 1,2,3 are known from the single-copy calculation [10], and r ij are found as solutions to the appropriate version of the two-copy equation (3).
Knowing the joint distribution p α1α2α3 (∆ 1 , ∆ 2 , ∆ 3 ), the "return probability" is easily computed as the conditional probability: For α 1 < α 2 < α 3 , this is the probability that a previously outcompeted species (alive at α 1 , dead at α 2 ) will again be present at equilibrium at a later time α 3 .

III. THE RETURN PROBABILITY AS A WINDOW INTO EVOLUTIONARY PROCESS
is a Gaussian whose weight over the positive semi-axis is collected into a delta-peak corresponding to the survivors [10]. (b) For k = 2, the 2-point distribution is a 2-dimensional Gaussian restricted to the negative quadrant, with excess probability accumulating at the quadrant boundaries. The projection of this distribution on any of the two axes takes the shape shown in (a), as it must. Knowing the marginals leaves one free parameter to be determined through a replica-theoretic calculation, corresponding to the correlation of ∆ (1) and ∆ (2) . (c) For k = 3, the shape of the distribution is again a multi-dimensional Gaussian restricted to negative ∆ (1),(2), (3) . Importantly, for all k ≥ 3 such a distribution is fully determined by its known 2-dimensional projections (gray lines).
For comparison, the panel also shows the 2-point conditional probability p(alive at α 3 | dead at α 2 ); this corresponds to not specifying any information at α 1 , and so this curve is labeled "α 1 = ". The analytical results are in good agreement with simulations (Fig. 3a). The observed deviations are a small-N effect of simulations. Although we used a value of N as large as computationally feasible (N = 200), estimating many-point correlation functions from simulations is especially demanding in sample size; but note the excellent match with the α 1 = curve, which is a two-point function.
This agreement allows us to use our analytical results to investigate the large-argument behavior of p return (α 3 | α 1 , α 2 ). Information about a species' past modifies the expected probability of its survival: knowing it was present or absent a particular time in the past makes our chances to observe it at a future time, respectively, higher or lower. Any consistent "selection pressure" (when species favored in the past continue to be favored in the future) should create persistent memory effects. Thus, the long-term behavior of p return provides important insight into how evolution proceeds in this model.
When α 3 is large, p return goes to zero and it is convenient to multiply this probability by α 3 . Indeed, the number of survivors at equilibrium is of order N (in fact, asymptotes to exactly N ; see [10]), while the number of species in the pool is α 3 N . The probability of a randomly drawn species to be present at equilibrium is therefore 1/α 3 . At large α 3 , this is a natural baseline against which the return probability should be evaluated. The ratio of p return to 1/α 3 quantifies the "effect of knowing a species' past", and is shown in Fig. 3b.
Consider the dependence of these curves on parameter α 1 . Qualitatively, the behavior is easy to understand. For fixed {α 2 , α 3 }, the return probability increases with α 1 , and the α 1 = curve is the lowest of all (Fig. 3a,b). This is the expected behavior: knowing that a species was previously alive increases its chances of being found alive again, and the effect is stronger if the information is recent.
Quantitatively, the behavior is much more interesting. Keeping α 2 = 10, let us fix α 3 to some large value, for instance α 3 = 10 3 , and consider p return as a function of α 1 (the transect in Fig. 3b). We find that the return probability exhibits a marked transition in behavior: knowing that a species was present at α 1 can greatly increase its likelihood of survival at α 3 -but only if α 1 exceeds a critical value α 1 > α C (Fig. 3b, inset). This behavior stems from the V/S phase transition described previously in this model, and indicates that the "selection pressure" is very different in the two phases. Indeed, one can show that for α below α C , a species' cost has no bearing whatsoever on its survival, in stark contrast with the intuition . Also shown is the 2-point conditional probability of a species absent at α2 to appear at α3 (dotted line labeled α1 = ). (b) Same as (a), normalized by 1/α3 (the asymptotic probability of a randomly drawn species to be present at α3) to show long-term behavior. Knowing that a species was absent at α2 makes it less likely to be present at α3 than a randomly drawn species (dotted line, always below 1). Additional information of being alive at α1 enhances the chances [colored lines as in (a)]. For a large α3 (transect) the enhancement remains substantial only if α1 exceeds a critical value αC ≈ 3.7 (inset; black line is the theoretical prediction). This critical value corresponds to the V/S phase transition in the model [10], and signals a change in how selection pressure acts in the two phases. (c) The correlation of species cost x with resource surplus ∆ (and therefore survival) as a function of α. Shown are three curves for different heterogeneity of resource supply δR 2 , a parameter that shifts the location of the phase transition (the critical alpha αC = 2.0, 2.7 and 3.7, respectively). For α < αC , cost and resource surplus are uncorrelated. In this phase cost is irrelevant for survival, in stark contrast to the intuition of Fig. 1. of Fig. 1a. This surprising behavior is demonstrated in Fig. 3c which shows the correlation coefficient between a species' cost x µ and its resource surplus at equilibrium ∆ µ . The theoretical curves (derived in Appendix B) predict zero correlation for α < α C , again in excellent agreement with simulations. Knowing that a species was alive at α 1 < α C tells us nothing about its cost, and so does not modify its likelihood of survival at a large α 3 , explaining the behavior observed in Fig. 3b (inset).
As we can see, the phase transition exhibited by our model endows it with a rich behavior. Although highly interesting, this transition is specific to this particular model, and we will return to its detailed exploration elsewhere. Here, we will focus instead on the general lesson, namely the failure of the intuition suggested by the lowdimensional picture of Fig. 1. In the final section, we will trace the origin of this failure to a generic property of high-dimensional geometry.

IV. IMPROVEMENT VS. INNOVATION
What exactly was wrong with the intuition of Fig. 1a? The low-dimensional picture places significant emphasis on cost. However, a simple argument illustrates why this emphasis is misleading at high dimension.
Consider a pool of αN species, whose strategy vec-tors σ are drawn uniformly from the (N − 1)-dimensional sphere, all with cost exactly 1. Together, they define a certain unsustainable region Ω. For any other species, let us call it viable if there exists a point in Ω that can support its growth. Viability of a species means that there exist circumstances (a vector of resource supply R * ) under which it would be able to invade. For a given, externally fixed resource supply R, viability is required, but not sufficient for invasion.
Drawing a random new strategy σ * uniformly from the same sphere, what cost χ should we give it to ensure a 90% probability that this species is viable? Specifically, how high can we go above χ = 1?
For N = 2, the unsustainable region lies almost entirely within the unit circle (Fig. 4a), and the maximum tolerated δχ ≡ χ − 1 is exceedingly small (Fig. 4b). However, as dimensionality increases, the relative volume of the shaded regions in Fig. 4a explodes, and the cost of the strategy becomes essentially irrelevant for viability. As an example, let us take α = 10. At N = 2 ensuring 90% viability requires a cost of at most 1.00008, so cost is extremely important: any successful invader is required to be as efficient as the existing community members. But at N = 20, the threshold becomes 1.25, and so any new species with a cost close to 1 is viable. In this regime, the cost has virtually no role in determining invasion success: 4. (a) The unsustainable region Ω defined by species with strategies drawn from a unit sphere, all of cost 1. For N = 2, most of the region Ω lies within the unit circle, and any successful invader must have a cost very close to 1. However, as dimensionality increases, the relative volume of the shaded "corners" lying outside of the circle grows exponentially, and higher cost differences can be tolerated. (b) The tolerated cost difference δχ ≡ χ − 1 for a random new strategy to be "viable" (see text), i.e. to intersect the unsustainable region Ω defined by αN random strategies of cost 1 as in (a). For each N and α, 1000 trials were simulated (sampling 100 random invaders, in 10 independent replicates); shown is the δχ sufficient in 90% of the cases. Error bars not shown to reduce clutter. At N = 2 the tolerated δχ is very small, but quickly becomes of order 1 as the dimension increases. (c) The two paths to invade a community.
whether a species can invade the community depends exclusively on its metabolic strategy.
To put this in intuitive terms, there are two paths to invade a given community (Fig. 4c). One is to do what someone is already doing, but better. This corresponds to slicing off a region of Ω by a cut parallel to an existing plane, and we call this path "improvement". The alternative is to slice off one of the many corners of Ω along a new plane, introducing a strategy distinct from any existing ones (and quite possibly more costly): a path we might call "innovation". In dimension 1 (the classic onedimensional fitness), the innovation path does not exist, but in high dimension it can easily become the dominant mode of invasion, as we have seen in our model.

V. DISCUSSION
In this work, we have used an exactly solvable model with eco-evolutionary feedback to demonstrate how the intuition derived from a low-dimensional picture fails at high dimension (in our model, the dimension is the number N of resources for which organisms are competing). To the extent that metaphors are useful, we propose that the high-dimensional evolutionary process considered here can be better described through the metaphor of innovation, rather than improvement. Borrowed from the business literature, the term "innovation" also seems to carry the connotation of enhancement, but it is often stressed that innovation is rarely focused on cost reduction, but rather on finding new markets (niches), and is often costly, but necessary for survival of an enterprise [27]. Importantly, innovation is also distinct from invention: while some evolutionary "discoveries" result in a qualitative change of lifestyle or physiology (e.g., the evolution of photosynthesis or flight), most need not have that character. As far as metaphors go, the analogy seems appropriate.
In real life, the reasons for the improvement metaphor to be violated are plentiful, and frequently cited. For instance, coevolution of predator and prey can proceed in circles, prey continuously changing its strategy to avoid the predator, and the predator adapting. Other mechanisms include non-transitivity (rock-paper-scissors scenarios; see e.g. Ref. [28] in this issue), or non-adaptive evolutionary mechanisms such as inherent stochasticity, or hitchhiking mutations. For all these reasons, as is often discussed, it would certainly be naive to expect real ecological or evolutionary dynamics to ever be a simple gradient ascent. It is therefore important to stress once again that all the phenomena described in this work were studied in a model whose global dynamics does constitute a gradient ascent, making the initial observation of Fig. 1 all the more surprising. Our analysis highlights that in high dimension, the intuition promoted by a onedimensional notion of fitness is structurally incorrect even if none of the additional mechanisms are at play.
Appendix A: Replica calculation of the 2-point joint distribution Consider a system of α 1 N species. We would like to understand what will happen when the number of species is increased to α 2 N with α 2 > α 1 . To do so, we formally consider two ecosystems, one with α 1 N species, and another with α 2 N , the first α 1 N of which are identical to the species in the ecosystem #1. Denoting Z 1 and Z 2 the partition functions for each of the two copies, we seek to compute log(Z 1 Z 2 ) xµ,σµi , where the angular brackets denote the averaging over "disorder", namely the particular realization of the organisms' costs x µ and metabolic strategies σ µi . Using the replica trick: we replace this problem by the (much easier) task of computing Z n 1 Z n 2 xµ,σµi , which we will do for an integer n, but then reinterpret n as being a real number one can send to zero. The validity of this procedure, as always, is justified a posteriori by the agreement with simulations.
The calculation we are about to undertake is a more complex (two-copy) version of the calculation performed in [10] and we will therefore keep the notations consistent with that computation, continuing from section S7.1 of the supplementary material to that paper. Now that we formally have two coupled ecosystems, we denote the resource availability in each as h ik and the resource surplus of each species as ∆ µk , where the index k takes two values, k = 1 or 2.
Since the cost per pathway of all species is approximately 1, we expect h ik ≈ 1 at equilibrium, at it is therefore convenient to work with the shifted variables g i ≡ 1 − h i . Note that this definition corrects the unfortunate choice of scaling made in Ref. [10], which introduced an extra factor of 1/N (Ref. [10], section S7.1), only to later remove it with another change of variables (section S7.6, ibid ). Except for this modification, the other notations have been kept consistent.

Splitting the partition function into parts: Ai, B1 and B2
Just as in Ref. [10], we start from the observation that our ecosystem adopts a state that maximizes a certain simple function F over a complicated region (which is exactly the unsustainable region Ω introduced in the main text). The function being maximized has the form i F i (g i ), and encodes the choice of the resource supply model, which could in principle be kept fairly general. The choice made in the main text corresponds to setting F i (g i ) = −R i g i + 1 2 g 2 i (compare to section S5.3, ibid ), but for now we will keep the more general (and, more importantly, concise) notation F i (g i ).
With these notations, the individual partition functions for the two systems are given by: Introducing replicas, and considering separately the species shared by the two systems, and those unique to the second ecosystem: Performing the averaging over disorder is no different than in the single-copy calculation [10], and we find: a all valid (µ, k) Our next step is to decouple indices i and µ by introducing new variables: Here, q 1 , q 2 as well as m 1 , m 2 are quantities that refer only to a single copy; when it comes to solving saddle-point equations, we will be able to simply substitute the (already established [10]) single-copy values expected for the corresponding α. Note that unlike q 1 , q 2 , the overlap matrix r ab is not manifestly symmetric in its indices.
Introducing new variables requires inserting a corresponding delta-function: We can now factor our integral as follows: with A i , B 1 and B 2 given by: At a fully replica-symmetric saddle point, we must have m a k = m * k (similarly form a k ). For two-index quantities we distinguish the diagonal and off-diagonal parts, e.g.: This yields: Introduce rescaled notations: These are similar to the ones used in [10], except with new variables (r, ρ, γ and δ) added to describe the coupling between our two ecosystems. Note that the somewhat non-orthodox scaling of some of the variables is retained to preserve consistency with Ref. [10]. In these notations, and taking the limit n → 0: The calculation of A i and B 2 in this expression is a straightforward generalization of the one performed in [10], and we will do this first (sections A 2, A 3). The final ingredient, namely the expression for B 1 , will take slightly more effort, and will be the subject of sections A 4 through A 6.

a. Symmetric resources
First, let us assume for simplicity that the supply of all resources is identical R i ≡R = 1. We then have A i ≡ A: Recall now the expression for the resource supply function F (g) = −g − 1 2 g 2 to find: We now write (leaving the subscript k implicit): Substituting this into the expression for A: In rescaled notations: In this expression, the first terms have the desired form a (. . . ) (uncoupled replicas). To deal with the last two, we use Feynman's trick of introducing extra Gaussian variables: (the curly D denotes the standard Gaussian measure with variance 1). This lets us write: As in Ref. [10], we shift the variablem k by settingm k ≡ −1 + δm k √ N (note that the "1" in this substitution is actually the average resource supplyR, which we set to 1). This scaling ansatz will later be verified by the extremum condition for δm.
where g ≡ (g 1 , g 2 ) and log Dz x n = log Dz (1 + n log x + . . . ) = log 1 + n Dz log x + . . . = n Dz log x + . . . Therefore: This is a Gaussian integral of a quadratic form, and is easily computed: As in the single-copy calculation [10], the extremum conditions with respect to δm k set δm k = 0, and the expression for the partition function simplifies to: The variables a k , b k , γ, δ do not appear in B 1 or B 2 . Consequently, the extremum conditions for these variables can already be computed, and these variables eliminated: The partition function is now a function of only m k , q k , x k , r, and ρ: In the interest of clarity, the calculation above assumed fully symmetric resources. Restoring full generality renders the intermediate expressions slightly more complicated, but is fairly straightforward, so we summarize it briefly. The exact same steps lead us to the following expression: As a sanity check, note that setting R i ≡ 1, with δR i ≡ 0 yields the expression for the homogeneous case. Taking the Gaussian integral is again straightforward, at which point we recall that what we need is not A i individually, but the sum i log Ai n . After summation over i, the term linear in δm cancels to leading order in N ; the expansion starts with the quadratic term. The extremum condition will therefore set δm k = 0 as before. Also to leading order in N , in the expression for M i the difference between R i andR = 1 is negligible. Integrating away a k , b k , γ, δ as before, we obtain our final expression, which now includes two new terms, each proportional to δR 2 :

Computing log B2
Recall the definition of B 2 : This is a single-copy expression: it pertains to the species that are present only in one of the ecosystems. Not surprisingly, therefore, this B 2 is an exact match to an expression already computed in Ref. [10], where it was called simply B with no subscripts. We can therefore write the result directly: where ψ 2 ≡ p(1 − p)q 2 + 2 , λ 2 ≡ pm 2 /ψ 2 and I(λ) is given by:

Computing log B1 (and introducing X )
So far our calculation was a straightforward generalization of Ref. [10]. The only thing that remains is B 1 , and this is the part that requires some effort: Under replica-symmetric assumptions: By analogy with ψ ≡ p(1 − p)q + 2 , introduce a yet another notation: The exponent in B 1 becomes: Decoupling replicas by introducing auxiliary Gaussian fields: where ∆ ≡ (∆, ∆ ) (both with and without the "hats"), and As a result, and neglecting additive terms of lower order: To avoid carrying around the overall multiplication factor, let's rewrite the partition function a little bit, and state this as an intermediate result: The X in the last term is given by: and it would perhaps be helpful to remind the reader of our notations: The full expression for X is rather complicated, but we will now make an important observation that will simplify our life dramatically.

Simplifying observation: ρ = 0
Recall the original definition of r ab : Note, however, that the numbering of the two sets of replicas (for one ecosystem and for the other) is arbitrary, and we can reorder one set without changing the other. It follows that at the saddle point, all the entries in this matrix must be the same (this is the parameter we called r). In particular, there can be no difference between the diagonal and the off-diagonal entries, and therefore the variable we called ρ ≡ β N (r D − r O ) must vanish at the saddle point. This observation means that we only need to calculate X up to linear order in ρ. Indeed, to determine r all we need is the saddle-point condition for ρ, which is computed at ρ = 0: log(Z 1 Z 2 ) = 0 ⇒ will determine r.
6. Final ingredient: computing X to linear order in ρ Let us rotate the basis of the Gaussian noise variables to introduce: (note that these are not orthogonal, and therefore these Gaussian variables are no longer independent). We choose the third basis vector to be of unit length, and orthogonal to both ξ 1 and ξ 2 ; its explicit expression need not be specified. Since the integrand only involves ξ 1 and ξ 2 , this third Gaussian noise variable can be integrated out, and we are left with: Here ξ ≡ ξ 1 ξ 2 and the integration measureD ξ is that of two correlated Gaussian variables: with the correlation matrix: , where ψ k = p(1 − p)q k + 2 as before.
Minimizing a quadratic form is not hard; the only complication is the negativity constraint on ∆ 1 , ∆ 2 . Globally, the global minimum is always zero, but depending on the values of ξ 1 , ξ 2 , the point where it is achieved may lie outside of the allowed quadrant, in which case we will have to content ourselves with the smallest value at the quadrant boundary (where either ∆ 1 or ∆ 2 is zero). We must therefore consider several cases, and the expression for the integrand will be different. Specifically, we find that the integration plane splits into 4 regions we label A, B, C and D: Integrand: Integrand: 0 Region C: Region D: Integrand: Integrand: To get some reassurance we are on the right track, note that setting ρ = 0 decouples the integrals over ξ 1 and ξ 2 , and we immediately find: This is precisely the expected result. Indeed, substituting this into our expression (A4) yields which is reassuring, as we of course expect log Z 1 Z 2 = log Z 1 + log Z 1 .
To determine r, we now need to calculate the first-order correction to (A5) and find the relevant saddle-point (. . . ) = 0. A non-zero ρ brings two differences: first, the integration region is not quite the same, and second, there is an extra term in the integrand. Let us write this as follows: (Here we add and subtract integration regions in the obvious sense of adding/subtracting the integrals taken over them.) The terms in the last line all vanish at ρ = 0: the first two, because the integration region vanishes, and the final one, because it is explicitly multiplied by ρ. Conveniently, therefore, only the very last term contributes to the derivative ∂ ∂ρ ρ=0 , and we find: Recalling Eq. (A4), we write the saddle-point equation: Changing the sign of the integration variable, we find the equation for r: at the saddle point (see Eq. (S17) in Ref. [10]). The version quoted in the main text assumed δR 2 = 0 for simplicity. 7. The joint distribution p(∆1, ∆2) Having computed r, we will now demonstrate how this parameter is related to the shape of the joint distribution p(∆ 1 , ∆ 2 ) (for species that are common to both ecosystems). For this, we go back to the expression for B 1 given in Eq. (A3). Instead of directly sending n → 0, we separate out one of the replicas, leaving n − 1 remaining ones, and only then send n to zero (for a more detailed explanation of this trick, see Ref. [10], section S9). We find: The key observation is that the conditional distribution p(∆ 1 , ∆ 2 |w 1 , w 2 , u) is manifestly a properly normalized probability distribution on the negative quadrant ∆ 1,2 < 0. For any {w 1 , w 2 , u}, the limit of p(∆ 1 , ∆ 2 |w 1 , w 2 , u) as β → ∞ is therefore easy to compute: it is always a (properly normalized) delta-function, concentrated where the quadratic form in the exponent reaches its minimum, which could be either within the quadrant, or at its boundary. Let us determine the portion of p(∆ 1 , ∆ 2 ) for purely negative ∆ 1 , ∆ 2 . As β → ∞, the conditional distribution p(∆ 1 , ∆ 2 |w 1 , w 2 , u) concentrates into a 2-dimensional delta-function: p(∆ 1 , ∆ 2 |w 1 , w 2 , u) β→∞ −→ δ(∆ 1 + pm 1 + w 1 φ 1 + uζ) δ(∆ 2 + pm 2 + w 2 φ 2 + uζ) Therefore, for strictly negative ∆ 1 , ∆ 2 : p(∆ 1 , ∆ 2 ) = Dw 1 Dw 2 Du δ(∆ 1 + pm 1 + w 1 φ 1 + uζ) δ(∆ 2 + pm 2 + w 2 φ 2 + uζ).
Introduce new noise variables ξ 1 = w 1 φ 1 + uζ and ξ 2 = w 2 φ 2 + uζ. The noise direction orthogonal to both can be integrated away and gives 1. What remains is a 2-dimensional integral: p(∆ 1 , ∆ 2 ) = Dξ 1 Dξ 2 δ(∆ 1 + pm 1 + ξ 1 ) δ(∆ 2 + pm 2 + ξ 2 ), where the measure is that of two correlated Gaussian variables: (C depends on r through ζ.) The result, obviously, is a double Gaussian with the exact same correlation structure, for variables ∆ 1,2 + pm 1,2 : This expression is valid only for strictly negative ∆ 1,2 . We find that, as claimed in the main text, the variable r indeed controls the degree of correlation between ∆ 1 and ∆ 2 (the resource surplus of the same species between the two ecosystems). Following the same procedure, one can verify that the probability weight at the quadrant boundary is precisely the weight of the same double-Gaussian expression integrated over the forbidden (positive) range of ∆, also as claimed. For instance: The simplest way to prove this is to observe that the distribution we computed must have the marginal computed in Ref. [10], and shown in Fig. 2(a).
Grouping together the terms with∆ and taking that integral: Omitting the constants that cancel out, recalling that q D − q O = N β x, and rescaling the integration variable by a factor p(1 − p)x, we can write: where the numerator and the denominator are given by: Here ∆ 0 ≡ pm+w √ p(1−p)q+ w √ p(1−p)x . If ∆ 0 is positive, then as β → ∞, the exponent tends to a delta function (centered at −∆ 0 ), and the ratio tends to −w ∆ 0 p(1 − p)x. If, however, ∆ 0 is negative, then the exponent starts concentrating at ∆ = 0, where the numerator vanishes. As a result: x∆ = Dw Dw θ pm + w p(1 − p)q + w p(1 − p)x × (−w )(pm + w p(1 − p)q + w ).
Denoting F (u) ≡ u θ(u) and integrating by parts (note that w Dw is a full derivative): Rotating the (w, w ) basis to integrate away the one of the noise directions, we are left with a 1-d integral in which we recognize the function E(λ) ≡ 1 2 erfc λ √ 2 , concluding with a very simple formula: The parameter combination λ ≡ pm/ψ plays an important role throughout our model; for instance, the number of survivors is αE(λ) [10]. Here we see that it also controls the correlation of species' cost and resource surplus (and thus their survival as well). To get the correlation coefficient as plotted in Fig. 3(c), we should normalize this result by the standard deviations. The standard deviation of x is 1. As for ∆, its distribution is shown in Fig. 3(a); the standard deviation of the Gaussian part is ψ. Neglecting the effect of the delta-shaped tail (whose total weight is at most 1/α, and whose effect on the variance quickly becomes negligible as α increases), we find: corr(x, ∆) ≡ x∆ Consider the behavior of this expression as → 0; this is the limit where V/S becomes a true phase transition. As established in Ref. [10], for α below a critical value, ψ remains of order 1, and we find that the correlation vanishes in the V phase, as promised. For α above the critical value, ψ was shown to go to zero linearly with , so that their ratio is of order 1 and the correlation becomes non-trivial. For small, but finite the behavior retains these qualitative features and agrees with simulations, as shown in Fig. 3(c).