Abstract
Microbial communities feature an immense diversity of species and this diversity is linked to outcomes ranging from ecosystem stability to medical prognoses. Yet the mechanisms underlying microbial diversity are under debate. While simple resourcecompetition models don't allow for coexistence of a large number of species, it was recently shown that metabolic tradeoffs can allow unlimited diversity. Does this diversity persist with more realistic, intermittent nutrient supply? Here, we demonstrate theoretically that in serial dilution culture, metabolic tradeoffs allow for high diversity. When a small amount of nutrient is supplied to each batch, the serial dilution dynamics mimic a chemostatlike steady state. If more nutrient is supplied, community diversity shifts due to an 'earlybird' effect. The interplay of this effect with different environmental factors and diversitysupporting mechanisms leads to a variety of relationships between nutrient supply and diversity, suggesting that real ecosystems may not obey a universal nutrientdiversity relationship.
eLife digest
In most environments, organisms compete for limited resources. The number and relative abundance of species that an ecosystem can host is referred to as ‘species diversity’. The competitiveexclusion principle is a hypothesis which proposes that, in an ecosystem, competition for resources results in decreased diversity: only species best equipped to consume the available resources thrive, while their less successful competitors die off. However, many natural ecosystems foster a wide array of species despite offering relatively few resources.
Researchers have proposed many competing theories to explain how this paradox can emerge, but they have mainly focused on ecosystems where nutrients are steadily supplied. By contrast, less is known about the way species diversity is maintained when nutrients are only intermittently available, for example in ecosystems that have seasons.
To address this question, Erez, Lopez et al. modeled communities of bacteria in which nutrients were repeatedly added and then used up. Depending on conditions, a variety of relationships between the amount of nutrient supplied and community diversity could emerge, suggesting that ecosystems do not follow a simple, universal rule that dictates species diversity. In particular, the resulting communities displayed a higher diversity of microbes than the limit imposed by the competitiveexclusion principle.
Further observations allowed Erez, Lopez et al. to suggest guiding principles for when diversity in ecosystems will be maintained or lost. In this framework, ‘earlybird’ species, which rapidly use a subset of the available nutrients, grow to dominate the ecosystem. Even though ‘latebird’ species are more effective at consuming the remaining resources, they cannot compete with the increased sheer numbers of the ‘earlybirds’, leading to a ‘richgetricher’ phenomenon.
Oceanic plankton, arctic permafrost and many other threatened, resourcepoor ecosystems across the world can dramatically influence our daily lives. Closer to home, shifts in the microbe communities that live on the surface of the human body and in the digestive system are linked to poor health. Understanding how species diversity emerges and changes will help to protect our external and internal environments.
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 (LloydPrice et al., 2016; Ladau et al., 2013; Weigel and Pfister, 2019). These observations appear to clash with a prediction of some resourcecompetition models, known as the competitiveexclusion principle – namely, that steadystate coexistence is possible for only as many species as resources (Levin, 1970; Armstrong and McGehee, 1980). This conundrum is familiarly known as the ‘paradox of the plankton’ (Hutchinson, 1961). Solving this paradox may provide one key to predicting and controlling outcomes ranging from ecosystem stability to successful cancer treatments in humans (Ptacnik et al., 2008; van Elsas et al., 2012; Taur et al., 2014; Stein et al., 2013). Chesson, 2000 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 crossfeeding or antibiotic production and degradation (Goyal and Maslov, 2018; Kelsic et al., 2015), (ii) spatial heterogeneity (Murrell and Law, 2003; Tilman, 1994), (iii) persistent nonsteadystate dynamics (Hutchinson, 1961), and (iv) predation (Thingstad, 2000). Equalizing mechanisms have been studied through neutral theory, in which all species are assumed to have equal fitness (Hubbell, 2005), and recent work has proposed resourcecompetition models that selforganize to a neutral state (Posfai et al., 2017). 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 (Palmer, 1994). 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 (Chang, 2003) to the microbiota of some human populations (Smits et al., 2017). 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 (Chesson, 1994). But do these conclusions extend to equalizing mechanisms? To address this question, we consider a known resourcecompetition model that permits high diversity at steady state due to the equalizing effects of metabolic tradeoffs, which assume that microbes have a limited enzyme production capacity they must apportion. Here, we investigate the equalizing effect of metabolic tradeoffs 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 wellestablished as an experimental approach. For example, the bacterial populations in the Lenski longterm evolution experiment (Lenski and Travisano, 1994), experiments on community assembly (Goldford et al., 2018), and antibiotic crossprotection (Yurtsev et al., 2016) were all performed in serial dilution. While previous models of serial dilution have characterized competition between small numbers of species with tradeoffs in their growth characteristics (Stewart and Levin, 1973; Smith, 2011), the theoretical understanding of diversity in serial dilution is much less developed than for chemostatbased steadystate growth.
Here, we show that under serial dilution, metabolic tradeoffs can still support high diversity communities, but that this coexistence is now sensitive to environmental conditions. Interestingly, seasonality can both increase and decrease 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 ‘earlybird’ 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 (Tilman, 1982; Abrams, 1995; Leibold, 1996): 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 (Mittelbach et al., 2001; Waide et al., 1999; Adler et al., 2011), though recently there has been increased focus on microbial systems (Bienhold et al., 2012; Bernstein et al., 2017). In microbial experiments, some evidence has supported the ‘humpshaped’ unimodal trend predicted by many theories (Kassen et al., 2000). However, a metaanalysis by Smith, 2007 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 crossfeeding, 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 Figure 1A (see Appendix 1—table 1). 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,\mathrm{\dots},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 7—figure 1). 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 nutrientspecific 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}(\mathrm{\infty})\approx 0$ (we consider the case of incomplete depletion in Appendix 7—figure 2). Thus, at the end of a batch, the total biomass density of cells is ${\sum}_{\sigma}{\rho}_{\sigma}(\mathrm{\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 σ is defined by its unique enzyme strategy ${\overrightarrow{\alpha}}_{\sigma}=({\alpha}_{\sigma ,i},\mathrm{\dots},{\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 (KovárováKovar and Egli, 1998), though this assumption may not hold for all microbial species. Specifically, we assume that species σ consumes nutrient i at a rate ${j}_{\sigma ,i}$ (per unit biomass) that depends on nutrient availability c_{i} and on its enzymeallocation strategy ${\alpha}_{\sigma ,i}$ according to
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:
Since the level of one enzyme inevitably comes at the expense of another, we model this tradeoff via an approximately fixed total enzyme budget E. Formally, we take ${\sum}_{i}{\alpha}_{\sigma ,i}=E+\epsilon {\xi}_{\sigma}$, where ${\xi}_{\sigma}$ is a zeromean and unitvariance Gaussian variable. Without loss of generality we take $E=1$; initially we set $\epsilon =0$, which we call exact tradeoffs. This allows us to visualize the strategies ${\overrightarrow{\alpha}}_{\sigma}$ as points on a simplex, depicted as colored circles embedded in: (i) the interval $[0,1]$ for two nutrients (Figure 1B), or (ii) a triangle for three nutrients (Figure 1C), etc. One can plot the nutrient bolus composition ${c}_{i}/{c}_{0}$ on the same simplex, as depicted by the black diamonds in Figure 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 (Posfai et al., 2017).
For completeness, we now briefly describe the chemostat results from Posfai et al., 2017. In the presence of metabolic tradeoffs, the chemostat can support a higher species diversity than prescribed by the competitive exclusion principle as we demonstrate theoretically in Appendix 4. 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 Figure 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 steadystate levels. Such a selforganized neutral state is an attractor of the chemostat dynamics (Posfai et al., 2017) 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 4).
Thus, in the chemostatlimit of the cases shown in Figure 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 Figure 1B) the number of surviving species would be strictly less than the number of nutrients, consistent with competitive exclusion. To understand the convexhull rule, note that a state of arbitrarily high coexistence can only occur if the chemostat selforganizes 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 wellknown phenomenon of critical slowing down (Dakos and Bascompte, 2014). Though in principle, critical slowing down is not a simulation artifact and could manifest in similar realworld 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), that is, ${\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 4),
Dividing both sides by t_{f} and defining,
we reach the ${c}_{0}/K\ll 1$ steadystate condition for the serial dilution system:
Averaged over a batch, s_{i} is the average rate that nutrient i is supplied, and $\stackrel{~}{\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 $\stackrel{~}{\delta}$ which would be the rate all nutrients are continuously supplied per unit biomass, would need to equal the dilution rate of the chemostat, δ, 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 4.
Effect of total nutrient bolus on coexistence
In the chemostat limit, increasing the nutrient supply rate simply proportionally increases the steadystate population abundances. However, away from this limit we find that the magnitude of the nutrient bolus can qualitatively affect the steadystate outcome of serial dilutions. To understand this effect, we first consider a simple case of two nutrients and two species as depicted in Figure 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 speciesspecific 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 3).
Since the remapped coexistence boundaries depend on the nutrient bolus size c_{0}, changing bolus size can qualitatively change the steadystate 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, Figure 2B (left) shows that all species present achieve steadystate coexistence. This follows because the nutrient supply (black diamond) lies inside the convex hull. When c_{0} is increased to ${c}_{0}\approx K$ (Figure 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 predatorprey systems (Rosenzweig, 1971). 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 steadystate coexistence of all species is recovered.
What causes the remapping of the coexistence boundaries inwards as ${c}_{0}/K\to 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 7—figure 3). In essence, as ${c}_{0}/K$ increases it is more difficult for the invader to grow because the resident gains an ‘earlybird’ 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 Equation 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 Figure 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 5).
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 ($\epsilon =0$), this condition occurs when the time integrals of the nutrient Monod functions within a batch are all equal, that is,
To understand this condition for coexistence beyond competitive exclusion, note that the instantaneous rate of growth of a species σ is ${\sum}_{i}{\alpha}_{\sigma ,i}{c}_{i}/({K}_{i}+{c}_{i})$, so that the fold increase of a species during a batch is $\mathrm{exp}({\overrightarrow{\alpha}}_{\sigma}\cdot \overrightarrow{I})$. This fold increase will be equal for all species if and only if Equation 7 holds. When there are two nutrients, Equation 7 holds at steady state whenever the supply is inside the convex hull of the coexistence boundaries of the species present (details in Appendix 3). 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 growthdilution batch procedure results in a steady state where the populations at the beginning of a batch do not change from batch to batch. The steadystate populations depend on the initial populations, with the set of all possible steadystate populations defining a coexistence manifold.
Steadystate diversity
As is apparent in Figure 2C, not all strategies are remapped to the same extent. In Figure 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}\to 0$ and ${c}_{0}\to \mathrm{\infty}$ limits (see also Appendix 7—figure 4). The extent of remapping also depends on the inoculum size ${\rho}_{0}$ as shown in Figure 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 steadystate species diversity? A useful summary statistic for quantifying diversity (Jost, 2006) is the effective number of species ${m}_{e}={e}^{S}$ with the Shannon diversity $S={\sum}_{\sigma}{P}_{\sigma}\mathrm{ln}{P}_{\sigma}$ and ${P}_{\sigma}={\rho}_{\sigma}^{*}(0)/{\rho}_{0}$, with ${\rho}_{\sigma}^{*}(0)$ the steadystate species abundances at the beginning of a batch. Diversity as measured by m_{e} is shown in Figure 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 steadystate abundances of all species (Figure 3D, top). Conversely, if Nutrient 1 comprises only 5% of supplied nutrient (Figure 3C, bottom curve, cyan), the number of effective species m_{e} is lower than maximal even in the chemostatlimit 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 steadystate abundances of strategies that favor Nutrient 1 (Figure 3D, bottom). Two different effects underlie this change in community structure. The first is the earlybird 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 wellknown property of single nutrient competition and can be viewed as a ‘singlenutrient’ earlybird 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 earlybird 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 chemostatlimit. Though here we focused on the case of two nutrients, these results extend to more nutrients (for three nutrients see Appendix 7—figure 5).
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}\ne K$, and different nutrient yields ${Y}_{i}\ne Y$. This is followed by a model that allows crossfeeding of metabolites. Finally, we consider population bottlenecks and what happens when the fixedenzymebudget constraint is relaxed. We show that in all these cases, the dependence on bolus size can be understood as manifestations of the earlybird 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 (Robinson and Tiedje, 1984; Wandrey and Aivasidis, 1983). How would such a large difference in the K_{i} values impact diversity in our serial dilution ecosystem? In Figure 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 steadystate 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 nutrientdiversity 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 earlybird 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 earlybird advantage. Why does this change in the early bird’s identity change the form of the nutrientdiversity 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 earlybird effect, increasing the biomass of already abundant species and lowering diversity. By contrast, if the low abundance nutrient has a low K_{i}, the earlybird effect will have the opposite effect on diversity. Now the earlybird 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 Figure 4B. The change in the identity of the early bird can also explain more complex relationships between diversity and bolus size (see Appendix 7—figure 6).
In addition to unequal enzyme affinities, it is possible for different nutrients to have different yields, Y_{i}. In Figure 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 symmetryrelated boluscomposition 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 earlybird 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 earlybird 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 earlybird 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 Figure 4D. Note that for both variable K_{i} and variable Y_{i}, these trends are also reflected in the remapping (see Appendix 7—figure 7).
Despite the large variation in relationships between diversity and bolus size, these phenomena can all be understood as consequences of the earlybird effect. As the model becomes more complex there are additional factors to consider in determining which nutrient will provide an earlybird advantage, but the fundamental mechanism of exploiting early growth advantages remains.
Crossfeeding
It is possible to extend Equations 2 and 3 beyond a single trophic layer, allowing for consumption of metabolic byproducts. This is a form of crossfeeding, which has generally been found to promote diversity (Goyal and Maslov, 2018) and stable community structure (Goldford et al., 2018). Here, crossfeeding is introduced through the byproduct matrix ${\mathrm{\Gamma}}_{i,{i}^{\prime}}^{\sigma}$, which converts the consumption of nutrient ${i}^{\prime}$ to production of nutrient i such that,
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 ${\mathrm{\Gamma}}_{i,{i}^{\prime}}^{\sigma}=\left(\begin{array}{cc}\hfill 0\hfill & \hfill \mathrm{\Gamma}\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right)$ for all species. When $\mathrm{\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}^{\mathrm{\infty}}{\sum}_{\sigma}{\rho}_{\sigma}{j}_{\sigma ,1}dt=\mathrm{\Gamma}{c}_{2}(0)$ which allows a direct comparison between the unitrophic and bitrophic regimes: starting with ${c}_{2}(0)$ results in $(\mathrm{\Gamma}+1){c}_{2}(0)$ total nutrient, and hence the Nutrient 1 fraction is $\frac{\mathrm{\Gamma}}{1+\mathrm{\Gamma}}$ of the total.
How does crossfeeding influence diversity in our serial dilution model? In Figure 5A we compare bitrophic diversity for six values of $\mathrm{\Gamma}$ to their unitrophic equivalents (in Figure 3C). We note that: (i) bitrophy still supports diversity greater than the competitiveexclusion 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}\to 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 ($\mathrm{\Gamma}=1$), bitrophy leads to lower than maximal m_{e} outside the chemostat limit. These features are clarified in Figure 5B, which shows steadystate species abundances for $\mathrm{\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 earlybird 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 earlybird 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$ (with constant ${c}_{0}/{\rho}_{0}$), since this allows the earlybird species more time to grow before the byproduct accumulates to high enough levels to be significantly consumed (Appendix 7—figure 8). Note that this effect is dependent on metabolite byproducts being also consumed by their producer. If the species in each trophic layer are singlenutrient specialists, then changes in ${c}_{0}/K$ have no impact on community diversity.
The behavior of the model with crossfeeding shows that the earlybird effect extends beyond simple metabolic tradeoffs. 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 earlybird 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 earlybird species to outcompete other species even when consuming nutrients it is not specialized to consume. Despite its consumption inefficiencies, through sheer numbers the earlybird 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 semiisolated communities. For these communities, fluctuations can play an important role. In particular, population bottlenecks can lead to large demographic changes (Abel et al., 2015). 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 1). 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 nearcomplete loss of diversity. For small realworld populations, however, diversity may be maintained by migration. To model such migration we augmented the population at each dilution with a ‘spikein’ from a global pool of species, in the spirit of MacArthur’s theory of island biogeography (MacArthur and Wilson, 2001). Specifically, in the spikein procedure, to prevent extinctions caused by sampling fluctuations, every new batch is inoculated with a small number of the original, founder species.
In Figure 6A we show results of spikein 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. Figure 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% spikein dominates sampling noise. As expected, for a balanced nutrient supply at any c_{0}, all species have the same average abundance (Figure 6B top). By contrast, when Nutrient 1’s fraction is low (Figure 6A cyan and 6B bottom), increasing c_{0} increases the abundance gaps between the species, reflecting the uneven competition for Nutrient 2. Overall, the spikein protocol leads to higher diversity at low c_{0} than the deterministic case (starting from equal species abundances but with no spikein, Figure 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 spikein 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 tradeoffs to achieve diversity within a resourcecompetition model, the tradeoffs present among real microorganisms will not be exact. For the serial dilution protocol with spikeins, diversity is maintained by migration and so it is possible to relax the constraint of exact tradeoffs. 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 $\epsilon =0.1$, that is, a standard deviation of 10%, and plotted effective number of species m_{e} in Figure 6C. As in the $\epsilon =0$ limit (Figure 6A), at sufficiently small c_{0} the spikein 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 competitiveexclusion limit). Examining the species abundances in Figure 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}\gg K$. Reducing Nutrient 1’s fraction to 0.05 results in a shifting abundance hierarchy (Figure 6D, bottom): at low c_{0} the highest abundance species is the one that consumes only Nutrient 2, as in the equivalent $\epsilon =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 earlybird effect: as the amount of growth in a batch increases, the advantage of a larger enzyme budget further compounds. In short, for spikein 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 (Smith, 2011; Smith, 2007). 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 can both weaken and strengthen coexistence in this model. This occurs as the result of an ‘earlybird’ 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 earlybird effect can both restrict and expand the range of environments in which communities can selforganize to a neutral state. We show that even when metabolic tradeoffs are combined with stabilizing mechanisms, the impact of the earlybird effect remains. For example, in the case of crossfeeding, the community diversity falls as a function of ${c}_{0}/K$ due to the earlybird 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 (Barton and Charlesworth, 1984; Brown,, 1957), to the best of our knowledge this is the first demonstration of the influence of the earlybird 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 timedependent manner. Interestingly, while the earlybird effect plays a large role in our model, it is not the only bolusdependent effect that influences diversity. We also observe another effect that can be viewed as a ‘singlenutrient’ version of the earlybird effect. This effect arises from a wellstudied property of competition: as growth proceeds, a superior competitor for a nutrient gains an exponential advantage over inferior competitors for that nutrient. Like the earlybird effect, this shifts the system’s biomass towards species more specialized in initially abundant nutrients, particularly for large but nonsaturating nutrient bolus sizes. This singlenutrient effect can cooccur with the earlybird effect, for example in competition for an abundant nutrient between two earlybird 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 earlybird effect and singlenutrient competition will likely differ between different forms of seasonality. For example, we show in Appendix 7—figure 1 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 singlenutrient competition effect initially observed in our mixedbolus 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 longstanding goal of theoretical ecology (Tilman, 1982; Abrams, 1995; Leibold, 1996). We found that in our model the form of the nutrientdiversity 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 (Posfai et al., 2017). 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 earlybird effect and other mechanisms. Even with a single trophic layer, the NDR can be Ushaped, humpshaped, 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 deepsea sediments found an increasing trend between energy input and richness (Bienhold et al., 2012), while a study on photosynthetic microbial mats found a negative relationship between energy input and richness (Bernstein et al., 2017). A metaanalysis of aquatic microbial ecosystems found examples of both monotonic and nonmonotonic NDRs, with no single form dominating (Smith, 2007). 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 onesizefitsall theory, but to analyze the role of different factors such as crossfeeding, tradeoffs, 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 (Adler et al., 2011).
We found that the stringency of metabolic tradeoffs 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 (Basan et al., 2015). However, while it is true that microbes have limited biosynthetic capacity, it is unclear how strict are the resulting tradeoffs. For this reason, we characterized versions of the model with both exact and inexact tradeoffs. Our results show that the form of an ecosystem’s NDR can depend on the stringency of metabolic tradeoffs. This finding is not exclusive to the serial dilution model. The stringency of tradeoffs was also important in the original chemostat setting: in a birthdeathimmigration framework, small violations of the enzyme budget still allowed for high levels of coexistence, but large violations disrupted coexistence (Posfai et al., 2017). These results suggest that an experimental characterization of the stringency of metabolic tradeoffs among microbes would provide a valuable ecological parameter. Note that metabolic tradeoffs are only one of the many types of tradeoffs microbes are subject to; other types of tradeoffs, such as constraints between biomass yield and growth rate (Wortel et al., 2018), 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}\forall i,j$). However, biomass yields can vary between nutrients and between species, which we explored in Figure 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. 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 tradeoffs, 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 tradeoffs that may impact diversity, such as the aforementioned tradeoff between yield and growth rate (Wortel et al., 2018). We also explored the effects of unequal Monod constants for different nutrients (cf. Figure 4A–B). We found that if a lowabundance nutrient also has a low K_{i}, the earlybird effect favors species that were disadvantaged in the chemostat limit, thus reversing the equalK_{i} NDR and leading to humpshaped NDR curves. Indeed, large differences in K_{i} values can lead to a multipeaked NDR as shown in Appendix 7—figure 6.
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 (Aldén et al., 2001; Demoling et al., 2007). 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 earlybird 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 6. As is apparent in Appendix 6—figure 1, the growth dynamics of E. coli at low nutrient levels is well described by our modeling framework. The experiments we compared 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, nutrientdiversity 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.
Appendix 1
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 (Erez, 2020; copy archived at https://github.com/elifesciencespublications/SeasonalEcosystem).
Deterministic dynamics
We numerically solve the ODEs within each batch using a custom MATLABcoded fourthorder RungeKutta 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 endofbatch population of ${\rho}_{0}+{c}_{0}$. If there are noninteger populations at the end of a batch (as can occur with deterministic dynamics), they are rounded up if ${\rho}_{\sigma}\text{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 2
General form of the model
The most general form of the model considered in this manuscript includes variable nutrient yield Y_{i}, Monod halfvelocity constant ${\stackrel{~}{K}}_{i}$, and enzyme cost w_{i}:
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 ${\stackrel{~}{\alpha}}_{\sigma i}=(E/{w}_{i}){\alpha}_{\sigma i}$ and $\stackrel{~}{{c}_{i}}=(E/{w}_{i}){c}_{i}$. This rescaling leads to a new effective Monod halfvelocity constant and yield such that ${\stackrel{~}{K}}_{i}=(E/{w}_{i}){K}_{i}$ and ${\stackrel{~}{Y}}_{i}=({w}_{i}/E){Y}_{i}$. The simplified equations are therefore:
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}:
Appendix 3
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 ${\overrightarrow{\alpha}}_{1}$ is to the left of ${\overrightarrow{\alpha}}_{2}$ on the 1simplex. Species 2 can invade Species 1 if the nutrient supply is to the right of ${\overrightarrow{\alpha}}_{1}$. Species 1 can invade Species 2 if the nutrient supply is to the left of ${\overrightarrow{\alpha}}_{2}$. Therefore, the two species can mutually invade and coexist if and only if the nutrient supply lies between ${\overrightarrow{\alpha}}_{1}$ and ${\overrightarrow{\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 4
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 nutrientdilute conditions as found in some marine microbial ecosystems.
We define a perturbation expansion with respect to the small parameter $\varphi ={c}_{0}/K$,
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,
Accordingly, in the kinetic equation for c_{i},
substituting the expansion in Equation 19 and keeping the leading order, ${c}_{i}^{(1)}$, gives,
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,
Substituting the $\varphi $ expansion gives, to leading order,
Taking the longtime limit, $t\gg \frac{1}{{\gamma}_{i}}$, we obtain,
Focusing on the leading order, we conclude that,
Substituting for ${\gamma}_{i}$ and for ${c}_{i}^{(1)}(0)={c}_{i}(0)/\varphi $, we have,
Explicitly stating the batch number d, at the end of the batch, that is, at time $t={t}_{f}\gg {\gamma}_{i}^{1}$, the biomass density is,
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,
At steady state, we require that ${\rho}_{\sigma}(d+1,0)={\rho}_{\sigma}(d,0)$:
Our calculation, to order $\varphi $, gives,
Dividing both sides by t_{f} and defining,
we finally reach the ${c}_{0}/K\ll 1$ steadystate condition for the serial dilution system:
Averaged over a batch, s_{i} is the average rate that nutrient i is supplied, and $\stackrel{~}{\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 $\stackrel{~}{\delta}$, which would be the rate all nutrients are continuously supplied per unit biomass, would need to equal δ, the dilution rate of the chemostat to maintain steady state. Indeed, Equation 32 is precisely the steadystate condition for the chemostat (Equation 4 from Posfai et al., 2017) with s_{i} and $\stackrel{~}{\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.
Secondorder corrections to remapping of the coexistence boundaries for ${c}_{\mathrm{0}}\mathrm{/}K\mathrm{\ll}\mathrm{1}$
We have demonstrated above that the leading terms in an expansion for small nutrient supply retrieve the steadystate solution of the chemostat model. However, we know from numerical simulations, that as ${c}_{0}/K\equiv \varphi $ 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 $\varphi $ as proved above, to capture the remapping we expand to order ${\varphi}^{2}$ in ${c}_{i}(t)$. To this end, we return to the $\varphi $ expansion and extract the ${\varphi}^{2}$ contribution,
which gives,
Now we can solve for the growthfunction integrals I_{i} for the case of a single species growing in isolation. We expand the growthfunction integral,
Substituting the order $\varphi $ from Equation 22 and integrating, gives:
For a single species growing in isolation, ${\gamma}_{i}=\frac{1}{K}{\alpha}_{\sigma ,i}{\rho}_{\sigma}(0)$. Thus, to order $\varphi $, the chemostat limit, we obtain, (in the chemostat limit, for a single species),
Thus, to satisfy the coexistence boundary conditions: $\forall i:{I}_{i}=\text{const}$, to order $\varphi $ it must be that $\forall i:{c}_{i}(0)/{\alpha}_{\sigma ,i}=\text{const}$. This is precisely the coexistence condition for the chemostat, explained in Appendix 3. To obtain the remapping of the coexistence boundaries, we must expand I_{i} to order ${\varphi}^{2}$.
To order ${\varphi}^{2}$ we substitute Equation 22 for ${c}_{i}^{(1)}$ and Equation 34 for ${c}_{i}^{(2)}$ and integrate, giving:
which upon substituting ${\gamma}_{i}$ for a single species simplifies to:
Collecting terms to order ${\varphi}^{2}$ gives
As before, the coexistence boundaries are defined by ${I}_{i}=\text{const}$. To order ${\varphi}^{2}$, Equation 40 can be used to solve for this remapping analytically. Equation 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})=\text{const.}$, meaning that the secondorder I_{i} will all be equal if the firstorder I_{i} are all the same. Moreover, the term ${I}_{j}^{(1)}\propto {\rho}_{\sigma}^{1}(0)$ multiplies the order ${\varphi}^{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 4—figure 1. 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.
Secondorder corrections to the steadystate abundance manifold for ${c}_{\mathrm{0}}\mathrm{/}K\mathrm{\ll}\mathrm{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 δ and nutrient supply rate s_{i} (Posfai et al., 2017),
We use the asterisk in ${\rho}_{{\sigma}^{\prime}}^{*}$ to make explicit that the abundances are the steadystate 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 $\stackrel{~}{\delta}=\frac{{c}_{0}}{{\rho}_{0}{t}_{f}}$ and ${s}_{i}=\frac{{c}_{i}(0)}{{t}_{f}}$, to first order in $\varphi \equiv {c}_{0}/K\ll 1$, the steadystate condition for the serial dilution system is identical to the steady state of the equivalent chemostat. Thus, to order $\varphi $,
In the serialdilution framework, to leading order in $\varphi $, the chemostat limit of the serialdilutions steady state is,
We note that $K{\gamma}_{i}$ from Equation 43 is a solution of Equation 42, with ${\sum}_{i}{\alpha}_{\sigma ,i}=E$. Having established Equation 43 as the leading order (chemostat limit) term in an expansion in $\varphi $, we proceed to calculate $K{\gamma}_{i}$ to order ${\varphi}^{2}$ to obtain corrections to the chemostat limit.
From Equation 3, we obtain ${\rho}_{\sigma}(t)={\rho}_{0}(0){e}^{{\sum}_{i}{\alpha}_{\sigma ,i}{I}_{i}}$, so that,
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 steadystate species abundances ${\rho}_{\sigma}^{*}(0)$ obey
Indeed, ${\rho}_{\sigma}^{*}(0)$ cancels, yielding the serialdilutions steadystate to order ${\varphi}^{2}$,
We have added a tilde sign, as in ${\stackrel{~}{I}}_{i}^{(1)}$, to stress that we use leading order in $\varphi $, having already explicitly accounted for $\varphi $ in the expansion, and so, take $\varphi {\stackrel{~}{I}}_{i}^{(1)}=\frac{{c}_{0}}{{\rho}_{0}E}$. Substituting for ${\stackrel{~}{I}}_{i}^{(2)}$ similarly, reduces Equation 46 to,
Comparing Equation 47 with Equation 43 we note the small, order ${c}_{0}/{\rho}_{0}$ , corrections to the chemostat limit. We overlay the perturbation theory solution, Equation 47, on the numerical solution, showing good agreement for ${c}_{0}/K<1$, plotted in Appendix 4—figure 2. 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 5
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
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:
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 growthfunction integral ${I}_{i}={\int}_{0}^{\mathrm{\infty}}\frac{{c}_{i}}{{K}_{i}+{c}_{i}}\mathit{d}t$ 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:
Note that the coexistence boundaries are defined by $\forall i:{I}_{i}=\text{const}$ which is satisfied when the fraction of nutrients in the initial bolus matches the strategies,
This is precisely the result in Appendix 3, indicating that in the ${c}_{0}/K\gg 1$ limit the coexistence boundaries return to their chemostat values.
Appendix 6
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 wellknown form of nutrient utilization in the microbiology literature is sequential utilization, where a preferred sugar (often glucose) is consumed before others (Monod, 1942). 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 (Münster, 1993; Flourie et al., 1986). At such low concentrations, simultaneous utilization of multiple substitutable sugars is observed (KovárováKovar and Egli, 1998; Egli et al., 1993; Lendenmann et al., 1996). We therefore compared our modeling for singlespecies 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., 1996 who measured the steadystate 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 (Posfai et al., 2017) and constrained the fit with previously measured values of the Monod constants K_{i} for this strain (Lendenmann and Egli, 1998). 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 6—figures 1A and B, the resulting model matches the data well with strategies, measured in $(\mathrm{m}\mathrm{g}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{s}\mathrm{u}\mathrm{g}\mathrm{a}\mathrm{r})(\mathrm{m}\mathrm{g}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s})}^{1}{\mathrm{h}}^{1$, of ${\alpha}_{\text{gluc}}=1.96\pm 0.12$, ${\alpha}_{\text{fruc}}=2.04\pm 0.11$, and ${\alpha}_{\text{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 bestfit model and the data occurs for two fructose steady states. These deviations would be corrected if $K}_{\mathrm{f}\mathrm{r}\mathrm{u}\mathrm{c}$ was larger, suggesting that the $K}_{\mathrm{f}\mathrm{r}\mathrm{u}\mathrm{c}$ used here may not reflect the actual value in the experiment. The model also accurately predicts the resulting steadystate 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 steadystate 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., 1993. 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 glucoselimited chemostat (we also compared our model to a batch seeded with E. coli from a galactoselimited chemostat, see Appendix 6—figure 2). For this data we used our serial dilution model with Monod kinetics and K_{i} values from measurements on the same strain (Lendenmann and Egli, 1998). As shown in Appendix 6—figure 1C, the agreement between the bestfit model and the experimental data is generally quite good over the entire time course. The estimated ${\alpha}_{i}$, measured in units of (mg sugar)(mg biomass)^{–1}h^{–1}, were 0.46 ± 0.04 and 0.41 ± 0.03 for glucose and galactose, respectively. The estimated yield was 0.42 ± 0.03, similar to the experimentally measured yield used in the chemostat model of 0.45 (mg biomass)(ms 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}_{\text{gal}}(t)=\text{const}\rho (t)Y{c}_{\text{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 (Kim et al., 2009), and has even been observed in nonprokaryotic organisms. For example, the eukaryote Kloeckera sp. 2201 has been shown to simultaneously utilize methanol and glucose as carbon sources (KovárováKovar and Egli, 1998). Similarly, the methanogenic archaeon Methanosarcina barkeri can simultaneously utilize methanol and acetate in batch culture (Scherer and Sahm, 1981). 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 (Scott et al., 2010), 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., 1996, we first digitally extracted the steadystate data points from the experimental figures. We used the model from Posfai et al., 2017 with Monod kinetics. The K_{i} of glucose, ribose, and fructose were taken as 73, 132, and 125 μg/L sugar, respectively (Lendenmann and Egli, 1998). 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(\text{g biomass}){(\text{g sugar})}^{1}$, as measured experimentally (Lendenmann et al., 1996). 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 δ 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 steadystate 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., 1993 and used the MATLAB curve fitting toolbox. The biomass was reported as ${\text{OD}}_{546}$ and was converted to mg/L using a conversion factor measured for the same strain (Lendenmann et al., 1996). Two sugar data points that were taken before the first biomass measurement were removed so that the initial conditions of the system would be welldefined. We estimated the parameters of the serial dilution model developed in this paper assuming Monod kinetics (Equation 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 μg/L, respectively (Lendenmann and Egli, 1998). 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 7
Supplemental figures
In this section we present supplemental figures that support the main text. Each figure’s caption contains all pertinent information.
Data availability
No datasets were generated in this work. All code and data used in this manuscript available at: https://github.com/AmirErez/SeasonalEcosystem (copy archived at https://github.com/elifesciencespublications/SeasonalEcosystem).
References

Analysis of bottlenecks in experimental models of infectionPLOS Pathogens 11:e1004823.https://doi.org/10.1371/journal.ppat.1004823

Rapid method of determining factors limiting bacterial growth in soilApplied and Environmental Microbiology 67:1830–1838.https://doi.org/10.1128/AEM.67.4.18301838.2001

Genetic revolutions, founder effects, and speciationAnnual Review of Ecology and Systematics 15:133–164.https://doi.org/10.1146/annurev.es.15.110184.001025

Multispecies competition in variable environmentsTheoretical Population Biology 45:227–276.https://doi.org/10.1006/tpbi.1994.1013

Mechanisms of maintenance of species diversityAnnual Review of Ecology and Systematics 31:343–366.https://doi.org/10.1146/annurev.ecolsys.31.1.343

Comparison of factors limiting bacterial growth in different soilsSoil Biology and Biochemistry 39:2485–2495.https://doi.org/10.1016/j.soilbio.2007.05.002

Kinetics of microbial growth with mixtures of carbon sourcesAntonie Van Leeuwenhoek 63:289–298.https://doi.org/10.1007/BF00871224

Diversity, stability, and reproducibility in stochastically assembled microbial ecosystemsPhysical Review Letters 120:158102.https://doi.org/10.1103/PhysRevLett.120.158102

Relaxed control of sugar utilization in Lactobacillus brevisMicrobiology 155:1351–1359.https://doi.org/10.1099/mic.0.0246530

Growth kinetics of suspended microbial cells: from singlesubstratecontrolled growth to mixedsubstrate kineticsMicrobiology and Molecular Biology Reviews 62:646–666.https://doi.org/10.1128/MMBR.62.3.646666.1998

Global marine bacterial diversity peaks at high latitudes in winterThe ISME Journal 7:1669–1677.https://doi.org/10.1038/ismej.2013.37

Kinetics of the simultaneous utilization of sugar mixtures by Escherichia coli in continuous cultureApplied and Environmental Microbiology 62:1493–1499.https://doi.org/10.1128/AEM.62.5.14931499.1996

Dynamics of mixed substrate growth of Escherichia coli in batch culture: the transition between simultaneous and sequential utilisation of carbon substratesActa Universitatis Carolinae Environmentalica 14:21–30.

Community equilibria and stability, and an extension of the competitive exclusion principleThe American Naturalist 104:413–423.https://doi.org/10.1086/282676

Concentrations and fluxes of organic carbon substrates in the aquatic environmentAntonie Van Leeuwenhoek 63:243–274.https://doi.org/10.1007/BF00871222

Variation in species richness: towards a unification of hypothesesFolia Geobotanica Et Phytotaxonomica 29:511–530.https://doi.org/10.1007/BF02883148

Metabolic TradeOffs promote diversity in a model ecosystemPhysical Review Letters 118:028103.https://doi.org/10.1103/PhysRevLett.118.028103

Influence of sulphurcontaining compounds on the growth of Methanosarcina barkeri in a defined mediumEuropean Journal of Applied Microbiology and Biotechnology 12:28–35.https://doi.org/10.1007/BF00508115

Microbial diversityproductivity relationships in aquatic ecosystemsFEMS Microbiology Ecology 62:181–186.https://doi.org/10.1111/j.15746941.2007.00381.x

Bacterial competition in serial transfer cultureMathematical Biosciences 229:149–159.https://doi.org/10.1016/j.mbs.2010.12.001

Partitioning of resources and the outcome of interspecific competition: a model and some general considerationsThe American Naturalist 107:171–198.https://doi.org/10.1086/282825

The relationship between productivity and species richnessAnnual Review of Ecology and Systematics 30:257–300.https://doi.org/10.1146/annurev.ecolsys.30.1.257

Continuous anaerobic digestion with Methanosarcina barkeriAnnals of the New York Academy of Sciences 413:489–500.https://doi.org/10.1111/j.17496632.1983.tb47929.x

Metabolic enzyme cost explains variable tradeoffs between microbial growth rate and yieldPLOS Computational Biology 14:e1006010.https://doi.org/10.1371/journal.pcbi.1006010
Decision letter

Detlef WeigelSenior Editor; Max Planck Institute for Developmental Biology, Germany

Alvaro SanchezReviewing Editor; Yale University, United States

Antonella SuccurroReviewer

Lei DaiReviewer; Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, China
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
This paper provides an elegant modeling framework for improving our understanding of how nutrients influence the diversity of microbial communities. As such, it will be a valuable addition to the field.
Decision letter after peer review:
Thank you for submitting your article "Nutrient levels and tradeoffs control diversity in a model seasonal ecosystem" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Detlef Weigel as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Antonella Succurro (Reviewer #1); Lei Dai (Reviewer #3).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
We would like to draw your attention to changes in our revision policy that we have made in response to COVID19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.
Summary:
Overall, the reviewers liked the work, and they all agreed that this is an insightful paper that should be suitable for publication if the authors address a few issues, none of which requires additional experimentation. Notably, all reviewers thought that the manuscript was a bit difficult to read at times, and they agreed that it would be substantially improved if the authors could shorten the paper and provide a more streamlined narrative, with a focus on its central results. Finally, the reviewers were skeptical of the comparison to experiments and indicated that it should be removed from the Abstract and reframed as suggested below.
The authors set out to understand the conditions in which n (n > 1) species can coexist on m < n resources. This is important because such coexistence on few resources appears common in nature. Yet, some simple competition models suggest it is impossible (hence the name of the conundrum as "the paradox of the plankton"). They seek to understand the conditions leading to stable diversity in a specific environmental mode: seasonal batchcultures with bottlenecks, a la serially culturing microbial communities, as they rightfully claim that theory here is less developed than in chemostat.
The most interesting finding in the paper is that nutrient abundance ("bolus") at the onset of each "season" is an important determinant of diversity, whereas reservoir concentrations are less important in a chemostat. The reason is what they call an "earlybird" effect: medium nutrient boluses tend to favor species whose consumption vectors are most similar to the bolus vector, allowing these species to rapidly outpace other species and then preempt those other species from resources. In contrast, with lower or much higher nutrient boluses, and especially, when the bolus size is small enough that the model is akin to a chemostat, the fixed energy budget among species causes them to equalize and allows for coexistence. A second finding is that there is a simple mechanism (using their pictorial simplexes) to predict diversity among a set of species, given an understanding of their nutrient consumption vectors and the vector of the nutrient supply. The nonmonotonic relation between nutrient bolus size and remapping of strategies, which determines species coexistence at steady state, is a novel and interesting result that can be tested by experiments.
Essential revisions:
1) All reviewers agreed that the two central and most insightful results are (i) the relationship between the bolus size and diversity and the associated remapping of metabolic strategies and (ii) the connection between serial passaging and chemostat modes. The paper would benefit by making it shorter and with a more clear focus on those two results. In particular the second result ((ii) above) could be highlighted in the Abstract / Introduction / Discussion. The authors clearly show that under certain conditions (notably, c_{0} << K), a chemostat model is a special case of a serial dilution model. This is a fascinating result. However, more effort could be taken to detail some of the important assumptions in this result. The main one is that batch culture is allowed to reach steadystate and full nutrient depletion. While this is reasonable from the goal of mathematically connecting the two environmental regimes, it is difficult to consider whether it is ecologically reasonable without some discussion of the amount of time required for that to be reached. For example, are there any batch culture simulations that would require an unreasonably long time to reach steadystate, such that it is unlikely to occur, and therefore unlikely to be equivalent to chemostat steadystate conditions?
2) Many of the other results besides the two highlighted above could be organized in the appendix, or as subsections of a final Results section, where the authors examine how relaxing some of the simplifying assumptions made in the model (e.g. the exploration of a softer tradeoff imposed by unequal enzyme budgets, the regulation of enzyme expression, the possibility of crossfeeding, the case of unequal K_{i} and Y_{i}, or the effect of including the stochastic nature of dilutions) affect their conclusions. If the authors prefer to include it as a Results section in the main text, they should introduce this section by explicitly alerting the reader that its purpose is to check the effect of relaxing the simplifying assumptions they previously made on the main results of the paper, listing them before going into details. That will help the reader understand the connection of the various subsections to the main results.
3) Reviewers appreciated the Mutation/Selection section, but they all agreed that it does not add much to the narrative. Cutting it off the paper would help to maintain the focus and streamline the narrative. The authors could develop it further in a separate paper, as it is an interesting and relevant topic
4) The comparison to experimental results should be revised. The analysis performed (only one species, not an ecosystem) in Figure 7 does not work as much as a "validation" of the main predictions of the paper, but rather as "support of the model assumptions". The authors described this as: "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." The reviewers concurred that validation experiments of the model in this paper would require testing the predictions of remapping: e.g. how coexisting species change as a function of nutrient bolus size. The paper is valid and interesting as it is and without an experimental validation of this prediction, so the current discussion on experiments should be either eliminated or moved to the appendix and changed in tone to reflect its value as a check of some of the model assumptions.
https://doi.org/10.7554/eLife.57790.sa1Author response
Summary:
Overall, the reviewers liked the work, and they all agreed that this is an insightful paper that should be suitable for publication if the authors address a few issues, none of which requires additional experimentation. Notably, all reviewers thought that the manuscript was a bit difficult to read at times, and they agreed that it would be substantially improved if the authors could shorten the paper and provide a more streamlined narrative, with a focus on its central results. Finally, the reviewers were skeptical of the comparison to experiments and indicated that it should be removed from the Abstract and reframed as suggested below. […]
We have followed the reviewers’ suggestions and have shorted and streamlined the narrative. As requested, we have reframed the comparison to experiments, moving it to the Appendix (details below).
Essential revisions:
1) All reviewers agreed that the two central and most insightful results are (i) the relationship between the bolus size and diversity and the associated remapping of metabolic strategies and (ii) the connection between serial passaging and chemostat modes. The paper would benefit by making it shorter and with a more clear focus on those two results. In particular the second result ((ii) above) could be highlighted in the Abstract / Introduction / Discussion. The authors clearly show that under certain conditions (notably, c_{0} << K), a chemostat model is a special case of a serial dilution model. This is a fascinating result. However, more effort could be taken to detail some of the important assumptions in this result. The main one is that batch culture is allowed to reach steadystate and full nutrient depletion. While this is reasonable from the goal of mathematically connecting the two environmental regimes, it is difficult to consider whether it is ecologically reasonable without some discussion of the amount of time required for that to be reached. For example, are there any batch culture simulations that would require an unreasonably long time to reach steadystate, such that it is unlikely to occur, and therefore unlikely to be equivalent to chemostat steadystate conditions?
We agree with the reviewers that the connection between serial passaging and chemostat models is a central result. As the reviewers suggest, we now dedicate a short section specifically to this connection. We now also highlight the connection to the chemostat in the Abstract and Discussion.
There are certain regimes in our model that do require a large number of batches to reach steadystate. This occurs near the transition between the highdiversity and lowdiversity regimes and is an example of the wellknown phenomenon of “critical slowing down”. We now explicitly acknowledge this critical slowing down in a paragraph in the new section on connections to the chemostat. We also address the question of finite timescales by probing the behavior of our model when nutrients are not completely depleted in each batch. We show in a new figure (Appendix 7—figure 2) that coexistence holds even when substantial portions of the nutrients are left unconsumed (in the figure, coexistence holds until about 45% of the nutrients are unconsumed).
2) Many of the other results besides the two highlighted above could be organized either in the appendix, or as subsections of a final Results section, where the authors examine how relaxing some of the simplifying assumptions made in the model (e.g. the exploration of a softer tradeoff imposed by unequal enzyme budgets, the regulation of enzyme expression, the possibility of crossfeeding, the case of unequal K_{i} and Y_{i}, or the effect of including the stochastic nature of dilutions) affect their conclusions. If the authors prefer to include it as a Results section in the main text, they should introduce this section by explicitly alerting the reader that its purpose is to check the effect of relaxing the simplifying assumptions they previously made on the main results of the paper, listing them before going into details. That will help the reader understand the connection of the various subsections to the main results.
We have acted on the reviewers’ suggestion. We have created a new Results section where we explicitly alert the readers that the purpose of that section is to relax some of the simplifying assumptions made earlier in the manuscript. We have moved the analyses of unequal K_{i} and Y_{i}, crossfeeding, population bottlenecks and unequal enzyme budgets to this new section.
3) Reviewers appreciated the Mutation/Selection section, but they all agreed that it does not add much to the narrative. Cutting it off the paper would help to maintain the focus and streamline the narrative. The authors could develop it further in a separate paper, as it is an interesting and relevant topic
We have followed the reviewers’ suggestion and have cut the mutation/selection section from the manuscript, and intend to publish it separately. We have also cut the “adaptation” section for the same reason, since we consider the “adaptation” section the least relevant for the rest of the manuscript.
4) The comparison to experimental results should be revised. The analysis performed (only one species, not an ecosystem) in Figure 7 does not work as much as a "validation" of the main predictions of the paper, but rather as "support of the model assumptions". The authors described this as: "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." The reviewers concurred that validation experiments of the model in this paper would require testing the predictions of remapping: e.g. how coexisting species change as a function of nutrient bolus size. The paper is valid and interesting as it is and without an experimental validation of this prediction, so the current discussion on experiments should be either eliminated or moved to the appendix and changed in tone to reflect its value as a check of some of the model assumptions.
We accept the reviewers’ criticism and have moved the current discussion of experiments to the appendix, changing its tone to reflect its value as a check of some of the model assumptions. As a result, there is a new Appendix 6, which includes all the relevant details of the experimental comparison. This appendix is referred to in the last paragraph of the Discussion.
https://doi.org/10.7554/eLife.57790.sa2Article and author information
Author details
Funding
National Institutes of Health (GM082938)
 Amir Erez
 Yigal Meir
 Ned S Wingreen
National Science Foundation (DGE1656466)
 Jaime G Lopez
National Science Foundation (PHY1734030)
 Amir Erez
 Jaime G Lopez
 Benjamin G Weiner
 Yigal Meir
 Ned S Wingreen
National Science Foundation (PHY1748958)
 Yigal Meir
 Ned S Wingreen
National Institutes of Health (R25GM067110)
 Yigal Meir
 Ned S Wingreen
Gordon and Betty Moore Foundation (2919.01)
 Yigal Meir
 Ned S Wingreen
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Senior Editor
 Detlef Weigel, Max Planck Institute for Developmental Biology, Germany
Reviewing Editor
 Alvaro Sanchez, Yale University, United States
Reviewers
 Antonella Succurro
 Lei Dai, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, China
Publication history
 Received: April 12, 2020
 Accepted: August 16, 2020
 Version of Record published: September 11, 2020 (version 1)
Copyright
© 2020, Erez et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,452
 Page views

 212
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.