{'fraction_non_alphanumeric': 0.03899139225112585, 'fraction_numerical': 0.007638664563817052, 'mean_word_length': 4.449725587656622, 'pattern_counts': {'":': 1, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A collective phase in resource competition in a highly diverse ecosystem

Organisms shape their own environment, which in turn affects their survival. This feedback becomes especially important for communities containing a large number of species; however, few existing approaches allow studying this regime, except in simulations. Here, we use methods of statistical physics to analytically solve a classic ecological model of resource competition introduced by MacArthur in 1969. We show that the non-intuitive phenomenology of highly diverse ecosystems includes a phase where the environment constructed by the community becomes fully decoupled from the outside world.

Understanding the diversity of life forms on our planet is an age-old question. Recent technological advances uncovered that most habitats harbor hundreds of coexisting "species" (most of which are microbial [1][2][3]), and the problem of understanding such communities is currently at the forefront of medical and environmental sciences [4][5][6]. One of the key obstacles arises from the fact that ecological and evolutionary time scales are generally not separable, giving rise to a coupled "eco-evolutionary dynamics" [7][8][9]. The fitness of an organism depends on its environment, but this environment is not fixed: it includes all other organisms in the community, is shaped by their activity and changes on an ecological time scale. Understanding this feedback has long been recognized as an important question of community ecology [10].
A convenient example of such ecological feedback appears in models of resource competition [11]. The survival of an organism is determined by the availability of resources in its immediate environment. In quantitative theories of evolution (population genetics), we typically think of this environment as being fixed externally, but in an ecological setting an experimentalist can only set the conditions faced by the community as a whole, e.g. the overall influx of resources. The immediate environment of an individual is affected by the activity of all other organisms and is not under our direct control. For example, consider increasing the overall influx of maltose (a sugar) to a multi-species bacterial culture. This could lead to an increase of maltose in the medium, opening the community to invasion by a species that grows well on this sugar. Alternatively, this could enable existing maltose-consuming species to expand in population, driving maltose availability back to the same level, or perhaps even depleting it further. The relation between the resources supplied to the community and the immediate environment seen by individual organisms is non-trivial. Our control extends on the former, but organism survival and therefore community structure are determined by the latter.
The mechanisms by which organisms shape their environment (niche construction theory [12]) have been the subject of much research, both at equilibrium (e.g. resource competition models [11]) and out of equilibrium (e.g. in the study of ecological successions [13]). Perhaps the most progress was achieved in the problem of resource competition in a well-mixed community at equilibrium, introduced 50 years ago by MacArthur [14]. However, the geometric approach developed by Tilman in his classic work [15] allowed him to analyze only the cases with N = 1 and N = 2 resources. It is not clear to what extent the intuition derived from low-dimensional models applies to the high-dimensional case. Recently, a simulation-based study of a modestly larger number of resources (N = 10) exhibited a surprising effect whereby a community interacting with another community would exhibit an effective "cohesion" even in the absence of any cooperative interactions between its members, purely as a consequence of environmental feedback [16]. The number of metabolites at play in a complex microbial community in nature is even larger, of order N 100 [17,18]. It is an intriguing possibility that the phenomenology of highdiversity communities could contain qualitatively novel, non-intuitive regimes. However, few existing approaches allow studying niche construction or eco-evolutionary dynamics for a large number of interacting species, except in simulations.
In this work, we show that MacArthur's classic model of resource competition can be formulated as a problem of statistical physics of a disordered system, and solved analytically in the limit of large N . We observe a phase transition between two qualitatively distinct regimes. In one regime, changes of external conditions propagate to the immediate environment experienced by organisms, as expected. However, in the other regime, the immediate environment of individuals organisms becomes a collective property of the community, unaffected by the outside world. This regime, which only arises at sufficient diversity, documents the emergence of a collective behavior as a consequence of large dimensionality.
In defining our model, we follow Ref. [16], but allow for more generality. 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, the availability h i of each resource i in the immediate environment of individuals will determine the dynamics of n µ . The changes in species abundance will translate into changes in the total demand for resources, denoted T i . This total demand, in turn, will determine the resource availability h i . This feedback loop is the focus of our analysis.
A species is characterized by its requirement χ µ for the limiting element X , and the "metabolic strategy" {σ µi } it employs to try and meet this requirement. We think of σ µi as the investment of species µ into harvesting resource i (e.g., the expression level of the corresponding metabolic pathway). Specifically, for given resource availability {h i }, the population growth rate of species µ is determined by the resource surplus ∆ µ experienced by its individuals: The first term is the total harvest of X from all sources, and the second is the requirement an individual must meet to survive. The proportionality coefficient is not important, since we will only be concerned with the equilibrium state where dnµ dt = 0. Species abundances n µ determine the total resource demand T i ≡ µ n µ σ µi . This demand shapes resource availability h i . In the simplest model [16], organisms could be sharing a fixed total influx of resource R i : h i (T i ) = R i /T i . In his original formulation, MacArthur considered a more complex scenario of dynamical resources with renewal rate r i and maximal availability K i ; this would correspond to setting h i (T i ) = K i 1 − Ti ri , see eq. (3) in Ref. [14]. In the interest of generality, here we will say only that the availability of resource i is a decreasing function of this total demand: h i = H i (T i ), and allow the functions H i (·) to remain arbitrary, and possibly different for each resource.
This model admits a convenient geometric formulation, where we can think of the metabolic strategies {σ µi } as S vectors in the N -dimensional space of resource availability. Each hyperplane h· σ µ = χ µ separates this space into two regions (Fig. 1A). Above this hyperplane, a positive resource surplus allows species µ to multiply. Below this hyperplane (shaded), resources are insufficient to support species µ. The intersection of such regions over all competing strategies { σ µ , χ µ } defines the "unsustainable region" Ω: If 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 Ω, which we denote ∂Ω. The dynamics (1) possesses a Lyapunov function, which is convex and bounded from above, similar to the classic model of MacArthur of which this is a generalization (see SI). As a result, the equilibrium state always exists, is unique and stable, and can be found by solving a convex optimization problem over the region ∂Ω. At this equilibrium, each species is either extinct and cannot invade (n µ = 0, ∆ µ < 0), or is present and its resource balance is met (n µ > 0, ∆ µ = 0). Fig. 1B shows an example at N = 2. Here, a community of two specialists σ 1 = {1, 0} and σ 2 = {0, 1}, both with cost χ 0 , is exposed to a mixed strategy σ 12 = {x, 1 − x} with a cost slightly below χ 0 . The species σ 12 will be able to invade, and depending on resource supply, may coexist with one of the specialists (but not both). The equilibrium will harbor one or two species, corresponding to the equilibrium h being located either at an edge or at a vertex of ∂Ω.
The resource depletion rules H i (·) describe the external conditions: how much of each resource is supplied to the community as a whole. In contrast, h describes the availability of resources in the immediate environment of individuals, which ultimately dictates which species survive. Any set of competing strategies { σ µ , χ µ } defines a unique community equilibrium, and so implements a mapping from external conditions into the actual environment h. Our aim is to characterize the properties of A: If resource availability h lies above the line h· σµ = χµ, the species µ will multiply, depleting resources (arrow). B: Competition between S = 3 species; metabolic strategies indicated by arrows (two specialists and one mixed strategy). The equilibrium h is always located at the boundary (highlighted) of the "unsustainable region" Ω; one or two species may coexist. this mapping.
The geometric intuition described above was first developed by Tilman [15], who exhaustively analyzed the cases N = 1 and N = 2. In higher dimensions, however, the enumeration of co-existence regimes for a given set of strategies, like in Fig. 1B, quickly becomes a combinatorially difficult problem. In this work, we therefore adopt the statistical physics approach, and characterize the expected properties of a typical community, when the competing strategies are drawn out of some ensemble.
Specifically, for each species µ, we first pick its strategy as a random binary vector, where each component σ i µ is 1 with probability p, and 0 otherwise. The parameter p allows us to specify the location of a typical competitor on the specialist-generalist axis. We then draw a random cost χ µ = i σ µi + x µ , where is a parameter (the cost scatter, assumed small), and x µ is a Gaussian random variable of zero mean and unit variance. We set the total number of species to S ≡ αN .
The key simplification that makes the problem tractable analytically is the independence of σ µ and x µ : the strategy and its cost are effectively uncorrelated. This assumption is strong, but far from unreasonable. The species competing for the same resources in real communities differ in evolutionary history, lifestyle, and physiology. Modeling the cumulative effect of these differences as a random contribution to the species' likelihood to succeed is arguably a better null model than claiming that the single factor we explicitly consider (the species' metabolic preference) plays the dominant role in determining its intrinsic performance.
Note that setting h i = 1 satisfies the resource balance of all species within a quantity of order , so this cost model ensures that neither specialists not generalists have an obvious advantage [16]. To characterize the fluctuations of resource availability 1 − h i , we introduce: The resource surplus of a typical species is given by: (the angular brackets denote the mean over µ). Negative for most species, ∆ µ should hit zero for the lucky outliers who survive. We find that the spread of resource surplus values is given by ψ ≡ p(1 − p)q + 2 (see SI). Intuitively, this is because species differ in cost (variance 2 ), and their strategy ({σ µi } with variance p(1 − p)) picks out resources with different availability (variance q). For this reason, rather than using q and m directly, for our order parameters we choose ψ and the ratio λ ≡ pm ψ . Each particular set of competitors constitutes "frozen disorder", and the properties of a typical community can The phase transition at → 0. In the S-phase, above a critical α (dotted line), the fluctuations of internal resource availability ψ vanish, shown here on log scale to highlight the transition. B. The distribution of resource surplus at equilibrium. Black, the theoretical prediction; red, simulation data accumulated over 500 realizations at N = 50 and is shown for extinct species only (see SI for details). C. The number of surviving species at equilibrium as a function of α at δR 2 = 1 (cf. the arrow in panel A). Theoretical prediction (black); mean over 500 simulations at N = 50 (red); the deviation at = 10 −4 is an effect of small N . Standard error of the mean is too small to be visible. Dotted line at critical α; shading labels the two phases.
be computed using methods of statistical physics of disordered systems [19], as detailed in the SI. For simplicity, all the results will be quoted for the simplest supply Ti where each resource is characterized by a single parameter: its total supply R i (see SI for the general case). Our calculation yields explicit equations for the order parameters ψ and λ at equilibrium, in the thermodynamic limit N, S → ∞ at α held constant: Here δR 2 is the variance of resource supply R i , and 2π are auxiliary functions that can be expressed in terms of the error function erf.
To study these equations, consider first the limit → 0. In this limit, the parameter space separates into two phases ( Fig. 2A). One of these corresponds to the solution ψ = 1 − αE(λ) = 0 and will be called the S-phase; the other has ψ = 0 and will be called the V-phase. The critical line (dotted line in Fig. 2A) is described by: For δR 2 = 0 the transition occurs at α crit = 2, consistent with the perceptron phase transition [19,20].
To understand the physical meaning of these phases, consider first a community consisting of N perfect specialists with costs χ µ ≡ 1. This community constitutes an example of the S-phase, where the immediate environment of individuals is fully "shielded" from external conditions: faced with an uneven resource supply, species' abundance will adjust to drive resource availability to h i = 1 for all i, restoring symmetry. In general, a restricted set of species (small α) or a strongly heterogeneous resource supply (large δR 2 ) will prevent the community from exactly matching demand to the uneven supply, and the externally imposed asymmetry between resources will propagate into the organisms' actual environment h (the V-phase, "vulnerable" to external perturbations). However, as the community is exposed to new species (α is increased above the critical value; the arrow in Fig. 2A), the community transitions into the shielded phase where the environment h is fully symmetric (m = q = 0) and insensitive to external conditions.
To confirm this interpretation, consider the number of coexisting species at equilibrium. As we have seen, geometrically, this number is the co-dimension (N minus the dimension) of the region of ∂Ω where the equilibrium is located. Remarkably, this elusive quantity can also be computed analytically. Specifically, one can compute the distribution of the resource surplus ∆ of all αN species at equilibrium ( Fig. 2B; see SI): Here θ is the Heaviside function constraining ∆ to be negative. The delta-shaped peak at ∆ = 0 represents the fraction of species whose resource demand is met. The number of survivors is therefore αN E(λ), in excellent agreement with simulations (Fig. 2C). The S-phase where αE(λ) = 1 therefore harbors a complete set of exactly N species. If the perturbation of external conditions is small, no species will go extinct. Since the vectors h i and χ µ (µ running over N surviving species) are related by a full-rank matrix σ µi , this means that the resource availability at the new equilibrium will remain exactly the same, confirming our interpretation of this "shielded" phase. For a non-zero , the strict phase transition is replaced by a crossover (Fig. 3A). The role of in our model is to measure how strongly a species' fate is influenced by intrinsic, rather than environment-dependent (ecological) factors [16]. For large , community structure is no longer shaped by interactions between community members, but becomes dominated by species who outperform others in all circumstances, and the environmental feedback studied here becomes irrelevant. For small , however, the distinct features of the "shielded" and "vulnerable" phases remain clearly recognizable: the fluctuations of resource availability are, respectively, of order or much larger than (Fig. 3B).
This result has intriguing implications. Consider a community facing the strongly uneven resource supply shown in Fig. 3B (top panel). Define a species' individual performance as its growth rate when placed in this environment, with no other organisms present. One might expect this performance metric to be predictive of species' survival in a community setting: surely, increasing the supply of maltose to a community should favor organisms that grow well on maltose. In the more intuitive V-phase, this expectation is indeed correct. However, in the S-phase the internal environment becomes a collective property governed by the statistical properties of the species' pool, rather than by the external conditions ( Fig. 3B; bottom panel). As a result, the performance measured in external conditions becomes irrelevant: it no longer predicts whether a species will survive (Fig. S3).
In ecological terms, the model considered here was purely competitive: increasing the abundance of any species reduces the growth rates of everyone else, i.e. there are no "cooperative interactions". Nevertheless, we have shown that at high dimension, the parameter space of this classic resource competition model contains a strongly collective regime.
These conclusions were drawn in the context of a particular, highly simplified model. In particular, our analysis ignored spatial structure, assumed deterministic dynamics, and considered the equilibrium states only. It is clear that natural communities are never in steady state, and stochasticity and spatial structure are tremendously important in most contexts. Nevertheless, the goal of this work was to explore specifically the feedback of organisms onto their environment and identify the implications of large dimensionality. For this purpose, the simplified model adopted here provides a convenient starting point, and highlights the promise of applying statistical physics to gain analytical insight into the non-intuitive phenomenology of large-dimensional networks [21] and highly diverse ecosystems.
We thank Michael P. Brenner, Carl P. Goodrich, Alpha Lee, Emily Zakem and David Zwicker for helpful discussions, the Harvard Center of Mathematical Sciences and Applications, and the Simons Foundation. This work was completed at the Aspen Center for Physics, supported by National Science Foundation grant PHY-1066293.

SUPPLEMENTARY MATERIAL THE LYAPUNOV FUNCTION F
Recall that the dynamics of our model are given by This section will show that this dynamics possess a Lyapunov function: (S1) In other words, F increases on any trajectory of the dynamics above. In addition, we will show that F is convex and bounded from above.

Proposition 1: F increases on any trajectory
We first note that the derivative of F with respect to a species' abundance n µ is precisely the resource surplus ∆ µ : Therefore, F is indeed a Lyapunov function:

Proposition 2: F is bounded from above
To see this, recall that H i (·) was required to be a decreasing function of its argument; moreover, to forbid unbounded growth of any species, we required that for large enough demand T , the resource availability H i (T ) should go to zero. It follows that its integralĤ i (x) ≡ x H(T ) dT grows sub-linearly; in other words, for any λ > 0 we have We conclude that F ( n) goes to −∞ as the norm of the abundance vector increases (this precisely corresponds to forbidding infinite population growth). A continuous function defined on the positive quadrant {n µ ≥ 0} and going to −∞ at the boundary of this region is bounded from above, as claimed.

Proposition 3: F is convex
To see this, note that for any function f ( n), the following two operations leave its convexity invariant (M is an arbitrary matrix): 1. adding a linear function of its arguments: f ( n) → g( n) = f ( n) + M n; 2. performing a linear transformation of its arguments: Given these observations, convexity of F directly follows from the convexity ofĤ i (x) (which is an integral of a decreasing function).

LOCATING THE COMMUNITY EQUILIBRIUM: THE GEOMETRIC INTUITION
The main text shows that the equilibrium of community dynamics is always located at the boundary of the "unsustainable region" Ω defined in the text. Which boundary point is selected? Here we present an intuitive geometric argument, which will be formalized in the following section.
Let h * be the resource availability at community equilibrium. For concreteness, consider the case N = 2, and assume the equilibrium state harbors two species { σ 1 , χ 1 } and { σ 2 , χ 2 }, so that the point h * is the intersection of lines h · σ 1 = χ 1 and h · σ 2 = χ 2 . Consider now the vector of total demand T * at this equilibrium. By definition, it is a linear combination of the two strategy vectors: T * = n 1 σ 1 + n 2 σ 2 . Importantly, the coefficients here must be positive. We conclude that at equilibrium, the vector T * must point "strictly outward" relative to the region Ω, as in Fig. S1.
This property is sufficient to uniquely determine the equilibrium point. Indeed, consider a vector field T ( h), where to each point of the resource availability space h 0 we associate the vector of total demand T 0 that corresponds to such resource depletion, i.e. such that H( T 0 ) = h 0 . The intuitive argument above suggests that the equilibrium of community dynamics can be found by following this vector field. And indeed, this vector field is a gradient of a certain function, and therefore locating the equilibrium corresponds to maximizing this function. This is formally proven in the following section.
We stress that the vector field T ( h) does not describe the dynamics itself; it is merely a tool to find the equilibrium point. The trajectories of the system in the harvest space are not integral lines of this vector field. S1. The geometric intuition behind the selection of the equilibrium point: at equilibrium, the vector of total demand T * must be pointing "strictly outward" relative to the unsustainable region Ω. Here, at a two-species equilibrium of σ1 and σ2, T * must lie within the sector highlighted in blue. The equilibrium point can therefore be found by following the vector field T ( h) (in gray).

LOCATING THE COMMUNITY EQUILIBRIUM: THE FORMAL PROOF
Proposition 1: There exists a functionF defined on the harvest plane, such that its gradient at any point h is the demand vector T that corresponds to this resource availability vector: . This function has the property that ∂F ∂Ti = H i (T i ). The functionF can be explicitly constructed as the Legendre transform of F : It is easy to check that this function satisfies the desired requirement. Indeed, for each component i (and omitting this index for simplicity): Consider now an equilibrium community C with the total demand T * , and the resource availability vector h * . We already know that h * lies at the boundary of the unsustainable region Ω (see main text). To determine exactly which boundary point is selected, we make the following observation: Proposition 2: For any other vector h 1 ∈ Ω, we have Proof: Since h * is the equilibrium state, we can write: As for h 1 , it lies in the unsustainable region Ω, and therefore: Subtracting the former from the latter, and summing over µ, we conclude: This means that the equilibrium point is such that the value ofF cannot be further increased: any movement within the unsustainable region goes against the gradient field ofF . Corollary: The equilibrium community state corresponds to the maximum ofF over the unsustainable region Ω. It is worth contrasting our approach to other situations where community-level objective functions may appear, e.g. community-level flux balance analysis, or similar approaches. In certain contexts, an optimization-based framework is simply postulated, and serves as an exploratory tool to investigate the possible regimes of network performance: for instance, the total metabolic output of a consortium might be taken as a "plausible" global objective function for the community to optimize. Here, we stress that no community-level objectives are postulated; the fact that the ecological dynamics in this model take the form of a global optimization problem is a "lucky" consequence of explicitly specified dynamics of (purely "selfish") individual species. This special feature makes MacArthur's resource competition model an especially convenient starting point for investigating the consequences of high dimensionality in an ecological context.

THE ALGEBRA OF "PASSING INTO THE HARVEST SPACE"
Above, we have shown that locating the equilibrium of our ecological dynamics is in fact a convex optimization problem in the P -dimensional space of species abundances. We then described how this optimization problem can be formulated directly in the N -dimensional space of "harvests" h i . To build intuition, it is instructive to consider the following algebraic argument showing how the two optimization problems are mapped into each other.
In order to locate the maximum of the Lyapunov function F , we investigate the large-β limit of the partition function Z(β): To compute Z, we first introduce T i as convenient auxiliary variables. This allows integrating over n µ : We now focus on the integral over T i . For large β, it can be computed using saddle-point method. Denoting iθ i ≡ h i , we find that the saddle-point T * i is defined by the condition: which justifies our suggestive notation (we recognize h i as the substrate availability at equilibrium demand T * i ). This condition implicitly defines T * i as a function of h i , so all that remains is the N -dimensional integral over h: TheF in the exponent is precisely the Legendre transform ofĤ i : For large β, this is again a saddle-point integral. We are starting to recognize the problem of extremizingF ; however, here it is computed for purely imaginary arguments, and so a few more steps are needed. The integration contours cannot simply be rotated onto the real axes, since the integrand has a complicated pole structure. Instead, we can convert the integration contours into piecewise-linear shapes, two of which are purely imaginary, and one is purely The deformation of the integration contour is allowed only as long as the poles are not crossed, and the integrand has a pole whenever ∆ µ = 0 (we note that the denominator in (S4) is (−∆ µ ), the negative resource surplus of species µ). Thus in our N -dimensional integral, the shifting of each contour will depend on the exact values of all other variables. Thankfully, the integrand can have an extremum only if all h i are real, and whenever N − 1 variables h i take real values, the remaining one can vary (on its real-valued portion of the contour) from 0 to the highest value that can be reached without crossing any of the hyperplanes ∆ µ = 0. The region delimited by these hyperplanes is precisely the "unsustainable region" Ω defined in the main text. We conclude that for the purposes of the saddle-point calculation, our integral becomes: The certain lack of rigour in our description of the transition from Eq. S3 to Eq. S4 will not be a problem. The purpose of this section is to build additional intuition about the algebraic structure of the problem, and analyzing the expression (S4) will prove instructive. However, the following sections will only use the fact that community equilibrium maximizesF , a result that was rigorously obtained in section titled "Locating the community equilibrium: A formal proof". In expression (S4), the exponential term e −βF dominates the integrand everywhere, except in the immediate vicinity of the region boundary where 1 ∆µ diverges. If β is large, but finite, the extremum is achieved at a point h * lying strictly inside the region Ω, at a distance of order 1/β from the nearest bounding hyperplanes. In this "finite temperature" regime, all species n µ have non-zero abundance: since n µ enters into Z as e −βnµ |∆µ| , the observables n µ follow an exponential distribution with mean n µ = 1 β |∆µ| . In the zero-temperature limit (β = ∞), this expected abundance vanishes for all species except a select few, for which ∆ µ is precisely zero. Thus, as described in the main text, the extremum h * reaches the boundary of Ω. At this value of harvests, a finite set of species have resource surplus of precisely zero, corresponding to finite-abundance survivors. The resource surplus of all other species is negative, and they go extinct at equilibrium.
To make this argument more precise, we note that at large, but finite β, Eq. (S4) gives us For a large β, the sum over µ is dominated by only a few terms, those corresponding to the closest hyperplanes for which ∆ µ tends to zero. Denote their set S (for "survivors"). The extremum condition: This of course makes perfect sense given the definition ofF i as the Legendre transform ofĤ i . But if we are only interested in the identity of the species that survive at community equilibrium, it is wholly encoded in the location of the extremum h * at β = ∞. Since this h * is located at the boundary of Ω, the shape of the repulsive potential of interaction with the hyperplanes ∆ µ = 0 is irrelevant in this limit, and can be replaced by the Heaviside theta-function θ(−∆ µ ). In this limit, the problem reduces to computing the extremum ofF over the unsustainable region, as stated in the main text.

RESOURCE SUPPLY MODELS
The model of MacArthur Different models of resource supply correspond to different expressions of the functionF . The renewable resource of MacArthur, with renewal rate r and maximum resource availability K, is described by the following resource depletion rule (originally derived in Ref. [14]; see also the Supplementary section A in Ref. [16]): Integrating H(T ), we findĤ(T ) = KT − KT 2 2r , and therefore We conclude that for the resource model of MacArthur: The constant-supply model of Ref. [16] is a simpler model that postulates that a fixed amount of resource R is evenly divided among all competitors: H(T ) = R/T . In this model, we haveĤ(T ) = R log T and T * = R/h * .
Consider the close-to-symmetric scenario, where the supply of all resources is similar. If the cost of all strategies is close to χ 0 , then the availability of resources at equilibrium will be close to χ 0 as well. Linearizing around this point, a general resource supply model H(T ) can be characterized with two parameters. First, let τ be the value of demand at which resource availability hits χ 0 : by definition, H(τ ) = χ 0 . In the vicinity of this point, let γ be the "elasticity" of supply, describing how quickly resource is depleted by a small increase in demand δτ τ : Let us compute the Legendre transform ofĤ in this model. We have: The demand that corresponds to a particular value of resource availability close to χ 0 : In the vicinity of h = χ 0 it is convenient to work with shifted variables: h ≡ χ 0 − g N (the N in the denominator reminds that the deviation is small). After a little algebra we find: where the index i reminds us that parameters τ i and γ i could be different for different resources. Omitting the irrelevant global constant, we find the expression forF in this general cost model: Of course, the two particular models we considered above reduce to this same form in the vicinity of h i ≈ χ 0 . Specifically, for MacArthur's model of renewable resource: Similarly, for the constant supply model: Our calculation below will be for the general close-to-symmetric case where the supply of resources is similar: Here i δτ i = i δγ i = 0 by definition, and as N becomes large, δτ i and δγ i remain of order 1. Note that the functioñ F can be rescaled by a constant positive factor, leaving the maximization problem unchanged (we seek the location of the maximum, not its magnitude). Without restricting generality, therefore, we can setτ = 1.

THE COST MODEL: AN ILLUSTRATION
The main text made the argument that since competition is restricted to only a subset of species, the self-selected pool of low-cost outliers with similar costs χ µ , the details of the cost model matter only inasmuch as they determine the properties of this subset. To illustrate this point, consider a scenario at N = 2 where the cost of strategy {x, 1−x} is drawn out of a normal distribution with mean χ 0 (1 + sin πx) and width 1 2 χ 0 sin πx. One such realization for 20 equally spaced values of 0 ≤ x ≤ 1 is shown in Fig. S2. In this illustration, mixed strategies tend to be expensive; as a result, all three low-cost outliers (in bold) are close to specialists, and only these species are competitive (this is exactly the scenario depicted in Fig. 1B). However, the exact details of the cost model (the precise shape of the solid red curve in Fig. S2) are otherwise irrelevant for the coexistence problem.

THE REPLICA-THEORETIC CALCULATION
This section demonstrates how the geometrical problem formulated above can be solved using methods of statistical physics, specifically an approach termed "replica theory". An attempt is made to present this computation in a detailed and self-contained way, i.e. not assuming familiarity with statistical physics of disordered systems. For a more comprehensive introduction to this powerful technique, we refer the reader to Ref. [19].

The basic idea
We seek to compute: S2. The cost model: an illustration. Competition is restricted to a self-selected pool of low-cost outliers, so we only need to model this pool, where all strategies have similar costs. This illustration shows an example of a cost model at N = 2 yielding the competition scenario depicted in Fig. 1B of the main text (the cost of strategy {x, 1 − x} (at N = 2) is plotted as a function of x). Since metabolic diversification is penalized, all three competitive species (highlighted) are close to being specialists, but otherwise the details of the red curve have no effect on the coexistence problem.
The pool of competitors is modeled as follows. For each species µ, we pick σ µ as a random binary vector, where each component σ i µ is 1 with probability p, and 0 otherwise. We then draw a random cost χ µ = i σ i µ + x µ , where x µ is a Gaussian random variable of variance 1.
The argument of the Heaviside θ-function in Eq. (S8) is the negative resource surplus −∆ µ . Under the cost model described above, we have Change variables h i ≡ 1 − gi N (g i runs from −∞ to N ), and introduce ∆ µ as an explicit auxiliary integration variable: In this expression, x µ and σ i µ are "frozen disorder": they are drawn randomly, but are then kept fixed, while other variables relax to their equilibrium values. Computing this integral for a particular realization of the disorder, even if it were possible to do so, would not be very informative. Instead, we are interested in the behavior of the "typical" realization of the system. This means that we are interested in the typical free energy F = log Z (angular brackets denote averaging over disorder). This quantity is hard to compute directly, because the average is outside of the logarithm. The opposite case, the logarithm of the average, would be very simple to compute; unfortunately, unlike the free energy F , the partition function Z is dominated not by typical realizations of the disorder, but by extreme ones. The logarithm of the average would capture the behavior of the system in highly improbable extreme cases, which is of no use to us.
The gist of the "replica trick" is summarized in the following formula: This trick makes it possible to formally derive an expression for log Z by computing only expressions of type Z n , which is the partition function of n copies ("replicas") of the system, and then formally sending n to zero. The replicas are identical (have the same disorder), but independent (each has its own set of degrees of freedom). There are, of course, mathematical subtleties related to taking this limit, and we refer the reader to Ref. [19]. This reference also provide some intuition for the physical basis of the argument and the interpretation of the auxiliary variables that appear along the way. Here, the analytical result we derive will be validated by an excellent agreement with numerical simulations.
Proceeding with our argument, we write the partition function of n copies of our system. Each degree of freedom is copied n times, labeled by the replica index a running from 1 to n: We stress that all n replicas have the same disorder: quantities x µ and σ i µ have no index a. Thanks to our strategic choice of cost model, averaging Z n over these disorder variables is separable: To compute the second term, recall that σ i µ is either 1 or 0 with probabilities p and 1 − p: Here we used 1 + p + p 2 2 + · · · = exp(p + p(1 − p) 2 /2 + . . . ). Putting everything together: In order to make progress, we need to eliminate terms that directly couple indices µ and i. To achieve this, we introduce a yet another set of variables, using the same trick of inserting a delta-function into our integral: These auxiliary variables will become our order parameters, as they capture the mean and variance of the deviation of resource availability from 1. The indices i and µ are now decoupled, and we can split the integral accordingly. Recall that we have index i labels N different resources, while µ labels P different strategies. Recall also thatF is a sum over N terms:F ≡ N i=1F i .
with A i and B given by: After averaging over disorder, the problem becomes fully symmetric in indices µ (all strategies are drawn from the same distribution, so there is no inherent difference in how they contribute). The same would be true for resources, except in the interest of generality, we allowed the supply functions H i to be different for different resources.
Decoupling replicas: the replica-symmetric ansatz The idea now is to treat the integrals over m and q as saddle-point integrals. We will shortly introduce a rescaling of variables that will make β appear in the exponent to serve as the large parameter, making the saddle-point approximation appropriate. In this approximation, the integral is replaced by the value of the integrand at one location, the saddle point (up to a multiplicative prefactor, which, as explained below, is irrelevant for our purposes). When looking for the saddle point, we will make the assumption that it is symmetric under a permutation of replicas. This is the so-called "replica-symmetric ansatz". The validity of this assumption will be justified a posteriori by the fact that the saddle-point we will find is "well-behaved" and the analytical results match the numerical simulations.
At a fully replica-symmetric saddle point, all components of m a must coincide. As for the matrix q ab , all its diagonal elements must be equal, and all of the off-diagonal ones must be equal as well. The same holds for the conjugate variablesm a ,q ab . We therefore look for a saddle point of the following form: and similarly for q ab and m a . With these assumptions: log Z n = extr in q DqD + i n(n − 1) 2 q OqO + inm * m * + i log A i + P log B .
(Why did we take the logarithm? Recall from Eq. S9 that our ultimate goal is to compute Z n − 1 in the n → 0 limit. Conveniently, in this limit, subtracting 1 is the same as taking the logarithm.) We stress that all expressions need only be computed to the leading exponential order. In particular, multiplicative constants are irrelevant as they only amount to an additive constant under the logarithm. Below, such constants will be omitted, and the "equal" signs will mean "up to a constant multiplicative factor".
The limit n → 0 Note that in the expression we just found, n enters in a way that allows taking the formal limit n → 0, which is our final goal. It turns out that this limit also makes the expressions for log A i and log B somewhat easier to compute. Therefore, we write: For reasons that will become clear later, we note that this can also be written as Computing log Ai Recall the expression we denoted A i : We first write: Now use Feynman's trick of removing the square by introducing an extra Gaussian variable: (the curly D denotes the standard Gaussian measure with variance 1). This lets us write: At the price of introducing an extra Gaussian variable, all replicas are now fully decoupled. Plugging this into our expression for A i : Conveniently, for small n: log Dz x n = log Dz (1 + n log x + . . . ) = log 1 + n Dz log x + . . . = n Dz log x + . . . Therefore: Introduce rescaled variables as follows: In the new variables: Substitute the general form ofF i from (S6), for a close-to-symmetric resource supply (S7): In the limit β → ∞: Here y ≡ g N . The quadratic form −P y 2 + Qy is maximized at y = Q 2P , reaching the maximal value of Q 2 4P . In our case, therefore, for a given z the quadratic form reaches its maximum at Let us shift the variablem by definingm ≡ −τ − δm √ N . For now, we can treat this as a simple change of variables. We find: Recall thatm is one of the variables over which the extremum is computed in (S10). The shift ofm that we just did constitutes an assumption, namely that the extremum is located close tom ≈τ (i.e. that the difference is at most of order 1/ √ N ). We will check the consistency of this assumption below. Note that if we are correct to assume this, then y * i (z) is small, justifying the expansion of resource depletion functions to first order in y . The integral over z is now simple to compute: Finally, performing the sum over i and recalling that i δτ i = 0, we find a very simple final expression: where δτ 2 ≡ 1 N i (δτ i ) 2 is the variance of resource supply across i.
Computing log B Recall the definition of B: We then remove the square by introducing an extra Gaussian variable, making all replicas fully decoupled: Note that the integral over∆ inside the square brackets is a simple Gaussian integral, and we can write: 2π e −y 2 /2 . This is essentially the error function, up to a couple constants that would be a nuisance to carry around: We can then write: As above for A, we now introduce a rescaled variable x, and some convenient notations: In the new variables: The logarithm of E(x) ≡ ∞ x dy √ 2π e −y 2 /2 in the large-argument limit is very simple. Indeed: Therefore (omitting additive constants as always): Plugging this into our expression for log B, we find where I(λ) can be expressed in terms of the error function: Putting everything together Combining the results above, plugging them into (S10), and recalling that in the large-N limit, P also goes to infinity with P N ≡ α, we find: Recall that ψ ≡ p(1 − p)q + 2 , so that the extremum is taken over six variables: δm, a, b, q, λ and x. Consider the extremum condition for δm: As N → ∞, we therefore have δm → 0, demonstrating that the approximationm ≈τ was indeed self-consistent. Setting δm = 0, we find: The powers of N and β in our rescaled variables were chosen to ensure that this expression no longer depends on N and is proportional to β. Conveniently, the extremum conditions for variables a and b can also be solved, and these variables eliminated: Our final expression for the partition function: The extremum is to be computed over q, x and λ. Here α, p and (hidden in ψ ≡ p(1 − p)q + 2 ) are parameters characterizing the pool of competitors (number of strategies, typical functional sparsity, and intrinsic scatter cost, respectively). Parametersτ , δτ 2 andγ characterize resource supply (respectively: average capacity, variability across resources, and average "elasticity"). For the simplest single-parameter resource model used in the text, the influx of resource i is fixed at R i , whose average, without loss of generality, can be set to 1. In this case we haveτ = 1,γ = 1 and the remaining parameter is the variance of resource supply δτ 2 , denoted δR 2 in the main text.

THE SADDLE-POINT EQUATIONS
Simplifying the equations to solve them numerically Hidden in I(λ) is the error function erfc, which means that the extremum of (S12) cannot be found analytically. However, the equations can be simplified to a form where they can either be solved numerically, or investigated analytically in certain limits.
The extremum conditions: Using the first and the third equations, we write: Rearranging, we find a way to express ψ (and thus q) in terms of λ only: Plugging this into the first equation, and recalling the definition of I(λ), we find a very simple expression for x: .) This result makes it possible to eliminate x from the equations. Plugging all this into (S14), we find an equation that involves λ only, and can easily be solved numerically. Once λ is known, equations (S15) and (S16) determine q and x.
Investigating the limit → 0 To study the equations analytically in the limit → 0, we plug the expression for x into (S13), and reorganize the terms in equation (S14), putting it into the form cited in the main text: The easiest way to derive the expression for the critical line given in the main text is to observe the following. At = 0, we see that simultaneously setting ψ = 0 and 1 − αE(λ) = 0 yields a solution. One can check that in the vicinity of the transition, both go to zero linearly, so that their ratio remains well-defined. Omitting the negligible first terms in the right-hand sides of both equations, we rewrite them as follows: At the critical line, the system is degenerate, i.e. the two equations are proportional. We immediately read off the condition that must be satisfied: This is the expression quoted in the main text (with λ being fixed by the condition 1 − α crit E(λ) = 0).

