Nutrient levels and trade-offs control diversity in a serial dilution ecosystem

TL;DR

This study shows that in serial dilution microbial ecosystems, metabolic trade-offs and nutrient supply levels significantly influence species diversity, revealing complex relationships that challenge universal nutrient-diversity assumptions.

Contribution

It provides a theoretical framework demonstrating how metabolic trade-offs enable high microbial diversity in serial dilution cultures with varying nutrient supplies.

Findings

High diversity is possible with metabolic trade-offs in serial dilution cultures.

Diversity depends on nutrient amount due to an 'early-bird' effect.

Relationships between nutrient supply and diversity are complex and context-dependent.

Abstract

Microbial communities feature an immense diversity of species and this diversity is linked with outcomes ranging from ecosystem stability to medical prognoses. Yet the mechanisms underlying microbial diversity are under debate. While simple resource-competition models don't allow for coexistence of a large number of species, it was recently shown that metabolic trade-offs can allow unlimited diversity. Does this diversity persist with more realistic, intermittent nutrient supply? Here, we demonstrate theoretically that in serial dilution culture, metabolic trade-offs allow for high diversity. When a small amount of nutrient is supplied to each batch, the serial dilution dynamics mimic a chemostat-like steady state. If more nutrient is supplied, diversity depends on the amount of nutrient supplied due to an "early-bird" effect. The interplay of this effect with different environmental factors and diversity-supporting mechanisms leads to a variety of relationships between nutrient supply and diversity, suggesting that real ecosystems may not obey a universal nutrient-diversity relationship.

Tables

Table 1. Table 1: Annotation glossary

Symbol	Description
$t$	Time measured from the beginning of a batch
$p$	Number of nutrients
$m$	Number of species introduced at time $t=0$
$m_e$	Effective number of species at steady state
$i$	$\left(1\dots p\right)$ Latin index enumerating nutrients
$c_i\left(t\right)$	Time dependent concentration of nutrient $i$
$c_0$	${\sum }_{i=1}^p{c}_i\left(0\right)$; total nutrient concentration at time $t=0$
$K_i$	Monod half-velocity constant for nutrient $i$
$I_i$	${\int }_0^{\infty }\frac{c_i}{K_i+c_i}𝑑t$; nutrient Monod function time integral
$Y_i$	Biomass yield on nutrient $i$
$\Delta$	The fraction of nutrient remaining at the end of the batch
$s_i$	Average rate that nutrient $i$ is continuously supplied at the chemostat limit
$\tilde{\delta }$	Average rate all nutrients are continuously supplied at the chemostat limit
$\delta$	Continuous chemostat dilution rate at the chemostat limit
$\sigma ,{\sigma }^{\prime },\dots$	$\left(1\dots m\right)$ Greek indices enumerating species
${\rho }_{\sigma }\left(t\right)$	Species $\sigma$ biomass density at time $t$ since a start of the batch
${\overrightarrow{\alpha }}_{\sigma }$	$({\alpha }_{\sigma ,1},\dots ,{\alpha }_{\sigma ,p})$; enzyme allocation strategy for species $\sigma$
$\epsilon$	Standard deviation in enzyme budget
$E$	$E={\sum }_i{\alpha }_{\sigma ,i}=1$ for $\epsilon =0$; enzyme budget
${\Gamma }_{i,i^{\prime }}^{\sigma }$	Byproduct matrix converting nutrient $i^{\prime }$ to nutrient $i$
$j_{\sigma ,i}$	Nutrient $i$ consumption rate by species $\sigma$

Equations

$$j_{\sigma ,i}$$

$$\frac{d{c}_i}{dt}$$

$$\frac{d{\rho }_{\sigma }}{dt}$$

$$\frac{c_0}{{\rho }_0}=\sum _i\frac{{\alpha }_{\sigma ,i}{c}_i(0)}{{\sum }_{{\sigma }^{\prime }}{\alpha }_{{\sigma }^{\prime },i}{\rho }_{{\sigma }^{\prime }}(0)}+O{\left(\frac{c_0}{K}\right)}^2.$$

$$\tilde{\delta }=\frac{c_0}{{\rho }_0{t}_f},s_i=\frac{c_i(0)}{t_f},$$

$$\tilde{\delta }=\sum _i\frac{{\alpha }_{\sigma ,i}{s}_i}{{\sum }_{{\sigma }^{\prime }}{\alpha }_{{\sigma }^{\prime },i}{\rho }_{{\sigma }^{\prime }}(0)}.$$

$$I_i={\int }_0^{\infty }\frac{c_i}{K_i+c_i}dt=\text{mbox}const.$$

$$\frac{d{c}_i}{dt}=-\sum _{\sigma }{\rho }_{\sigma }\left(j_{\sigma ,i}-\sum _{i^{\prime }}{\Gamma }_{i,i^{\prime }}^{\sigma }{j}_{\sigma ,i^{\prime }}\right).$$

$$\frac{d{\rho }_{\sigma }}{dt}$$

$$\frac{d\widetilde{c_i}}{dt}$$

$$E$$

$$\frac{d{\rho }_{\sigma }}{dt}$$

$$\frac{d{c}_i}{dt}$$

$$1$$

$$\frac{d{\rho }_{\sigma }}{dt}$$

$$\frac{d{c}_i^{\prime }}{dt}$$

$$1$$

$${\rho }_{\sigma }(t)$$

$$c_i(t)$$

$${\rho }_{\sigma }^{(k>0)}(0)$$

$$c_i^{(k>1)}(0)$$

$$\frac{c_i}{c_i+K}$$

$$\dot{c_i}=-\frac{c_i}{K+c_i}\sum _{\sigma }{\alpha }_{\sigma ,i}{\rho }_{\sigma },$$

$${\dot{c}}_i^{(1)}$$

$$⟹c_i^{(1)}$$

$${\dot{\rho }}_{\sigma }={\rho }_{\sigma }\sum _i{\alpha }_{\sigma ,i}\frac{c_i}{K+c_i}.$$

$${\dot{\rho }}_{\sigma }^{(1)}$$

$$⟹{\rho }_{\sigma }^{(1)}$$

$${\rho }_{\sigma }^{(1)}(t\gg {\gamma }_i^{-1})={\rho }_{\sigma }(0)\sum _i\frac{{\alpha }_{\sigma ,i}{c}_i^{(1)}(0)}{K{\gamma }_i}.$$

$${\rho }_{\sigma }(t\gg {\gamma }_i^{-1})={\rho }_{\sigma }(0)+\varphi {\rho }_{\sigma }(0)\sum _i\frac{{\alpha }_{\sigma ,i}{c}_i^{(1)}(0)}{K{\gamma }_i}.$$

$${\rho }_{\sigma }(t\gg {\gamma }_i^{-1})\approx {\rho }_{\sigma }(0)+{\rho }_{\sigma }(0)\sum _i\frac{{\alpha }_{\sigma ,i}{c}_i(0)}{{\sum }_{{\sigma }^{\prime }}{\alpha }_{{\sigma }^{\prime },i}{\rho }_{{\sigma }^{\prime }}(0)}.$$

$${\rho }_{\sigma }(d,t_f)\approx \left(1+\sum _i\frac{{\alpha }_{\sigma ,i}{c}_i(0)}{{\sum }_{{\sigma }^{\prime }}{\alpha }_{{\sigma }^{\prime },i}{\rho }_{{\sigma }^{\prime }}(d,0)}\right){\rho }_{\sigma }(d,0).$$

$${\rho }_{\sigma }(d+1,0)$$

$$1+c_0/{\rho }_0=1+\sum _i\frac{{\alpha }_{\sigma ,i}{c}_i(0)}{{\sum }_{{\sigma }^{\prime }}{\alpha }_{{\sigma }^{\prime },i}{\rho }_{{\sigma }^{\prime }}(d,0)}.$$

$$\frac{c_0}{{\rho }_0}=\sum _i\frac{{\alpha }_{\sigma ,i}{c}_i(0)}{{\sum }_{{\sigma }^{\prime }}{\alpha }_{{\sigma }^{\prime },i}{\rho }_{{\sigma }^{\prime }}(0)}.$$

$$\tilde{\delta }=\frac{c_0}{{\rho }_0{t}_f},s_i=\frac{c_i(0)}{t_f},$$

Full text

Introduction

Microbial communities feature an immense diversity of organisms, with the typical human gut microbiota and a liter of seawater both containing hundreds of distinct microbial types [1, 2, 3]. These observations appear to clash with a prediction of some resource-competition models, known as the competitive-exclusion principle – namely, that steady-state coexistence is possible for only as many species as resources [4, 5]. This conundrum is familiarly known as the “paradox of the plankton” [6]. Solving this paradox may provide one key to predicting and controlling outcomes ranging from ecosystem stability to successful cancer treatments in humans [7, 8, 9, 10]. Chesson [11] classified mechanisms that purport to solve this paradox into two broad categories: stabilizing and equalizing. Stabilizing mechanisms prevent extinction by allowing species to recover from low populations, whereas equalizing mechanisms slow extinction by minimizing fitness differences between species.

