Consider N = 10 equiabundant substrates, and one random realization of organism costs with scatter ϵ = 10 - 3 . (MATLAB scripts (MATLAB, Inc.) performing simulations and reproducing all figures are available as Supplementary file 2). The numerical simulation of competition between all 1023 possible species, initialized at equal abundance, results in the equilibrium state depicted in Figure 2A, Supplementary file 2. In this example, it consists of nine species. It is natural to ask: for a given initial set of competitors, what determines the species that survive?

Figure 2 Download asset Open asset The individual performance rank of a species (its cost per pathway) is predictive of its survival and abundance in a community. (A) Community equilibrium for N = 10 substrates with abundance R i / χ 0 = 100 and one particular random realization of organism costs with scatter ϵ = 10 - 3 . Species are ordered by abundance and labeled by the pathways they carry. Also indicated is the individual performance rank; all surviving species were within the top 30 (out of 1023). (B) The median individual performance rank of survivors, weighted (dashed) or not weighted (solid) by abundance. Curves show mean over 100 random communities for each value of cost scatter ϵ ; the standard deviation across 100 instances is stable at approximately 40% of the mean for both curves, independently of ϵ (not shown to reduce clutter). https://doi.org/10.7554/eLife.15747.004

In the present model, the only intrinsic performance characteristic of a species is its cost per pathway. Consider an assay whereby a single individual of species σ → is placed in an environment with no other organisms present, and, for simplicity, an equal supply of all substrates R i = R . The initial population growth rate in this chemostat is given by:

d ⁢ n σ → d ⁢ t | t = 0 = 1 τ 0 ⁢ χ σ → ⁢ [ ∑ i R i ⁢ σ i - χ σ → ] = 1 τ 0 ⁢ [ R ⁢ | σ → | χ σ → - 1 ] ,

and abundance eventually equilibrates at n σ → = R ⁢ | σ → | / χ σ → . Both these quantities characterize performance of species σ → (the term 'fitness' is avoided as it is a micro-evolutionary concept that, strictly speaking, should be defined only within individuals of one species), and both are determined by the inverse cost per pathway. Define the 'individual' performance measure of species σ → as

(4) f σ → ≡ χ 0 ⁢ | σ → | χ σ → - 1.

This definition is convenient as it makes f σ → a dimensionless quantity of order ϵ . Under the cost model (2), the performance ranking of species is random, set by the random realization of the costs ξ : f σ → = 1 1 + ϵ ξ σ → − 1 ≈ − ϵ ξ σ → . This model was chosen so that no group of species has an obvious advantage. A different cost function would effectively reduce the pool of competitors, excluding certain (prohibitively expensive) species from the start.

Predictably, this performance ranking is correlated with the success of a species in a community, but not very well (Figure 2). The equilibrium depicted in panel A predominantly consists of top-ranked species, and the median performance rank of surviving species is consistently low across a range of values of the cost scatter ϵ (panel B). This median rank becomes even lower if the median is weighted by a species’ abundance at equilibrium, indicating that top-ranked species tend to be present at higher abundance (Davis et al., 1998; Birch, 1953). Still, at the equilibrium shown in Figure 2A, the species ranked fourth in intrinsic performance went extinct, but six others ranked as low as #29 remained present.

These observations reflect the well-known fact that the success of a species is context-dependent and observing a species in isolation does not measure its performance in the relevant environment (McGill et al., 2006; McIntire and Fajardo, 2014). In the model described here, the context experienced by all species is fully encoded in the vector of 'harvests', that is, the benefit an organism receives from carrying pathway i :

(5) H i ≡ R i / T i .

A growing demand for substrate i (increasing T i ) depletes its availability, in the sense that H i is reduced. Consider the three-substrate world depicted in Figure 1, and assume that A ⁢ B ¯ is the highest-performing species with a very low cost. As A ⁢ B ¯ multiplies, it depletes resources A and B . As a result, the final equilibrium is highly likely to include the specialist organism C ¯ , even if its cost is relatively high, and under other circumstances (if A ⁢ B ¯ were less fit) it would have yielded to A ⁢ C ¯ or B ⁢ C ¯ .

Conveniently, in the model described here, these complex effects studied by niche construction theory can be summarized in a single community-level objective function. The dynamics (3) possess a Lyapunov function, i.e. a quantity that is increasing on any trajectory of the system (compare to MacArthur, 1969):

(6) F = 1 R tot ⁢ ( ∑ i R i ⁢ ln ⁡ T i R i / χ 0 - ∑ σ → χ σ → ⁢ n σ → + R tot ) .

Here, R tot is a constant introduced for later convenience. Specifically, set R tot = ∑ i R i ; this choice ensures that close to community equilibrium, F is also of order ϵ (see Supplementary file 1, section C). This function has the property that R tot ⁢ ∂ ⁡ F ∂ ⁡ n σ → = Δ σ → (the resource surplus), and therefore