COMPUTING THE NUMBER OF SURVIVORS AT EQUILIBRIUM
To find the number of species that survive at equilibrium in our model, we set out to compute the distribution of the observable ∆. Recall that ∆ is the "resource surplus"; for a given set of competitors, there is a discrete set of values of ∆ µ , the resource surplus experienced by each species. However, after we average our partition function over disorder, ∆ becomes a random variable, drawn out of a certain distribution, whose shape is encoded in the partition function. This distribution is what we now set out to compute.
The ∆-dependent part of the partition function is fully contained in the expression for B. Recall the intermediate expression (S11) derived earlier, when computing B: Here we introduced β ≡ β p(1−p)N x to make the notations slightly less heavy. Let us put this expression in the following form, retaining one copy of the integral over ∆, while evaluating the remaining n − 1 copies as before: Once again, E(x) is a short-hand for 1 2 erfc(x/ √ 2). Sending n → 0 and recalling the notation λ ≡ pm ψ : From this we infer the distribution of ∆ for ∆ ≤ 0 (positive ∆ are forbidden): This is a complicated-looking expression, but its β → ∞ limit can in fact be computed very easily, using the following trick. By definition of the function E(x), this expression can be rewritten as follows: The key observation is that here, the "conditional distribution" p(∆|w) is (by inspection) a properly normalized distribution, for all w: ∀w : ∞ 0 d∆ p(∆|w) = 1 (S17) We now compute the limit lim β →∞ p(∆|w). For w > −λ, this is a Gaussian distribution of width 1 √ β → 0 centered at ∆ = −ψ(λ + w) < 0. For w < −λ, the probability density is highest at ∆ = 0 and goes to zero everywhere else (recall that positive ∆ are forbidden). Therefore, the normalization condition (S17) immediately tells us that: As a result: This is the expression quoted in the main text. The weight of the delta-shaped peak at ∆ = 0 corresponds to the species whose resource balance is met; these are the species who survive competition.