Many possible solutions of the paradox that rely on stabilizing mechanisms have been offered: (i) interactions between microbes, such as cross-feeding or antibiotic production and degradation [12, 13], (ii) spatial heterogeneity [14, 15], (iii) persistent non-steady-state dynamics [6], and (iv) predation [16]. Equalizing mechanisms have been studied through neutral theory, in which all species are assumed to have equal fitness [17], and recent work has proposed resource-competition models that self-organize to a neutral state [18]. Many proposed solutions for the paradox assume a chemostat framework wherein nutrients are continuously supplied and there is a continuous removal of biomass and unused nutrients [19]. However, in nature nutrients are rarely supplied in a constant and continuous fashion. In particular, seasonal variation is ubiquitous in ecology, influencing systems ranging from oceanic phytoplankton communities [20] to the microbiota of some human populations [21]. How does a variable nutrient supply influence diversity?

Existing literature on seasonality focuses on stabilizing mechanisms and generally finds that seasonality either promotes or has little effect on coexistence [22]. But do these conclusions extend to equalizing mechanisms? To address this question, we consider a known resource-competition model that permits high diversity at steady state due to the equalizing effects of metabolic trade-offs, which assume that microbes have a limited enzyme production capacity they must apportion. Here, we investigate the equalizing effect of metabolic trade-offs in the context of serial dilution, to reflect a more realistic variable nutrient supply.

Serial dilution, in which cultures of bacteria are periodically diluted and supplied with fresh nutrients, is well-established as an experimental approach. For example, the bacterial populations in the Lenski long-term evolution experiment [23], experiments on community assembly [24], and antibiotic cross-protection [25] were all performed in serial dilution. While previous models of serial dilution have characterized competition between small numbers of species with trade-offs in their growth characteristics [26, 27], the theoretical understanding of diversity in serial dilution is much less developed than for chemostat-based steady-state growth.

Here, we show that under serial dilution, metabolic trade-offs can still support high diversity communities, but that this coexistence is now sensitive to environmental conditions. Interestingly, seasonality reduces diversity in our model, contrasting what has been observed for stabilizing mechanisms. In particular, we find a surprising dependence between community diversity and the amount of nutrient provided to the community. These changes in diversity are driven by an “early-bird” effect in which species that efficiently consume nutrients that are initially more abundant gain a population advantage early in the batch. To our knowledge, this is the first time this effect has been identified as a major influencer of diversity in seasonal ecosystems.

This dependence between community diversity and the supplied nutrient concentration allowed us to explore an unresolved question in ecology [28, 29, 30]: what is the relationship between the amount of nutrient supplied and the resulting diversity of the community? Experimental studies of this question have mainly been performed in macroecological contexts [31, 32, 33], though recently there has been increased focus on microbial systems [34, 35]. In microbial experiments, some evidence has supported the “hump-shaped” unimodal trend predicted by many theories [36]. However, a meta-analysis by Smith [37] found no consistent trend across microbial experiments. What we observe here is concordant with Smith’s result: even in our highly simplified model, there is no general relationship between nutrient supply and diversity. Among the factors we find that influence this relationship are cross-feeding, relative enzyme budgets, differences in enzyme affinities, and differences in nutrient yields. That so much variation appears in a simple model suggests that real ecosystems are not likely to display a single universal relationship between nutrient supply and diversity.

Results

We employ the serial dilution model depicted in Fig. 1A. At the beginning of each batch ($t=0$), we introduce the inoculum, defined as a collection of species $\{\sigma\}$ with initial biomass densities $\rho_{\sigma}(0)$ in the batch such that the total initial biomass density is $\rho_{0}=\sum_{\sigma}\rho_{\sigma}(0)$. Together with the inoculum, we supply a nutrient bolus, defined as a mixture of $p$ nutrients each with concentration in the batch $c_{i}(0)$, $i=1,\ldots,p$ such that the total nutrient concentration is $c_{0}=\sum_{i=1}^{p}c_{i}(0)$ (we also consider the case of cycles of single nutrient boluses that approach a mixture distribution, cf. Appendix  F Fig. 11). It is assumed that all nutrients are substitutable (i.e. all nutrients are members of a single limiting class of nutrients, e.g. nitrogen sources). For simplicity, we assume ideal nutrient to biomass conversion, so that for a species to grow one unit of biomass density, it consumes one unit of nutrient concentration (we consider the case of nutrient-specific yields $Y_{i}$ in a later section). During each batch, the species biomass densities $\rho_{\sigma}(t)$ increase with time, starting at $t=0$, and growth continues until the nutrients are fully depleted, $\sum_{i=1}^{p}c_{i}(\infty)\approx 0$ (we consider the case of incomplete depletion in Appendix F Fig. 12). Thus, at the end of a batch, the total biomass density of cells is $\sum_{\sigma}\rho_{\sigma}(\infty)=\rho_{0}+c_{0}$. The next batch is then inoculated with a biomass density $\rho_{0}$ with a composition that reflects the relative abundance of each species in the total biomass at the end of the previous batch. This process is repeated until “steady state” is reached, i.e. when the biomass composition at the beginning of each batch stops changing.

In the model, a species $\sigma$ is defined by its unique enzyme strategy $\vec{\alpha}_{\sigma}=(\alpha_{\sigma,i},\ldots,\alpha_{\sigma,p})$ which determines its ability to consume different nutrients. We assume that each species can consume multiple nutrients simultaneously, in line with the behavior of microbes at low nutrient concentrations [38], though this assumption may not hold for all microbial species. Specifically, we assume that species $\sigma$ consumes nutrient $i$ at a rate $j_{\sigma,i}$ (per unit biomass) that depends on nutrient availability $c_{i}$ and on its enzyme-allocation strategy $\alpha_{\sigma,i}$ according to

[TABLE]

For simplicity, we take all Monod constants to be equal, $K_{i}=K$ (a more general form of the nutrient model is considered in a later section). During each batch, the dynamics of nutrient concentrations and biomass densities then follow from the rates $j_{\sigma,i}$ at which the species consume nutrients:

[TABLE]

Since the level of one enzyme inevitably comes at the expense of another, we model this trade-off via an approximately fixed total enzyme budget $E$. Formally, we take $\sum_{i}\alpha_{\sigma,i}=E+\varepsilon\,\xi_{\sigma}$, where $\xi_{\sigma}$ is a zero-mean and unit-variance Gaussian variable. Without loss of generality we take $E=1$; initially we set $\varepsilon=0$, which we call exact trade-offs. This allows us to visualize the strategies $\vec{\alpha}_{\sigma}$ as points on a simplex, depicted as colored circles embedded in: (i) the interval $[0,1]$ for two nutrients (Fig. 1B), or (ii) a triangle for three nutrients (Fig. 1C), etc. One can plot the nutrient bolus composition $c_{i}/c_{0}$ on the same simplex, as depicted by the black diamonds in Fig. 1B and C. In what follows, we focus on the case of two nutrients, though the main results extend to an arbitrarily large number of nutrients.

Connection between serial dilution and chemostat models

One can intuit that our serial dilution model at very low nutrient bolus size will mimic a chemostat. Adding a small nutrient bolus, letting it be consumed, then removing the additional biomass, and repeating is tantamount to operating a chemostat with a fixed nutrient supply and dilution rate. Indeed, the limit $c_{0}\ll K$ yields the same steady state as a chemostat. Thus, our results for serial dilution include and generalize those obtained for a closely related chemostat model [18].

For completeness, we now briefly describe the chemostat results from Ref. [18]. In the presence of metabolic trade-offs, the chemostat can support a higher species diversity than prescribed by the competitive exclusion principle as we demonstrate theoretically in Appendix  B. Specifically, if the nutrient supply lies within the convex hull of the strategies on the simplex (visualized by stretching a rubber band around the outermost strategies, see Fig. 1B-C), an arbitrarily large number of species can coexist at steady state. In the chemostat, such species coexistence is attained when the system organizes such that all nutrient levels are driven towards equality by consumption. Dynamically, if one nutrient level is high, the species that consume it increase in population, leading to faster consumption of that nutrient, thus acting to return the nutrients to equal steady-state levels. Such a self-organized neutral state is an attractor of the chemostat dynamics [18] and, correspondingly, of the $c_{0}\ll K$ limit of the serial dilution model. Note that the coexistence steady state is not a single fixed point, but rather a degenerate manifold of possible solutions (details in Appendix C).

Thus, in the chemostat-limit of the cases shown in Figs. 1B and C all the species will coexist. Conversely, if the supply lies outside the convex hull, (e.g., if we swapped the positions of the leftmost species and the supply in Fig. 1B) the number of surviving species would be strictly less than the number of nutrients, consistent with competitive exclusion. To understand the convex-hull rule, note that a state of arbitrarily high coexistence can only occur if the chemostat self-organizes to a “neutral” state in which the nutrient concentrations are all equal, and thus all strategies have the same growth rate. This state is achieved if and only if the total enzyme abundances lie along the same vector as the nutrient supply, which is achievable only if the supply lies within the convex hull of the strategies present.

As in the chemostat model, the serial dilution model can support either coexistence or competitive exclusion. However, if one chooses system parameters near the transition between these two states, it requires a very large number of batches for the simulation to reach steady state. This is a manifestation of the well-known phenomenon of critical slowing down [39]. Though in principle, critical slowing down is not a simulation artifact and could manifest in similar real-world systems, we expect that a variety of factors outside our modeling framework would preclude observation of this critical behavior.

We define the serial dilutions “steady state” to be reached when the relative species abundances after the nutrients are depleted (time $t_{f}$ after starting the batch) scale with the relative abundances at the beginning of that batch (time 0), i.e. $\rho_{\sigma}(t_{f})=\frac{\rho_{0}+c_{0}}{\rho_{0}}\rho_{\sigma}(0)$. We can expand the implicit equation for the steady state to first order in $c_{0}/K$ (details can be found in Appendix B),