d ⁢ F d ⁢ t = ∑ σ → ∂ ⁡ F ∂ ⁡ n σ → ⁢ d ⁢ n σ → d ⁢ t = ∑ σ → n σ → ⁢ ( Δ σ → ) 2 R tot ⁢ χ 0 ⁢ τ 0 > 0

Thus, F is indeed monotonically increasing as the system is converging to equilibrium. To illustrate this, Figure 3A shows 10 trajectories of ecological dynamics for the same system as in Figure 2A, starting from random initial conditions (with all species present; see Supplementary file 1, section H). Far from equilibrium, while most high-cost species are being eliminated by competitors, the mean intrinsic performance of surviving organisms and F increase together (Figure 3A, inset), confirming that intrinsic performance is a useful predictor. However, as equilibrium is approached, community-induced changes in substrate availability H i reduce the relevance of the original performance ranking, which was measured in the 'wrong' environment; previously successful species can be driven to extinction (Figure 3B). The set { H i } at equilibrium characterizes the environment the surviving species had 'carved' for themselves. The performance rank ordering will be all the more sensitive to the environment { H i } , the smaller the scatter of intrinsic organism costs ϵ . Therefore, the role of this parameter is to tune the relative magnitude of intrinsic and environment-dependent factors in determining a species’ fate. So far, ϵ was fixed at 10 - 3 ≈ 2 - N , and Figure 2B shows that for small cost scatter ϵ , the structure of the final equilibria does not significantly depend on this parameter (see Supplementary file 1, section D). The large- ϵ regime will be discussed later.

Figure 3 Download asset Open asset Community dynamics maximize a global objective function. (A) 10 trajectories of an example system, starting from random initial conditions and converging to the equilibrium depicted in Figure 2A. Direction of dynamics indicated by arrows. Far from equilibrium, mean intrinsic performance of members (weighted by abundance) and the community-level function F increase together (inset; data aspect ratio as in the main panel). Close to equilibrium, intrinsic performance loses relevance. (B) Time traces of species’ abundance for one community trajectory (thick red line in A). Arrowheads in panels A and B indicate matching time points. Species that eventually go extinct shown in red; many enjoy transient success. (C) The complex dynamics of panel B is driven by the simple objective to efficiently deplete all substrates simultaneously, encoded in F . Shown is mean availability of the 10 substrates, for each trajectory of panel A. https://doi.org/10.7554/eLife.15747.005

Each of the trajectories in Figure 3A converges to the same equilibrium (depicted in Figure 2A). This is because F is convex and bounded from above (see Supplementary file 1, section A). Therefore, for every set of species Ω , any community restricted to these species will always reach the same (stable) equilibrium, corresponding to the unique maximum of F within the subspace where only species of Ω are allowed non-zero abundance. This maximum will often be at the border of this subspace, corresponding to the extinction of some species.

Under the dynamics (3), no new species can 'appear' if their original abundance was zero. Imagine, however, that on each island, a rare mutation (or migration) occasionally introduces a random new species; if it can invade, the community transitions to a new equilibrium and awaits a new mutation. This process of adaptive dynamics defines the evolution of each island, and can be seen as a mesoscopic population genetics model for a multi-species community evolving through horizontal gene transfer (loss/acquisition of whole pathways). For each island, F is monotonically increasing throughout its evolution. Indeed, F is continuous and non-singular in all n σ → , so introducing an invader at a vanishingly small abundance will leave F unchanged, and the subsequent convergence to a new equilibrium is a valid trajectory of ecological dynamics on which F increases.

What is the intuitive meaning of the function F that is being optimized by the community? It is easy to show that ∑ σ → χ σ → n σ → = ∑ i R i at any equilibrium (total demand matches the total supply; see Supplementary file 1, section C). Therefore, for a community at equilibrium, F characterizes its ability to deplete substrates:

(7) F = - ∑ i R i ⁢ ln ⁡ H i R tot + const

If all substrates are equiabundant for simplicity, maximizing F is equivalent to minimizing ∑ i ln ⁡ H i . The optimization principle that appears in this model is therefore a generalization of Tilman’s R * rule (Tilman, 1982). In the classic form, this rule states that for a single limiting resource, the unique winning species is the one capable of depleting this resource to the lowest concentration. However, if the limiting resource can be harvested from multiple substrates, as considered here, multiple species may coexist; the winning community is the one that is most efficient at depleting all substrates simultaneously, weighted as described in Equation (7). This is illustrated in Figure 3C. While the time trajectories of individual species may be highly complex (Figure 3B), the net effect of these dynamics is to deplete substrate availability down to the lowest concentrations capable of sustaining a population (see also Figure S1 in Supplementary file 1, section H).