VALIDATION OF ANALYTICAL RESULTS AGAINST SIMULATIONS
Numerical simulations were performed in MatLab; a script generating all figures is available as Supplementary File 1. Random instances were generated by randomly drawing strategy vectors and costs as described in the main text. All simulations used N = 50 for the number of resources. The equilibrium of a community was determined by a direct N -dimensional numerical optimization ofF = i R i /T i subject to the linear constraints defining the "unsustainable region" Ω. The uneven resource supply was always implemented as a bimodal distribution shown in Fig. 3B (upper panel); the amplitude of the step was adjusted to match the required magnitude of δR 2 . Figure 2B 500 simulations were performed at α = 10, = 10 −3 and δR 2 = 1 (reusing the dataset computed for Fig. 3B). These parameters are comfortably in the S-phase, so the expected number of survivors is N . After numerical equilibration, all but the top N values of the resource surplus ∆ were recorded (the top N are within numerical error of 0 and correspond to survivors). Panel 2B shows a joint histogram of these values recorded over 500 instances, appropriately normalized to be comparable with the theoretical distribution. For visualization purposes, the delta-shaped peak is shown as a rectangle, its height is meaningless. The theoretical prediction for the number of survivors is verified in panel 2C.

Figure 2C
The theoretical curves are overlaid with simulations obtained as follows. For 10 values of α equispaced between 2 and 5, and for three values of indicated on the plot (0.1, 0.03 and 10 −4 ), communities were equilibrated using MatLab solver fmincon with numerical precision parameter 10 −10 . Species whose resource surplus was within 10 −8 of zeros were declared as survivors (for comparison, the resource surplus of the first extinct species was typically ≈ 10 −6 ). The panel shows the mean number of survivors over 500 instances for each (α, ) pair.   For = 10 −3 and δR 2 = 1, we performed 500 simulations at α = 2 for V phase, and at α = 10 for S-phase. The latter dataset was also used for generating the resource surplus histogram shown in 2B. . Define a species "success rank" by ordering all survivors by decreasing abundance at equilibrium, followed by all the species that went extinct, in order of increasing resource insufficiency |∆µ|. Now, consider two quantities that could potentially be used to predict this success. To obtain the first, we measure the growth rate of each species when it is placed into the externally supplied conditions, with no competitors present. Alternatively, the "null model" performance predictor is simply the species' cost per pathway: clearly, the high-cost species are less likely to survive. We call the latter quantity the "null model", because it includes no information about the environment to which the community is subjected. We now ask if either of the two quantities are in fact predictive of the true success of a species. We measure predictive power by the Spearman (rank-order) correlation coefficient between the predictor being tested and the true success rank. In the V-phase, the species' performance measured in external conditions is indeed predictive of the success rank (panel A.) A histogram of the correlation coefficients over 500 instances (panel B) confirms that this predictor significantly outperforms the null model. In contrast, in the S-phase, the environment-specific performance measured in the externally imposed conditions becomes irrelevant (panels C and D). Although the correlation observed in panel C is non-zero, this is due to the generic fact that low-cost organisms are generally more likely to survive that high-cost ones. This is demonstrated in panel D: unlike the more intuitive V-phase, in the S-phase the "environment-aware" predictor performs worse than the null model.