[TABLE]

Dividing both sides by $t_{f}$ and defining,

[TABLE]

we reach the $c_{0}/K\ll 1$ steady-state condition for the serial dilution system:

[TABLE]

Averaged over a batch, $s_{i}$ is the average rate that nutrient $i$ is supplied, and $\tilde{\delta}$ is the average rate that all the nutrients are supplied per unit inoculum biomass. In analogy to the chemostat model, one can think of $s_{i}$ as the rate nutrient $i$ is continuously supplied. Moreover, for a chemostat, the parameter $\tilde{\delta}$ which would be the rate all nutrients are continuously supplied per unit biomass, would need to equal the dilution rate of the chemostat, $\delta$, to maintain steady state. A detailed analysis of the effects of larger bolus size, to second order in $c_{0}/K$, can be found in Appendix  C.

Effect of total nutrient bolus on coexistence

In the chemostat limit, increasing the nutrient supply rate simply proportionally increases the steady-state population abundances. However, away from this limit we find that the magnitude of the nutrient bolus can qualitatively affect the steady-state outcome of serial dilutions. To understand this effect, we first consider a simple case of two nutrients and two species as depicted in Fig. 2A. The two species will coexist if each species is invasible by the other. In our example, we first determine the invasibility of species R (strategy indicated by red circle) by species with strategies lying to its left. To this end, we choose a nutrient supply and perform model serial dilutions until steady state is reached. For a particular finite bolus size, we find that for all supplies within the hatched region an infinitesimal inoculum of any species lying to the left of R will increase more than R during a batch, and therefore can invade R. Similarly, we determine the invasibility of species B (strategy indicated by blue circle) by any species with a strategy lying to its right, and find the second hatched region. The intersection of these hatched regions for which (1) B can invade R and (2) R can invade B is the supply interval of mutual invasibility where these two species will stably coexist. The coexistence interval is bounded by the red and blue triangles, and each of these coexistence boundaries is a unique property of its corresponding species. We call these species-specific boundaries remapped because they generally lie at different locations on the simplex than the strategies they originated from, with the extent of remapping depending on the nutrient bolus size. At a more technical level, the remapped boundary for a given species and bolus size is the nutrient supply for which, over the course of a batch, all nutrients are equally valuable and so, a species with any strategy can neutrally invade and persist. This equality of nutrient value is defined in terms of the Monod function integrals for each nutrient (for details see Appendix  B).

Since the remapped coexistence boundaries depend on the nutrient bolus size $c_{0}$, changing bolus size can qualitatively change the steady-state outcome of serial dilutions. Figure 2B-D depicts an example of how $c_{0}$ affects remapping, and the consequences for species coexistence. At low bolus size, $c_{0}\ll K$, corresponding to the chemostat limit, Fig. 2B (left) shows that all species present achieve steady-state coexistence. This follows because the nutrient supply (black diamond) lies inside the convex hull. When $c_{0}$ is increased to $c_{0}\approx K$ (Fig. 2C), the coexistence boundaries are remapped towards the center of the simplex (dashed arrow). In this example, the nutrient supply now lies outside the convex hull. This results in one winner species (the dark blue one nearest the supply), with all others decreasing exponentially from batch to batch. This loss of coexistence with increasing nutrient bolus size is reminiscent of Rosenzweig’s “paradox of enrichment” in predator-prey systems [40]. Strikingly, however, as bolus size is further increased to $c_{0}\gg K$, the coexistence boundaries are remapped back to their original positions, so that the nutrient supply once again lies within the convex hull, and so steady-state coexistence of all species is recovered.

What causes the remapping of the coexistence boundaries inwards as $c_{0}/K\rightarrow 1$? Let us consider a single species growing on two nutrients supplied in the same proportion as its strategy (i.e. the fraction of Nutrient 1 is equal to $\alpha_{\sigma,1}$). At low $c_{0}/K$, this marks the remapped coexistence boundary and both nutrients will be equally valuable over the course of a batch. The balance is achieved because the nutrient with a higher initial abundance is more rapidly exhausted, while the nutrient with lower initial abundance is consumed more slowly and is therefore available for a longer span of time. At small $c_{0}/K$, the more rapid initial consumption of the more abundant nutrient does not influence the consumption rate of the less abundant nutrient because the bolus size is small relative to the initial population size, so the population does not grow substantially during the batch. This changes as $c_{0}/K$ increases: the rapid initial consumption of the more abundant nutrient leads to a substantial increase in the total population. The remaining low initial abundance nutrient is now consumed more quickly and is less available to an invader with a more balanced enzyme strategy. The nutrients are no longer equally valuable on average, and remedying this requires a more equally balanced nutrient bolus. Thus, the remapped coexistence boundary moves inwards (see Appendix F Fig. 13). In essence, as $c_{0}/K$ increases it is more difficult for the invader to grow because the resident gains an “early-bird” advantage: its initial growth allows it to more effectively exhaust the nutrients.

Why does the coexistence boundary of a species map back to its original strategy in the limit of large bolus size, $c_{0}\gg K$? In this limit, the nutrient uptake functions in Eq. 1 will be saturated during almost the entire period of a batch. Each species will therefore consume nutrients strictly in proportion to its strategy $\alpha_{\sigma,i}$. For the case of two nutrients (e.g., as shown in Fig. 1B), if there is only a single species present then if the supply lies anywhere to the left of its strategy, at some time during the batch there will be some of Nutrient 2 remaining after the bulk of Nutrient 1 has been consumed. Thus a single species can be invaded by any strategy to its left, provided the supply also lies to its left. Similarly, a species can be invaded by any strategy to its right if the supply lies to its right. This is exactly the condition for the coexistence boundary of a species to coincide with its actual strategy (details in Appendix  D).

We have rationalized coexistence in our serial dilution model in terms of mutual invasibility, but have not explicitly stated the condition for an arbitrary number of species to coexist in steady state. In the chemostat, all species coexist when the concentrations of all nutrients are equal, implying the same growth rate for all strategies. However, for serial dilutions the nutrient concentrations are generally not equal and are not even constant in time. Instead, it is the integrated growth contribution of every nutrient that must be equal to allow for arbitrary coexistence. In the case of equal enzyme budgets ($\varepsilon=0$), this condition occurs when the time integrals of the nutrient Monod functions within a batch are all equal, i.e.,

[TABLE]

To understand this condition for coexistence beyond competitive exclusion, note that the instantaneous rate of growth of a species $\sigma$ is $\sum_{i}\alpha_{\sigma,i}c_{i}/(K_{i}+c_{i})$, so that the fold increase of a species during a batch is $\exp(\vec{\alpha}_{\sigma}\cdot\vec{I})$. This fold increase will be equal for all species if and only if Eq. 7 holds. When there are two nutrients, Eq. 7 holds at steady state whenever the supply is inside the convex hull of the coexistence boundaries of the species present (details in Appendix  B). For more nutrients, the corresponding condition is that the region of coexistence is bounded by contours that connect the outermost remapped nodes.

Given a fixed set of species and a choice of initial populations, repeating the growth-dilution batch procedure results in a steady state where the populations at the beginning of a batch do not change from batch to batch. The steady-state populations depend on the initial populations, with the set of all possible steady-state populations defining a coexistence manifold.

Steady-state diversity

As is apparent in Fig. 2C, not all strategies are remapped to the same extent. In Fig. 3A, we plot the remapping of coexistence boundaries as a function of nutrient bolus $c_{0}$. Note that: (i) the specialists $(0,1)$ and $(1,0)$ and the perfect generalist $(0.5,0.5)$ are not remapped at all; (ii) remapping is maximal for $c_{0}\approx K$; (iii) there is no remapping in both the $c_{0}\rightarrow 0$ and $c_{0}\rightarrow\infty$ limits (see also Appendix F Fig. 14). The extent of remapping also depends on the inoculum size $\rho_{0}$ as shown in Fig. 3B, which demonstrates that remapping is maximal for $\rho_{0}\ll K$ and vanishes for $\rho_{0}\gg K$.

How does nutrient bolus size influence steady-state species diversity? A useful summary statistic for quantifying diversity [41] is the effective number of species $m_{e}=e^{S}$ with the Shannon diversity $S=-\sum_{\sigma}P_{\sigma}\ln P_{\sigma}$ and $P_{\sigma}=\rho^{*}_{\sigma}(0)/\rho_{0}$, with $\rho^{*}_{\sigma}(0)$ the steady-state species abundances at the beginning of a batch. Diversity as measured by $m_{e}$ is shown in Fig. 3C for six choices of nutrient bolus composition. Notably, if the two nutrients are supplied equally (top curve, magenta), $m_{e}$ is independent of $c_{0}$ and coincides with the maximal possible diversity (dashed black line), namely equal steady-state abundances of all species (Fig. 3D, top). Conversely, if Nutrient 1 comprises only 5% of supplied nutrient (Fig. 3C, bottom curve, cyan), the number of effective species $m_{e}$ is lower than maximal even in the chemostat-limit of small bolus sizes $c_{0}\ll K$ and drops even further for $c_{0}\approx K$. This loss of diversity is due to the dramatically lowered steady-state abundances of strategies that favor Nutrient 1 (Fig. 3D, bottom). Two different effects underlie this change in community structure. The first is the early-bird effect described above: species specializing in more abundant nutrients gain a population advantage that allows them to rapidly consume less abundant nutrients that would otherwise support species with different enzyme specializations. The second effect is a well-known property of single nutrient competition and can be viewed as a “single-nutrient” early-bird effect. In this case, species that are superior competitors for a nutrient gain an exponential population advantage over inferior competitors, increasing their share of total nutrient beyond the ratio of initial consumption rates. Both of these effects increase in strength with larger bolus size because the early-bird advantage increases as growth proceeds. The combination of these effects results in the species specialized in consuming the most abundant nutrients consuming a larger fraction of all nutrients. However, for very large bolus sizes, saturation of nutrient uptake rates mitigates these two effects, leading to a lack of remapping for $c_{0}\gg K$ and diversity returning to its chemostat-limit. Though here we focused on the case of two nutrients, these results extend to more nutrients (for three nutrients see Appendix  F Fig. 15).

Models with fewer simplifying assumptions

In this final Results section, we consider the effects of relaxing some of our simplifying assumptions. We first assess the effect or different enzyme affinities, $K_{i}\neq K$, and different nutrient yields $Y_{i}\neq Y$. This is followed by a model that allows cross-feeding of metabolites. Finally, we consider populations bottlenecks and what happens when the fixed-enzyme-budget constraint is relaxed. We show that in all these cases, the dependence on bolus size can be understood as manifestations of the early-bird effect.

Unequal enzyme affinities and nutrient yields

We have thus far made the simplifying assumption that all enzymes have the same substrate affinity, such that $K_{i}\equiv K$. However, in nature different nutrients may have drastically different values of $K$. For example, the methanogen Methanosarcina barkeri has $K_{i}$ for the consumption of hydrogen and acetate that differ by approximately three orders of magnitude [42, 43]. How would such a large difference in the $K_{i}$ values impact diversity in our serial dilution ecosystem? In Fig. 4A we show diversity as a function of bolus size for a system with a large difference in $K_{i}$ ($K_{1}=10^{-3}$,$K_{2}=1$). Since the symmetry between nutrients is broken by the unequal $K_{i}$, we now show the entire range of nutrient proportions, not just the first half. In the chemostat limit, the diversity values are similar to those found in the system with equal $K_{i}$. This makes sense: in a chemostat the nutrients with higher $K_{i}$ can accumulate to higher levels to compensate for their slow consumption, leaving the steady-state behavior unchanged. However, outside the chemostat regime, differences in the $K_{i}$ have a drastic effect: when the nutrient with the lower $K_{i}$ is scant in supply, diversity generally increases with increasing $c_{0}$, while the opposite occurs when the nutrient with the lower $K_{i}$ is higher in supply.

We can understand these $K_{i}$-driven shifts in the nutrient-diversity relationship as due to changes in the identity of the early bird. In a model with equal $K_{i}$, the identity of the early bird is determined by which nutrient is more abundant: if the two nutrients have equal $K_{i}$, a species can gain an early-bird advantage by preferentially consuming the more abundant nutrient. This changes if one nutrient has a much lower $K_{i}$ than the other. In this case it may be advantageous to preferentially consume the nutrient with the lower $K_{i}$, even if it is the less abundant nutrient. If the nutrient with the lower $K_{i}$ is also the more abundant nutrient, this will intensify the early-bird advantage. Why does this change in the early bird’s identity change the form of the nutrient-diversity relationship? This change arises from a clash between optimal feeding behavior in the chemostat and seasonal regimes. In the chemostat, it is advantageous to focus on the most abundant nutrient, regardless of the value of $K_{i}$. Thus, in the chemostat limit, species focusing on the more abundant nutrient have an advantage. In the case of equal $K_{i}$ (or $K_{i}$ favoring the more abundant nutrient), this advantage is intensified by the early-bird effect, increasing the biomass of already abundant species and lowering diversity. By contrast, if the low abundance nutrient has a low $K_{i}$, the early-bird effect will have the opposite effect on diversity. Now the early-bird effect benefits species that were disadvantaged in the chemostat limit, leading to more equal abundances and higher diversity. This shift in abundances is shown in Fig. 4B. The change in the identity of the early bird can also explain more complex relationships between diversity and bolus size (see Appendix F Fig. 16).

In addition to unequal enzyme affinities, it is possible for different nutrients to have different yields, $Y_{i}$. In Fig. 4C we show the relationship between bolus and diversity for a system with $Y_{1}=10$ and $Y_{2}=1$. As expected, at low $c_{0}$ the diversity is similar to that in the case of equal $Y_{i}$. As $c_{0}/K$ increases, the diversity decreases initially and the symmetry-related bolus-composition cases (e.g. (0.2,0.8) and (0.8,0.2)) eventually diverge, with one’s diversity rising and the other’s continuing to fall. This behavior is explainable by the same logic as in the variable $K_{i}$ case: diversity rises or falls depending on whether the early-bird species was also favored in the chemostat limit. However, unlike the case of variable $K_{i}$, the diversity curves do not eventually return to the chemostat limit. Regardless of which nutrient the $Y_{i}$ favor, the diversity eventually begins decreasing monotonically as $c_{0}/K$ increases. This difference between the variable $K_{i}$ and variable $Y_{i}$ cases can be understood by considering what occurs when both nutrients are saturating. In the variable $K_{i}$ case, saturating nutrients are equal in value, implying a return to the chemostat limit as the early-bird effect weakens. In contrast, for variable $Y_{i}$, there remains a difference in the value of the two nutrients in the saturated regime, meaning that the early-bird effect will grow stronger and the early bird will take over the population. The beginning of this takeover can be seen at high bolus sizes in Fig. 4D. Note that for both variable $K_{i}$ and variable $Y_{i}$, these trends are also reflected in the remapping (see Appendix F Fig. 17).

Despite the large variation in relationships between diversity and bolus size, these phenomena can all be understood as consequences of the early-bird effect. As the model becomes more complex there are additional factors to consider in determining which nutrient will provide an early-bird advantage, but the fundamental mechanism of exploiting early growth advantages remains.

Cross-feeding

It is possible to extend Eqs. 2 and 3 beyond a single trophic layer, allowing for consumption of metabolic byproducts. This is a form of cross-feeding, which has generally been found to promote diversity [12] and stable community structure [24]. Here, cross-feeding is introduced through the byproduct matrix $\Gamma^{\sigma}_{i,i^{\prime}}$, which converts the consumption of nutrient $i^{\prime}$ to production of nutrient $i$ such that,

[TABLE]

In this framework, nutrient $i^{\prime}$ is converted to nutrient $i$ at no extra enzymatic cost, meaning that nutrient $i$ is simply a byproduct whenever nutrient $i^{\prime}$ is consumed for growth. (It would be straightforward to modify this framework so that nutrient conversion can be carried out independently from growth.) We focus on the simplest case: initially supplying only Nutrient 2, with Nutrient 1 solely derived as a metabolic byproduct via \Gamma^{\sigma}_{i,i^{\prime}}=\left(\begin{array}[]{cc}0&\Gamma\\ 0&0\end{array}\right) for all species. When $\Gamma=1$, upon consumption Nutrient 2 is perfectly converted to Nutrient 1, leading to an equal total supply of the two nutrients. More generally, $\int_{0}^{\infty}\sum_{\sigma}\rho_{\sigma}j_{\sigma,1}\,dt=\Gamma c_{2}(0)$ which allows a direct comparison between the unitrophic and bitrophic regimes: starting with $c_{2}(0)$ results in $(\Gamma+1)c_{2}(0)$ total nutrient, and hence the Nutrient 1 fraction is $\frac{\Gamma}{1+\Gamma}$ of the total.

How does cross-feeding influence diversity in our serial dilution model? In Fig. 5A we compare bitrophic diversity for six values of $\Gamma$ to their unitrophic equivalents (in Fig. 3C). We note that: (i) bitrophy still supports diversity greater than the competitive-exclusion limit; (ii) in the chemostat regime, $c_{0}\ll K$, the unitrophic and bitrophic schemes have identical values of $m_{e}$, and these drop as $c_{0}\rightarrow K$; (iii) but for bitrophy the $m_{e}$ does not recover for $c_{0}\gg K$; (iv) even when the total supply of both nutrients is equal ($\Gamma=1$), bitrophy leads to lower than maximal $m_{e}$ outside the chemostat limit. These features are clarified in Fig. 5B, which shows steady-state species abundances for $\Gamma$ values leading to a total Nutrient 1 supply fraction of 0.5 and 0.05, and highlights the lower diversity for bitrophy compared to unitrophy for large nutrient bolus size. This difference is due to an early-bird effect: the species consuming supplied nutrient early in the batch can build a sizable population before the competing species that rely on its byproduct. The early-bird population then outcompetes the others for byproduct consumption. As such, this effect increases with $c_{0}/\rho_{0}$. The effect also becomes stronger at low $c_{0}/K$, since this allows the early-bird species more time to grow before the byproduct accumulates to high enough levels to be significantly consumed (Appendix  F Fig. 18). Note that this effect is dependent on metabolite byproducts being also consumed by their producer. If the species in each trophic layer are single-nutrient specialists, then changes in $c_{0}/K$ have no impact on community diversity.

The behavior of the model with cross-feeding shows that the early-bird effect extends beyond simple metabolic trade-offs. More broadly, when species compete for multiple resources that are supplied in batches, a species’ survival depends on more than its ability to efficiently consume nutrients. An early-bird species, being more specialized in consuming the nutrients that are initially more abundant, gains a population advantage early in the batch. This population advantage may allow the early-bird species to out-compete other species even when consuming nutrients it is not specialized to consume. Despite its consumption inefficiencies, through sheer numbers the early-bird species can consume more of the remaining nutrients than its more specialized competitors.

Population bottlenecks

So far we have considered deterministic dynamics, which is appropriate for large populations. In natural settings, however, there are often small semi-isolated communities. For these communities, fluctuations can play an important role. In particular, population bottlenecks can lead to large demographic changes [44]. In our model, how does the nutrient supply affect diversity in such communities? To address this question, we applied discrete sampling of a finite population when diluting from one batch to the next (see Appendix Methods). With this protocol, an “extinction” occurs when sampling yields zero individuals of a species. For a long enough series of dilutions such extinctions would ultimately lead to near-complete loss of diversity. For small real-world populations, however, diversity may be maintained by migration. To model such migration we augmented the population at each dilution with a “spike-in” from a global pool of species, in the spirit of MacArthur’s theory of island biogeography [45]. Specifically, in the spike-in procedure, to prevent extinctions caused by sampling fluctuations, every new batch is inoculated with a small number of the original, founder species.

In Fig. 6A we show results of spike-in serial dilutions for a population bottleneck of 1008 cells. 95% of these cells are sampled from the previous batch, while 5% are sampled from a global pool, with equal abundances of 21 equally spaced strategies (cf. Fig. 3A). The resulting $m_{e}$ vs. $c_{0}$ curves have maximal $m_{e}$ for all six nutrient fractions in the regime $c_{0}\ll K$ where the 5% spike-in dominates sampling noise. As expected, for a balanced nutrient supply at any $c_{0}$, all species have the same average abundance (Fig. 6B top). By contrast, when Nutrient 1’s fraction is low (Fig. 6A cyan and 6B bottom), increasing $c_{0}$ increases the abundance gaps between the species, reflecting the uneven competition for Nutrient 2. Overall, the spike-in protocol leads to higher diversity at low $c_{0}$ than the deterministic case (starting from equal species abundances but with no spike-in, Fig. 3C). For large $c_{0}$, the $m_{e}$ vs. $c_{0}$ curves for these two protocols are indistinguishable. The only noticeable difference is that the spike-in maintains a higher level of the least competitive strains, but since these abundances are still low, this difference in not reflected in the $m_{e}$ values.

Unequal enzyme budgets

While we have assumed exact trade-offs to achieve diversity within a resource-competition model, the trade-offs present among real microorganisms will not be exact. For the serial dilution protocol with spike-ins, diversity is maintained by migration and so it is possible to relax the constraint of exact trade-offs. How does diversity depend on the nutrient supply if we allow species to have different enzyme budgets? We implemented random differences in species enzyme budgets by setting $\varepsilon=0.1$, i.e. a standard deviation of 10%, and plotted effective number of species $m_{e}$ in Fig. 6C. As in the $\varepsilon=0$ limit (Fig. 6A), at sufficiently small $c_{0}$ the spike-in procedure dominates both sampling noise and differential growth rates due to unequal enzyme budgets. Raising $c_{0}$ leads to a drop in $m_{e}$ (albeit still above the competitive-exclusion limit). Examining the species abundances in Fig. 6D, we note that differences in enzyme budget establish a fitness hierarchy even when nutrient fractions are equal (top), with those species with the highest budgets increasing in relative abundance as $c_{0}$ increases. The asterisk (*) marks the species with the highest total enzyme budget, which becomes the most abundant for $c_{0}>K$. Reducing Nutrient 1’s fraction to 0.05 results in a shifting abundance hierarchy (Fig. 6D, bottom): at low $c_{0}$ the highest abundance species is the one that consumes only Nutrient 2, as in the equivalent $\varepsilon=0$ case. However, increasing $c_{0}$ results in increased abundance for the species with the highest enzyme budget – which would ultimately lead to its domination for sufficiently large $c_{0}$. This increasing dominance of the species with the highest enzyme budget is another manifestation of the early-bird effect: as the amount of growth in a batch increases, the advantage of a larger enzyme budget further compounds. In short, for spike-in serial dilutions the influence of unequal enzyme budgets depends on the nutrient supply, such that the species with the largest budgets dominate for large, unbiased supplies.

Discussion

Natural ecosystems experience variations in the timing and magnitude of nutrient supply, and the impact of these variations on species diversity is not fully understood [27, 37]. To explore the impact of variable nutrient supply, we modeled resource competition in a serial dilution framework and analyzed the model’s steady states. We found that variable nutrient supply still allows for the high diversity seen in the continuous supply (“chemostat”) version of the model. Indeed, the serial dilution steady state mimics that of a chemostat when the amount of nutrients supplied in each batch is small. Surprisingly, however, supplying the nutrients as a bolus led to a dependence of diversity on the amount of supplied nutrients.

In contrast to existing literature on seasonality, we find that environmental fluctuations generally weaken coexistence in this model. This occurs as the result of an “early-bird” effect associated with supplying nutrients as large seasonal boluses instead of continuously. Some species can capitalize on rapid initial growth on an abundant nutrient to reach a large population size, which then allows them to deplete the remaining nutrients at the expense of their competitors. This early-bird effect restricts the range of environments in which communities can self-organize to a neutral state. We show that even when metabolic trade-offs are combined with stabilizing mechanisms, the impact of the early-bird effect remains. For example, in the case of cross-feeding, the community diversity falls as a function of $c_{0}/K$ due to the early-bird advantages gained by species at higher trophic levels.

While the idea of species gaining early advantages has been explored, such as in the literature on founder effects and speciation [46, 47], to the best of our knowledge this is the first demonstration of the influence of the early-bird effect on the diversity of seasonal ecosystems. We believe that this effect will occur in a variety of such ecosystems, as its only fundamental requirement is competition for multiple nutrients that are supplied in a time-dependent manner. Interestingly, while the early-bird effect plays a large role in our model, it is not the only bolus-dependent effect that influences diversity. We also observe another effect that can be viewed as a “single-nutrient” version of the early-bird effect. This effect arises from a well-studied property of competition: as growth proceeds, a superior competitor for a nutrient gains an exponential advantage over inferior competitors for that nutrient. Like the early-bird effect, this shifts the system’s biomass towards species more specialized in initially abundant nutrients, particularly for large but non-saturating nutrient bolus sizes. This single-nutrient effect can co-occur with the early-bird effect, for example in competition for an abundant nutrient between two early-bird species.

The form of seasonality we explore in this manuscript, where mixed boluses are supplied periodically, is only one possible form of seasonal nutrient supply. The impact of the early-bird effect and single-nutrient competition will likely differ between different forms of seasonality. For example, we show in Appendix F Fig. 11 that supplying cycles of single nutrient boluses that approach an equal distribution of nutrients results in lower diversity than supplying mixed equal nutrient boluses. While this form of seasonality differs from the one we characterized, we can still understand the loss of diversity as arising from the single-nutrient competition effect initially observed in our mixed-bolus models. We expect the principles gleaned from our models to be of use in understanding diversity in a variety of seasonal ecosystems.

Finding a general relation between the amount of nutrient supplied to a community and its diversity is a long-standing goal of theoretical ecology [28, 29, 30]. We found that in our model the form of the nutrient-diversity relation (NDR) can change based on model details. The model has two regimes: a low diversity and a high diversity regime. The former satisfies competitive exclusion (no more species coexisting than resources), whereas the latter exceeds competitive exclusion and occurs when the nutrient supply lies within the convex hull of the remapped metabolic strategies present [18]. At the bifurcation point between the two regimes, we observe critical slowing down in that the number of dilutions required to reach steady state diverges.

In the high diversity regime, the NDR can take several forms, resulting from the interplay of the early-bird effect and other mechanisms. Even with a single trophic layer, the NDR can be U-shaped, hump-shaped, monotonically decreasing, or have multiple peaks. These trends can then be further modified by the addition of more trophic layers, differences in enzyme budgets, etc.

Experimental studies that characterize the NDRs of microbial ecosystems have reached similarly variable conclusions. For example, one work studying bacterial communities in Arctic deep-sea sediments found an increasing trend between energy input and richness [34], while a study on photosynthetic microbial mats found a negative relationship between energy input and richness [35]. A meta-analysis of aquatic microbial ecosystems found examples of both monotonic and non-monotonic NDRs, with no single form dominating [37]. Our theoretical results, together with these experimental findings, indicate that there may be no single universal NDR in microbial ecosystems. This conclusion suggests that the best approach for characterizing the NDR of a given ecosystem is not to apply a one-size-fits-all theory, but to analyze the role of different factors such as cross-feeding, trade-offs, and immigration in determining that particular ecosystem’s NDR. While we have focused on microbial systems, the absence of a universal NDR is consistent with results from recent work in plants [33].

We found that the stringency of metabolic trade-offs has a large impact on community diversity. We imposed a metabolic enzyme budget on each species to reflect the reality that microbial cells have a finite capacity to synthesize proteins and must carefully apportion their proteome [48]. However, while it is true that microbes have limited biosynthetic capacity, it is unclear how strict are the resulting trade-offs. For this reason, we characterized versions of the model with both exact and inexact trade-offs. Our results show that the form of an ecosystem’s NDR can depend on the stringency of metabolic trade-offs. This finding is not exclusive to the serial dilution model. The stringency of trade-offs was also important in the original chemostat setting: in a birth-death-immigration framework, small violations of the enzyme budget still allowed for high levels of coexistence, but large violations disrupted coexistence [18]. These results suggest that an experimental characterization of the stringency of metabolic trade-offs among microbes would provide a valuable ecological parameter. Note that metabolic trade-offs are only one of the many types of trade-offs microbes are subject to; other types of trade-offs, such as constraints between biomass yield and growth rate [49], may also shape a community’s NDR.

In constructing a model, we made a number of assumptions about the way in which microbes consume and utilize nutrients. Some of these assumptions do not apply to all microbial communities, and the impact of relaxing these assumptions can affect the NDR. For example, we mostly focused on communities where all nutrients are equally valuable (i.e. $Y_{i}=Y_{j}\hskip 5.69054pt\forall i,j$). However, biomass yields can vary between nutrients and between species, which we explored in Fig. 4C-D. Notably, unequal yields create differences between nutrients even in the saturating regime ($c_{0}\gg K$), leading to a departure from the chemostat limit at large nutrient boluses. We show in Appendix F Fig.  17 that coexistence in the serial dilution model is robust to varying yield, as long as all species have the same yield on a given nutrient. The scenario where species have different biomass yields on the same nutrient is conceptually similar to the case of inexact trade-offs, since some species will have a strict advantage over others. Thus, it is likely that these unequal yields between species will lead to a reduction in community diversity. However, varying the yield in this manner also allows for the inclusion of new trade-offs that may impact diversity, such as the aforementioned trade-off between yield and growth rate [49]. We also explored the effects of unequal Monod constants for different nutrients (cf. Fig. 4A-B). We found that if a low-abundance nutrient also has a low $K_{i}$, the early-bird effect favors species that were disadvantaged in the chemostat limit, thus reversing the equal-$K_{i}$ NDR and leading to hump-shaped NDR curves. Indeed, large differences in $K_{i}$ values can lead to a multi-peaked NDR as shown in Appendix F Fig. 16.

Our model assumes that all nutrients are substitutable (i.e. only one of the multiple nutrients is required for growth). In real ecosystems, microbes can require multiple complementary nutrients to grow, e.g. sources of carbon, nitrogen, and phosphorus. In cases where one class of complementary nutrient is strongly limiting, a model with both complementary and substitutable resources would essentially reduce to the current model of only substitutable resources. This case is likely the more common one, e.g. as many soils are carbon limited [50, 51]. However, in cases where no single nutrient is strongly limiting, the presence of complementary nutrients would possibly lead to different NDRs, which will be an interesting direction for future study.

Our modeling predictions, e.g. the convex hull condition and the changes in diversity due to the early-bird effect, are in principle testable. To connect our modeling assumptions to real microbial systems, we compare our growth model of substitutable and simultaneous nutrient consumption to previously published experimental data from Escherichia coli growing in batch and chemostat conditions. We find that our modeling assumptions are consistent with both datasets and outline potential future experiments to test the model’s multispecies predictions, detailed in Appendix  E. As is apparent in Appendix E Fig.  9, the growth dynamics of E. coli at low nutrient levels is well described by our modeling framework. The experiments we compared to were performed with the same strain of E. coli, meaning that inclusion of different microbes would be needed to test the multispecies predictions. To determine the strategies of other microbes, including other strains of E. coli, the most practical approach would likely be batch culturing. Once strains with different strategies have been identified, nutrient-diversity relationships could then be obtained by competing strains in serial dilution culture and measuring the community diversity (e.g. via fluorescent tags or by 16S rRNA sequencing) as a function of the total concentration of multiple, substitutable nutrients provided at the start of each batch.

Acknowledgments

N.S.W. supervised the research; A.E. and J.G.L. wrote simulations and did analytic calculations; all authors interpreted results; A.E., J.G.L., and N.S.W. wrote the paper.

All code and data used in this manuscript available at:

https://github.com/AmirErez/SeasonalEcosystem.

Methods

This section describes the simulation methods used in this manuscript. All code and data used in this manuscript can be found at https://github.com/AmirErez/SeasonalEcosystem.

Deterministic dynamics

We numerically solve the ODEs within each batch using a custom MATLAB-coded fourth-order Runge-Kutta solver with adaptive step size. Step size at a given time step is chosen such that the relative change of all state variables is below a predetermined threshold.

Population bottleneck sampling

We implement discrete sampling when diluting from one batch to the next by picking without replacement $\rho_{0}$ individuals from a total end-of-batch population of $\rho_{0}+c_{0}$. If there are non-integer populations at the end of a batch (as can occur with deterministic dynamics), they are rounded up if $\rho_{\sigma}-\mbox{floor}(\rho_{\sigma})>U(0,1)$ where floor rounds down to the nearest integer and $U(0,1)$ is a uniform random variable between 0 and 1. For all simulations with stochastic bottlenecks, we allow the simulation to equilibrate for 10,000 dilutions and average over 10,000 further dilutions.

Appendix A General form of the model

The most general form of the model considered in this manuscript includes variable nutrient yield $Y_{i}$, Monod half-velocity constant $\tilde{K}_{i}$, and enzyme cost $w_{i}$:

[TABLE]

The enzyme costs $w_{i}$ and total enzyme budget $E$ of the original equations can be removed by rescaling the strategies and nutrient concentrations such that $\tilde{\alpha}_{\sigma i}=(E/w_{i})\alpha_{\sigma i}$ and $\tilde{c_{i}}=(E/w_{i})c_{i}$. This rescaling leads to a new effective Monod half-velocity constant and yield such that $\tilde{K}_{i}=(E/w_{i})K_{i}$ and $\tilde{Y}_{i}=(w_{i}/E)Y_{i}$. The simplified equations are therefore:

[TABLE]

A further rescaling with $c_{i}^{\prime}=c_{i}Y_{i}$, $K_{i}^{\prime}=K_{i}Y_{i}$, and $\alpha_{\sigma i}^{\prime}=\alpha_{\sigma i}Y_{i}$ reveals the impact of $Y_{i}$:

[TABLE]

Appendix B Mutual invasibility condition for coexistence beyond competitive exclusion

In our model, coexistence of an unlimited number of species can be traced back to the conditions for the coexistence of a smaller number of species. This is because in a system with $p$ nutrients, once $p$ species coexist they create an environment where all nutrients are equally valuable and all species can coexist. For example, understanding the conditions for unlimited coexistence in two nutrient competition requires us to examine the conditions that allow two species to coexist. In order for two species to coexist, they must be able to invade each other. This means that in an environment dominated by Species 1, Species 2 will have higher fitness and vice versa.

Under what nutrient supplies, $c_{i}(0)/c_{0}$, can two species invade each other? In the chemostat version of the model, these invasibility conditions are simple to determine. Consider two species where $\vec{\alpha}_{1}$ is to the left of $\vec{\alpha}_{2}$ on the 1-simplex. Species 2 can invade Species 1 if the nutrient supply is to the right of $\vec{\alpha}_{1}$. Species 1 can invade Species 2 if the nutrient supply is to the left of $\vec{\alpha}_{2}$. Therefore, the two species can mutually invade and coexist if and only if the nutrient supply lies between $\vec{\alpha}_{1}$ and $\vec{\alpha}_{2}$. This is precisely the convex hull condition, with no remapping.

For the same pair of species, how do we determine the nutrient supplies for which Species 1 can be invaded by Species 2 in the serial dilution version of the model? The fitness of a species in this model is the growth exponent $\sum_{i}\alpha_{\sigma,i}I_{i}$, meaning that Species 2 can invade Species 1 at nutrient supplies where Species 1 creates an environment such that $I_{2}>I_{1}$. The nutrient supply at which Species 1 creates an environment where $I_{1}=I_{2}$ therefore bounds the region of nutrient supplies for which Species 2 can invade. By the same logic, the border for the region where Species 1 can invade Species 2 is the nutrient supply at which Species 2 creates an environment where $I_{1}=I_{2}$. Therefore, the mutual invasibility region is now defined by the nutrient supplies where each species growing in isolation creates an environment where $I_{1}=I_{2}$. These points are what we refer to as the ”remapped coexistence boundaries” and, unlike in the chemostat version of the model, these generally do not correspond to the species’ strategies.

Appendix C Perturbation theory for $c_{0}/K\ll 1$

In the main text, we provide an explanation why in the limit of small nutrient bolus size $c_{0}/K$, the serial dilution model effectively becomes a chemostat. In this section, we prove this chemostat limit using a perturbation expansion to first order in $c_{0}/K$. Essentially, the Monod constant $K$ acts as the unit of nutrient and biomass in the system, which are measured in dimensionless units $c_{0}/K$ and $\rho_{0}/K$, respectively. Alternatively, one might choose to expand around a third ratio, $c_{0}/\rho_{0}$, which would be useful to model extremely nutrient-dilute conditions as found in some marine microbial ecosystems.

We define a perturbation expansion with respect to the small parameter $\phi=c_{0}/K$,

[TABLE]

We note that at $O(1)$ we have $\rho_{\sigma}(t)=\rho_{\sigma}(0)$ and $c_{i}(t)=0$ as expected. We begin by expanding the Monod function,

[TABLE]

Accordingly, in the kinetic equation for $c_{i}$,

[TABLE]

substituting the expansion in Eq. C and keeping the leading order, $c_{i}^{(1)}$, gives,

[TABLE]

We next solve for $\rho_{\sigma}^{(1)}$ using $c_{i}^{(1)}$, and then we will use $\rho_{\sigma}^{(1)}$ to solve for $c_{i}^{(2)}$. It is possible but not necessary for our purposes to iterate further. The kinetic equation for the biomass density $\rho_{\sigma}$ is,

[TABLE]

Substituting the $\phi$ expansion gives, to leading order,

[TABLE]

Taking the long-time limit, $t\gg\frac{1}{\gamma_{i}}$, we obtain,

[TABLE]

Focusing on the leading order, we conclude that,

[TABLE]

Substituting for $\gamma_{i}$ and for $c_{i}^{(1)}(0)=c_{i}(0)/\phi$, we have,

[TABLE]

Explicitly stating the batch number $d$, at the end of the batch, i.e. at time $t=t_{f}\gg\gamma_{i}^{-1}$, the biomass density is,

[TABLE]

In the serial dilution model with complete consumption of all nutrients $c_{0}$ and initial biomass $\rho_{0}$, the inoculum populations in batch $d+1$ can be computed from the populations at the time of complete nutrient consumption, $t_{f}$, in batch $d$,

[TABLE]

At steady state, we require that $\rho_{\sigma}(d+1,0)=\rho_{\sigma}(d,0)$:

[TABLE]

Our calculation, to order $\phi$, gives,

[TABLE]

Dividing both sides by $t_{f}$ and defining,

[TABLE]

we finally reach the $c_{0}/K\ll 1$ steady-state condition for the serial dilution system:

[TABLE]

Averaged over a batch, $s_{i}$ is the average rate that nutrient $i$ is supplied, and $\tilde{\delta}$ is the average rate that all the nutrients are supplied per unit inoculum biomass. If this were a chemostat rather than a serial dilution model, then one could think of $s_{i}$ as the rate nutrient $i$ is continuously supplied. Moreover, for a chemostat, the parameter $\tilde{\delta}$, which would be the rate all nutrients are continuously supplied per unit biomass, would need to equal $\delta$, the dilution rate of the chemostat to maintain steady state. Indeed, Eq. 32 is precisely the steady-state condition for the chemostat (Eq. 4 from Ref. [18]) with $s_{i}$ and $\tilde{\delta}=\delta$ interpreted as above.

Thus we complete the proof that in $c_{0}/K\ll 1$, the steady state of our serial dilution model is identical to the steady state of the equivalent chemostat model.

2nd-order corrections to remapping of the coexistence boundaries for $c_{0}/K\ll 1$

We have demonstrated above that the leading terms in an expansion for small nutrient supply retrieve the steady-state solution of the chemostat model. However, we know from numerical simulations, that as $c_{0}/K\equiv\phi$ is increased, the coexistence boundaries become remapped, away from the enzyme strategies. Since there is no remapping in the chemostat limit, equivalent to an expansion to order $\phi$ as proved above, to capture the remapping we expand to order $\phi^{2}$ in $c_{i}(t)$. To this end, we return to the $\phi$ expansion and extract the $\phi^{2}$ contribution,

[TABLE]

which gives,

[TABLE]

Now we can solve for the growth-function integrals $I_{i}$ for the case of a single species growing in isolation. We expand the growth-function integral,

[TABLE]

Substituting the order $\phi$ from Eq. C and integrating, gives:

[TABLE]

For a single species growing in isolation, $\gamma_{i}=\frac{1}{K}\alpha_{\sigma,i}\rho_{\sigma}(0)$. Thus, to order $\phi$, the chemostat limit, we obtain, (in the chemostat limit, for a single species),

[TABLE]

Thus, to satisfy the coexistence boundary conditions: $\forall i:I_{i}=\mbox{const}$, to order $\phi$ it must be that $\forall i:c_{i}(0)/\alpha_{\sigma,i}=\mbox{const}$. This is precisely the coexistence condition for the chemostat, explained in Appendix B. To obtain the remapping of the coexistence boundaries, we must expand $I_{i}$ to order $\phi^{2}$.

To order $\phi^{2}$ we substitute Eq. C for $c_{i}^{(1)}$ and Eq. 34 for $c_{i}^{(2)}$ and integrate, giving:

[TABLE]

which upon substituting $\gamma_{i}$ for a single species simplifies to:

[TABLE]

Collecting terms to order $\phi^{2}$ gives

[TABLE]

As before, the coexistence boundaries are defined by $I_{i}=\mbox{const}$. To order $\phi^{2}$, Eq. 40 can be used to solve for this remapping analytically. Eq. 40 also clarifies why perfect generalists ($\forall i,j:\alpha_{\sigma,i}=\alpha_{\sigma,j}$) do not get remapped, as stated in the main text. This is because for a generalist $\forall i,j:\alpha_{\sigma,j}^{2}/(\alpha_{\sigma,i}+\alpha_{\sigma,j})=\mbox{const.}$, meaning that the second-order $I_{i}$ will all be equal if the first-order $I_{i}$ are all the same. Moreover, the term $I_{j}^{(1)}\propto\rho_{\sigma}^{-1}(0)$ multiplies the order $\phi^{2}$ correction to the coexistence boundary, and as a result, the larger $\rho_{\sigma}(0)$, the smaller the remapping. A comparison between the analytic form of the remapping at small $c_{0}/K$ and its numerical form is shown in Appendix B Fig. 7. As is apparent, the agreement is excellent and extends to higher $c_{0}/K$ as $\rho_{0}/K$ increases because at high $\rho_{0}/K$ the remapping is small.

2nd-order corrections to the steady-state abundance manifold for $c_{0}/K\ll 1$

Intuitively, one expects that at steady state, the capacity to consume a nutrient will match the total amount of nutrient supplied. Indeed, this was previously shown in the chemostat, and we will extend this statement beyond the chemostat limit, to second order in $c_{0}/K$. First, we review the results for a chemostat with dilution rate $\delta$ and nutrient supply rate $s_{i}$ [18],

[TABLE]

We use the asterisk in $\rho^{*}_{\sigma^{\prime}}$ to make explicit that the abundances are the steady-state abundances, after the system has moved from its initial conditions. As a result, the steady state constrains the abundances $\rho_{\sigma}^{*}$, such that they must lie on a manifold of solutions that satisfy $\sum_{\sigma}\alpha_{\sigma,i}\rho^{*}_{\sigma}=\frac{E}{\delta}s_{i}$. This is precisely the requirement that the total nutrient consumption rate matches the nutrient supply rate.

As demonstrated earlier in this section, by identifying $\tilde{\delta}=\frac{c_{0}}{\rho_{0}t_{f}}$ and $s_{i}=\frac{c_{i}(0)}{t_{f}}$, to first order in $\phi\equiv c_{0}/K\ll 1$, the steady-state condition for the serial dilution system is identical to the steady state of the equivalent chemostat. Thus, to order $\phi$,

[TABLE]

In the serial-dilution framework, to leading order in $\phi$, the chemostat limit of the serial-dilutions steady state is,

[TABLE]

We note that $K\gamma_{i}$ from Eq. 43 is a solution of Eq. 42, with $\sum_{i}\alpha_{\sigma,i}=E$. Having established Eq. 43 as the leading order (chemostat limit) term in an expansion in $\phi$, we proceed to calculate $K\gamma_{i}$ to order $\phi^{2}$ to obtain corrections to the chemostat limit.

From Eq. 3, we obtain $\rho_{\sigma}(t)=\rho_{0}(0)e^{\sum_{i}\alpha_{\sigma,i}I_{i}}$, so that,

[TABLE]

At steady state, dilution at the end of the batch brings the system to its initial conditions, such that for some time $t_{f}$ when the nutrients in a batch have been consumed, the steady-state species abundances $\rho^{*}_{\sigma}(0)$ obey

[TABLE]

Indeed, $\rho^{*}_{\sigma}(0)$ cancels, yielding the serial-dilutions steady-state to order $\phi^{2}$,

[TABLE]

We have added a tilde sign, as in $\tilde{I}_{i}^{(1)}$, to stress that we use leading order in $\phi$, having already explicitly accounted for $\phi$ in the expansion, and so, take $\phi\tilde{I}_{i}^{(1)}=\frac{c_{0}}{\rho_{0}E}$. Substituting for $\tilde{I}_{i}^{(2)}$ similarly, reduces Eq. 46 to,

[TABLE]

Comparing Eq. C with Eq. 43 we note the small, order $c_{0}/\rho_{0}$ , corrections to the chemostat limit. We overlay the perturbation theory solution, Eq. C, on the numerical solution, showing good agreement for $c_{0}/K<1$, plotted in Appendix B Fig. 8. We note that for the case of balanced nutrients, $\forall i:c_{i}(0)=c_{0}/p$, we have $X_{i}=\frac{1}{2}$, and therefore, when the nutrient supply is balanced, there are no second order corrections to the chemostat solution. Moreover, for balanced supply, the exact numerical solution also does not deviate from the chemostat limit.

Appendix D Remapping of the coexistence boundaries for $c_{0}/K\gg 1$

Here we show that at high $c_{0}/K$ the coexistence boundaries remap to their chemostat positions. When a large nutrient bolus is present, the growth function is effectively always saturated such that

[TABLE]

where $\sum_{i}\alpha_{\sigma,i}=E$ is in units of 1/time, without loss of generality $E$ can be set to unity, but we keep it here to make the units explicit. Solving for $\rho(t)$ yields $\rho(t)=\rho_{0}e^{Et}$. The assumption that $c_{i}\gg K$ can then be applied to the nutrient dynamics, yielding:

[TABLE]

Solving for the nutrient dynamics leads to $c_{i}(t)=c_{i}(0)+\frac{\alpha_{\sigma,i}}{E}\rho_{0}(1-e^{Et})$. Since the growth function is nearly always saturated (giving an integrand value of 1), the growth-function integral $I_{i}=\int_{0}^{\infty}\frac{c_{i}}{K_{i}+c_{i}}\,dt$ approximately equals the time of nutrient exhaustion. Thus for a given nutrient $i$, the time when that nutrient is depleted $t_{i,f}$ is given by:

[TABLE]

Note that the coexistence boundaries are defined by $\forall i:I_{i}=\mbox{const}$ which is satisfied when the fraction of nutrients in the initial bolus matches the strategies,

[TABLE]

This is precisely the result in Appendix B, indicating that in the $c_{0}/K\gg 1$ limit the coexistence boundaries return to their chemostat values.

Appendix E Comparison of growth model assumptions with experimental data

We have explored the implications of our model in a variety of contexts, but our modeling framework drastically simplifies bacterial growth, ignoring many factors relevant for microbial coexistence, e.g., lag times for the recovery of growth and various responses to stresses including starvation. In this section we compare certain key aspects of our growth model assumptions with experimental data.

A well-known form of nutrient utilization in the microbiology literature is sequential utilization, where a preferred sugar (often glucose) is consumed before others [52]. However, this mechanism applies to high sugar levels (on the order of grams per liter), such as those found in laboratory media. Many natural environments, such as marine systems and feces, contain low concentrations of sugars [53, 54]. At such low concentrations, simultaneous utilization of multiple substitutable sugars is observed [38, 55, 56]. We therefore compared our modeling for single-species growth to previously published data from chemostat and batch experiments on E. coli supplied with multiple sugars at low concentrations.

We use chemostat data from Lendenmann et al. [56] who measured the steady-state concentrations of biomass and sugars, with E. coli continuously supplied with mixtures of glucose, fructose, and ribose. We applied the chemostat version of the model [18] and constrained the fit with previously measured values of the Monod constants $K_{i}$ for this strain [57]. From the fit, we estimated the consumption strategies $\alpha_{i}$ for glucose, fructose, and ribose, with the rest of the parameters being defined experimentally (see the end of this section for details of the fitting procedure). As shown in Appendix E Figs. 9A and B, the resulting model matches the data well with strategies, measured in $\textrm{(mg sugar)}\textrm{(mg biomass)}^{-1}\textrm{h}^{-1}$, of $\alpha_{\textrm{gluc}}=1.96\pm 0.12$, $\alpha_{\textrm{fruc}}=2.04\pm 0.11$, and $\alpha_{\textrm{ribo}}=1.41\pm 0.01$. This corresponds to a normalized strategy of (0.36, 0.38, 0.26). The only notable deviations between the best-fit model and the data occurs for two fructose steady states. These deviations would be corrected if $K_{\textrm{fruc}}$ was larger, suggesting that the $K_{\textrm{fruc}}$ used here may not reflect the actual value in the experiment. The model also accurately predicts the resulting steady-state biomass concentrations, which are a constant $47$ mg/L in the experiment and approximately constant at $45$ mg/L in our model. This agreement suggests that our growth model assumptions are consistent with the behavior of E. coli growing at low nutrient concentration with a continuous nutrient supply. Despite being supplied with a variety of different sugar mixtures, E. coli maintains a constant steady-state biomass in these experiments because all of the carbon sources are substitutable.

To explicitly test growth dynamics, though for a single species only, we compared our model to batch growth data from Egli et al. [55]. In this experiment, timecourses of biomass and nutrient concentrations were measured in a culture of E. coli supplied with a mixture of glucose and galactose. The E. coli used to seed this culture came from a glucose-limited chemostat (we also compared our model to a batch seeded with E. coli from a galactose-limited chemostat, see Appendix E Fig. 10). For this data we used our serial dilution model with Monod kinetics and $K_{i}$ values from measurements on the same strain [57]. As shown in Appendix E Fig. 9C, the agreement between the best-fit model and the experimental data is generally quite good over the entire time course. The estimated $\alpha_{i}$, measured in units of $\textrm{(mg sugar)}\textrm{(mg biomass)}^{-1}\textrm{h}^{-1}$, were $0.46\pm 0.04$ and $0.41\pm 0.03$ for glucose and galactose, respectively. The estimated yield was $0.42\pm 0.03$ $\textrm{(mg biomass)}\textrm{(mg sugar)}^{-1}$, similar to the experimentally measured yield used in the chemostat model of $0.45$ $\textrm{(mg biomass)}\textrm{(mg sugar)}^{-1}$. Our model captures the glucose and biomass trends very well, but some galactose data points fall outside of the confidence interval. In addition to possible experimental noise, this may be due to small variations in yield during growth, as a constant yield would imply that accurately modeling biomass and glucose would necessarily also accurately capture galactose ($c_{\textrm{gal}}(t)=\textrm{const}-\rho(t)-Yc_{\textrm{gluc}}(t)$). Qualitatively, the data support the assumptions of substitutable and simultaneous nutrient consumption, while the strong quantitative fit to our model supports the assumption that enzyme strategies do not vary significantly within a batch.

While we only compare our model to data from E. coli, substitutable and simultaneous growth on multiple nutrients has been observed in other bacteria such as Lactobacillus brevis [58], and has even been observed in non-prokaryotic organisms. For example, the eukaryote Kloeckera sp. 2201 has been shown to simultaneously utilize methanol and glucose as carbon sources [38]. Similarly, the methanogenic archaeon Methanosarcina barkeri can simultaneously utilize methanol and acetate in batch culture [59]. However, it should be noted that our simple model cannot describe the growth kinetics of all microbes in all conditions. For example, the inferred strategy of E. coli for glucose varied between the batch and chemostat experiments examined here, suggesting that the total metabolic enzyme budget of microbes changes in different conditions. Such variation is likely due to other cell functions, such as ribosome synthesis [60], consuming different fractions of the cell’s total material and energy budget, something we do not explicitly model. Indeed, our goal is not to precisely model all microbial growth phenomenon, but rather to construct a widely applicable approximation of microbial growth in order to better understand ecological dynamics.

Fitting to experimental data

To fit our model to the experimental data in Lendenmann et al. [56], we first digitally extracted the steady-state data points from the experimental figures. We used the model from Posfai et al. [18] with Monod kinetics. The $K_{i}$ of glucose, ribose, and fructose were taken as 73, 132, and 125 $\mu$g/L sugar, respectively [57]. The model was fit to the data and the standard error of parameters were estimated using the MATLAB curve fitting toolbox. The only parameters estimated were the $\alpha_{i}$ of glucose, ribose, and fructose. It was assumed that all sugars had a biomass yield of $Y=0.45$ $(\mbox{g biomass})(\mbox{g sugar})^{-1}$, as measured experimentally [56]. The supply rates for a given simulation were computed as $S_{i}=c_{f,i}\delta$, where $S_{i}$ is the nutrient supply rate of nutrient $i$, $c_{f,i}$ is the concentration in the feed of nutrient $i$, and $\delta$ is the dilution rate of the chemostat. The fitting process minimized only the sum of squared errors between the model and the nutrient concentration data, since steady-state biomass within the model is approximately constant and determined by measured parameters ($\rho_{ss}\approx\frac{Y\sum_{i}S_{i}}{\delta}$). Confidence intervals for parameters were estimated using MATLAB’s confint function, which computes the interval using an estimate of the diagonal elements of the covariance matrix of the coefficients multiplied by the inverse of the Student’s $t$ distribution. The prediction bounds (shaded regions) are calculated using MATLAB’s predint function, which uses the estimated covariance matrix and the Jacobian of the fitted values to the parameters to predict the bounds.

The data fitting procedure for the batch experimental data was similar to that employed for the chemostat experimental data. We digitally extracted the data points from the figure in Egli et al. [55] and used the MATLAB curve fitting toolbox. The biomass was reported as $\textrm{OD}_{546}$ and was converted to mg/L using a conversion factor measured for the same strain [56]. Two sugar data points that were taken before the first biomass measurement were removed so that the initial conditions of the system would be well-defined. We estimated the parameters of the serial dilution model developed in this paper assuming Monod kinetics (Eq. 1). It was further assumed that both sugars had the same yield, $Y$. The yield was not measured in the experimental study and was therefore left as a fitting parameter. The $K_{i}$ for glucose and galactose were 73 and 98 $\mu$g/L, respectively [57]. The three fitting parameters were the yield $Y$ and the strategies $\alpha_{i}$ for glucose and galactose. The data points of the sugar and biomass measurements were taken at slightly different times, so the effective biomass for the experimental data was obtained by linear interpolation of the data points. Confidence intervals and prediction bounds were estimated using the same methods as for the chemostat model.

Appendix F Supplemental Figures

In this section we present supplemental figures that support the main text. Each figure’s caption contains all pertinent information.