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

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 resource-competition models don't allow for coexistence of a large number of species, it was recently shown that metabolic trade-offs can allow unlimited diversity. Does this diversity persist with more realistic, intermittent nutrient supply? Here, we demonstrate theoretically that in serial dilution culture, metabolic trade-offs allow for high diversity. When a small amount of nutrient is supplied to each batch, the serial dilution dynamics mimic a chemostat-like steady state. If more nutrient is supplied, community diversity shifts due to an 'early-bird' effect. The interplay of this effect with different environmental factors and diversity-supporting mechanisms leads to a variety of relationships between nutrient supply and diversity, suggesting that real ecosystems may not obey a universal nutrient-diversity relationship.


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 (Lloyd-Price et al., 2016;Ladau et al., 2013;Weigel and Pfister, 2019). These observations appear to clash with a prediction of some resource-competition models, known as the competitive-exclusion principlenamely, that steady-state 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 cross-feeding or antibiotic production and degradation (Goyal and Maslov, 2018;Kelsic et al., 2015), (ii) spatial heterogeneity (Murrell and Law, 2003;Tilman, 1994), (iii) persistent non-steady-state 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 self-organize 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 trade-offs, which assume that microbes have a limited enzyme production capacity they must apportion. Here, we investigate the equalizing effect of metabolic trade-offs in the context of serial dilution, to reflect a more realistic variable nutrient supply.
Serial dilution, in which cultures of bacteria are periodically diluted and supplied with fresh nutrients, is well-established as an experimental approach. For example, the bacterial populations in the Lenski long-term evolution experiment (Lenski and Travisano, 1994), experiments on community assembly (Goldford et al., 2018), and antibiotic cross-protection (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 trade-offs in their growth characteristics (Stewart and Levin, 1973;Smith, 2011), the theoretical understanding of diversity in serial dilution is much less developed than for chemostat-based steady-state growth.
Here, we show that under serial dilution, metabolic trade-offs can still support high diversity communities, but that this coexistence is now sensitive to environmental conditions. Interestingly, seasonality 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 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 competitive-exclusion 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 competitive-exclusion 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, 'early-bird' species, which rapidly use a subset of the available nutrients, grow to dominate the ecosystem. Even though 'late-bird' species are more effective at consuming the remaining resources, they cannot compete with the increased sheer numbers of the 'early-birds', leading to a 'rich-get-richer' phenomenon.
Oceanic plankton, arctic permafrost and many other threatened, resource-poor 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. diversity and the amount of nutrient provided to the community. These changes in diversity are driven by an 'early-bird' effect in which species that efficiently consume nutrients that are initially more abundant gain a population advantage early in the batch. To our knowledge, this is the first time this effect has been identified as a major influencer of diversity in seasonal ecosystems.
This dependence between community diversity and the supplied nutrient concentration allowed us to explore an unresolved question in ecology (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 'hump-shaped' unimodal trend predicted by many theories (Kassen et al., 2000). However, a meta-analysis 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 cross-feeding, relative enzyme budgets, differences in enzyme affinities, and differences in nutrient yields. That so much variation appears in a simple model suggests that real ecosystems are not likely to display a single universal relationship between nutrient supply and diversity.

Results
We employ the serial dilution model depicted in 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 fsg with initial biomass densities s ð0Þ in the batch such that the total initial biomass density is 0 ¼ P s s ð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; . . . ; p such that the total nutrient concentration is c 0 ¼ P p i¼1 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 Figure 1. Illustration of serial dilution resource-competition model. (A) Serial dilution protocol. Each cycle of batch growth begins with a cellular biomass density 0 and total nutrient concentration c 0 . The system evolves according to Equations 2-3 until nutrients are completely consumed. A sample of the total biomass is then used to inoculate the next batch again at density 0 . (B) Representation of particular enzyme-allocation strategies fa s g (colored circles) and nutrient supply composition c i =c 0 (black diamond) on a 2-nutrient simplex, where the right endpoint corresponds to c 1 =c 0 ¼ 1. (C) Representation of particular strategies (circles) and nutrient supply (black diamond) on a 3-nutrient simplex. Dashed blue -the convex hull of the enzyme-allocation strategies. Here, the nutrient supply (black diamond) is inside the convex hull, implying coexistence of all species in the chemostat limit (see text). 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 s ðtÞ increase with time, starting at t ¼ 0, and growth continues until the nutrients are fully depleted, P p i¼1 c i ð¥Þ » 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 P s s ð¥Þ ¼ 0 þ c 0 . The next batch is then inoculated with a biomass density 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 s is defined by its unique enzyme strategyã s ¼ ða s;i ; . . . ; a s;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 s consumes nutrient i at a rate j s;i (per unit biomass) that depends on nutrient availability c i and on its enzyme-allocation strategy a s;i according to For simplicity, we take all Monod constants to be equal, K i ¼ K (a more general form of the Figure 2. The nutrient bolus size c 0 affects the relative abundance of species and even their coexistence at steady state. (A) Schematic of the mutual invasibility condition for two species and two nutrients. Top: The red species can be invaded by any species with a strategy to its left if the supply lies in the region marked by the hatched rectangle. Middle: Similarly, showing the supplies for which blue can be invaded by any species with a strategy to its right. Bottom: The intersection defines a mutual invasibility region of supplies for which the two species red and blue will coexist. Triangles mark the boundaries of this coexistence region. (B-D) Example of the effect of c 0 on coexistence for more than two species: the approach to steady state, showing s versus batch number (left column) with the corresponding c 0 -dependent remapping of coexistence boundaries (right column). (B) For the chemostat limit c 0 ( K, where K is the Monod constant for nutrient uptake, the triangles marking coexistence boundaries coincide with the species' strategies, a s . (C) For c 0 » K the triangles are remapped towards the center of the simplex compared to the strategies fa s g. In this example the nutrient supply (black diamond) ends up outside the coexistence boundaries, so only one species survives. (D) For c 0 ) K the triangles again coincide with the strategies fa s g, leading again to coexistence. 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 s;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 P i a s;i ¼ E þ " s , where s is a zero-mean and unit-variance Gaussian variable. Without loss of generality we take E ¼ 1; initially we set " ¼ 0, which we call exact trade-offs. This allows us to visualize the strategiesã s 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 ( 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 trade-offs, 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 steady-state levels. Such a self-organized neutral state is an attractor of the chemostat dynamics (Posfai et al., 2017) and, correspondingly, of the c 0 ( 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 chemostat-limit 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 convex-hull rule, note that a state of arbitrarily high coexistence can only occur if the chemostat self-organizes to a 'neutral' state in which the nutrient concentrations are all equal, and thus all strategies have the same growth rate. This state is achieved if and only if the total enzyme abundances lie along the same vector as the nutrient supply, which is achievable only if the supply lies within the convex hull of the strategies present.
As in the chemostat model, the serial dilution model can support either coexistence or competitive exclusion. However, if one chooses system parameters near the transition between these two states, it requires a very large number of batches for the simulation to reach steady state. This is a manifestation of the well-known phenomenon of critical slowing down (Dakos and Bascompte, 2014). Though in principle, critical slowing down is not a simulation artifact and could manifest in similar real-world systems, we expect that a variety of factors outside our modeling framework would preclude observation of this critical behavior.
We define the serial dilutions 'steady state' to be reached when the relative species abundances after the nutrients are depleted (time t f after starting the batch) scale with the relative abundances at the beginning of that batch (time 0), that is, s ðt f Þ ¼ 0þc0 0 s ð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,d we reach the c 0 =K ( 1 steady-state condition for the serial dilution system: : Averaged over a batch, s i is the average rate that nutrient i is supplied, andd 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 parameterd which would be the rate all nutrients are continuously supplied per unit biomass, would need to equal the dilution rate of the chemostat, d, 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 steady-state outcome of serial dilutions. To understand this effect, we first consider a simple case of two nutrients and two species as depicted in 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 species-specific boundaries remapped because they generally lie at different locations on the simplex than the strategies they originated from, with the extent of remapping depending on the nutrient bolus size. At a more technical level, the remapped boundary for a given species and bolus size is the nutrient supply for which, over the course of a batch, all nutrients are equally valuable and so, a species with any strategy can neutrally invade and persist. This equality of nutrient value is defined in terms of the Monod function integrals for each nutrient (for details see Appendix 3).
Since the remapped coexistence boundaries depend on the nutrient bolus size c 0 , changing bolus size can qualitatively change the steady-state outcome of serial dilutions. Figure 2B-D depicts an example of how c 0 affects remapping, and the consequences for species coexistence. At low bolus size, c 0 ( K, corresponding to the chemostat limit, Figure 2B (left) shows that all species present achieve steady-state coexistence. This follows because the nutrient supply (black diamond) lies inside the convex hull. When c 0 is increased to c 0 » 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 predator-prey systems (Rosenzweig, 1971). Strikingly, however, as bolus size is further increased to c 0 ) K, the coexistence boundaries are remapped back to their original positions, so that the nutrient supply once again lies within the convex hull, and so steady-state coexistence of all species is recovered.
What causes the remapping of the coexistence boundaries inwards as c 0 =K ! 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 a s;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 'early-bird' advantage: its initial growth allows it to more effectively exhaust the nutrients.
Why does the coexistence boundary of a species map back to its original strategy in the limit of large bolus size, c 0 ) 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 a s;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 (" ¼ 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 s is P i a s;i c i =ðK i þ c i Þ, so that the fold increase of a species during a batch is expðã s ÁĨÞ. 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 growth-dilution batch procedure results in a steady state where the populations at the beginning of a batch do not change from batch to batch. The steady-state populations depend on the initial populations, with the set of all possible steady-state populations defining a coexistence manifold.

Steady-state diversity
As is apparent in 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 » K; (iii) there is no remapping in both the c 0 ! 0 and c 0 ! ¥ limits (see also Appendix 7-figure 4). The extent of remapping also depends on the inoculum size 0 as shown in Figure 3B, which demonstrates that remapping is maximal for 0 ( K and vanishes for 0 ) K.
How does nutrient bolus size influence steady-state 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 ¼ À P s P s ln P s and P s ¼ Ã s ð0Þ= 0 , with Ã s ð0Þ the steady-state 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 steady-state 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 chemostat-limit of small bolus sizes c 0 ( K and drops even further for c 0 » K. This loss of diversity is due to the dramatically lowered steady-state abundances of strategies that favor Nutrient 1 ( Figure 3D, bottom). Two different effects underlie this change in community structure. The first is the early-bird effect described above: species specializing in more abundant nutrients gain a population advantage that allows them to rapidly consume less abundant nutrients that would otherwise support species with different enzyme specializations. The second effect is a well-known property of single nutrient competition and can be viewed as a 'single-nutrient' early-bird effect. In this case, species that are superior competitors for a nutrient gain an exponential population advantage over inferior competitors, increasing their share of total nutrient beyond the ratio of initial consumption rates. Both of these effects increase in strength with larger bolus size because the early-bird advantage increases as growth proceeds. The combination of these effects results in the species specialized in consuming the most abundant nutrients consuming a larger fraction of all nutrients. However, for very large bolus sizes, saturation of nutrient uptake rates mitigates these two effects, leading to a lack of remapping for c 0 ) K and diversity returning to its chemostat-limit. Though here we focused on the case of two nutrients, these results extend to more nutrients (for three nutrients see Appendix 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 6 ¼ K, and different nutrient yields Y i 6 ¼ Y. This is followed by a model that allows cross-feeding of metabolites. Finally, we consider population bottlenecks and what happens when the fixed-enzyme-budget constraint is relaxed. We show that in  Figure 3C, but with K 1 ¼ 10 À3 and K 2 ¼ 0 ¼ 1. Colors correspond to different nutrient supply compositions, solid curves for c 1 =c 0 2 ½0; 0:5 and dashed curves for c 1 =c 0 2 ½0:5; 1. Dashed black line: maximum diversity (equal species abundances) is no longer attained when nutrient composition is ð0:5; 0:5Þ. (B) Steady-state species abundances f s g for nutrient composition ð0:5; 0:5Þ (top) and ð0:05; 0:95Þ (bottom), as in Figure 3D. (C) Steadystate effective number of species m e as a function of bolus size c 0 =K, as in Figure 3C, but with Y 1 ¼ 10 and Y 2 ¼ 0 ¼ 1. Note that the nutrient compositions are normalized to yield such that ð0:5 Ã ; 0:5 Ã Þ is actually ð0:5=Y 1 ; 0:5=Y 2 Þ. Colors the same as in A. (D) Steady-state species abundances f s g for nutrient composition ð0:5 Ã ; 0:5 Ã Þ (top) and ð0:2 Ã ; 0:8 Ã Þ (bottom). all these cases, the dependence on bolus size can be understood as manifestations of the early-bird effect.

Unequal enzyme affinities and nutrient yields
We have thus far made the simplifying assumption that all enzymes have the same substrate affinity, such that K i 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 steady-state behavior unchanged. However, outside the chemostat regime, differences in the K i have a drastic effect: when the nutrient with the lower K i is scant in supply, diversity generally increases with increasing c 0 , while the opposite occurs when the nutrient with the lower K i is higher in supply.
We can understand these K i -driven shifts in the nutrient-diversity relationship as due to changes in the identity of the early bird. In a model with equal K i , the identity of the early bird is determined by which nutrient is more abundant: if the two nutrients have equal K i , a species can gain an 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 early-bird advantage. Why does this change in the early bird's identity change the form of the nutrient-diversity relationship? This change arises from a clash between optimal feeding behavior in the chemostat and seasonal regimes. In the chemostat, it is advantageous to focus on the most abundant nutrient, regardless of the value of K i . Thus, in the chemostat limit, species focusing on the more abundant nutrient have an advantage. In the case of equal K i (or K i favoring the more abundant nutrient), this advantage is intensified by the early-bird effect, increasing the biomass of already abundant species and lowering diversity. By contrast, if the low abundance nutrient has a low K i , the early-bird effect will have the opposite effect on diversity. Now the early-bird effect benefits species that were disadvantaged in the chemostat limit, leading to more equal abundances and higher diversity. This shift in abundances is shown in 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 symmetry-related bolus-composition cases (e.g. [0.2,0.8] and [0.8,0.2]) eventually diverge, with one's diversity rising and the other's continuing to fall. This behavior is explainable by the same logic as in the variable K i case: diversity rises or falls depending on whether the early-bird species was also favored in the chemostat limit. However, unlike the case of variable K i , the diversity curves do not eventually return to the chemostat limit. Regardless of which nutrient the Y i favor, the diversity eventually begins decreasing monotonically as c 0 =K increases. This difference between the variable K i and variable Y i cases can be understood by considering what occurs when both nutrients are saturating. In the variable K i case, saturating nutrients are equal in value, implying a return to the chemostat limit as the early-bird effect weakens. In contrast, for variable Y i , there remains a difference in the value of the two nutrients in the saturated regime, meaning that the early-bird effect will grow stronger and the early bird will take over the population. The beginning of this takeover can be seen at high bolus sizes in 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 early-bird 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.

Cross-feeding
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 cross-feeding, which has generally been found to promote diversity (Goyal and Maslov, 2018) and stable community structure (Goldford et al., 2018). Here, cross-feeding is introduced through the byproduct matrix G s i;i 0 , which converts the consumption of nutrient i 0 to production of nutrient i such that, In this framework, nutrient i 0 is converted to nutrient i at no extra enzymatic cost, meaning that nutrient i is simply a byproduct whenever nutrient i 0 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 G s i;i 0 ¼ 0 G 0 0 for all species. When G ¼ 1, upon consumption Nutrient 2 is perfectly converted to Nutrient 1, leading to an equal total supply of the two nutrients. More generally, R ¥ 0 P s s j s;1 dt ¼ Gc 2 ð0Þ which allows a direct comparison between the unitrophic and bitrophic regimes: starting with c 2 ð0Þ results in ðG þ 1Þc 2 ð0Þ total nutrient, and hence the Nutrient 1 fraction is G 1þG of the total. How does cross-feeding influence diversity in our serial dilution model? In Figure 5A we compare bitrophic diversity for six values of G to their unitrophic equivalents (in Figure 3C). We note that: (i) bitrophy still supports diversity greater than the competitive-exclusion limit; (ii) in the chemostat regime, c 0 ( K, the unitrophic and bitrophic schemes have identical values of m e , and these drop as c 0 ! K; (iii) but for bitrophy the m e does not recover for c 0 ) K; (iv) even when the total supply of both nutrients is equal (G ¼ 1), bitrophy leads to lower than maximal m e outside the chemostat limit. These features are clarified in Figure 5B, which shows steady-state species abundances for G values leading to a total Nutrient 1 supply fraction of 0.5 and 0.05, and highlights the lower diversity for bitrophy compared to unitrophy for large nutrient bolus size. This difference is due to an early-bird effect: the species consuming supplied nutrient early in the batch can build a sizable population before the competing species that rely on its byproduct. The early-bird population then outcompetes the others for byproduct consumption. As such, this effect increases with c 0 = 0 . The effect also becomes stronger at low c 0 =K (with constant c 0 = 0 ), since this allows the early-bird 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 single-nutrient specialists, then changes in c 0 =K have no impact on community diversity.
The behavior of the model with cross-feeding shows that the early-bird effect extends beyond simple metabolic trade-offs. More broadly, when species compete for multiple resources that are supplied in batches, a species' survival depends on more than its ability to efficiently consume nutrients. An early-bird species, being more specialized in consuming the nutrients that are initially more abundant, gains a population advantage early in the batch. This population advantage may allow the early-bird species to out-compete other species even when consuming nutrients it is not specialized to consume. Despite its consumption inefficiencies, through sheer numbers the 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 semi-isolated 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 near-complete loss of diversity. For small real-world populations, however, diversity may be maintained by migration. To model such migration we augmented the population at each dilution with a 'spike-in' from a global pool of species, in the spirit of MacArthur's theory of island biogeography (MacArthur and Wilson, 2001). Specifically, in the spike-in procedure, to prevent extinctions caused by sampling fluctuations, every new batch is inoculated with a small number of the original, founder species.
In Figure 6A we show results of spike-in serial dilutions for a population bottleneck of 1008 cells. 95% of these cells are sampled from the previous batch, while 5% are sampled from a global pool, with equal abundances of 21 equally spaced strategies (cf. Figure 3A). The resulting m e vs. c 0 curves have maximal m e for all six nutrient fractions in the regime c 0 ( K where the 5% spike-in dominates sampling noise. As expected, for a balanced nutrient supply at any c 0 , all species have the same average abundance ( 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 spike-in protocol leads to higher diversity at low c 0 than the deterministic case (starting from equal species abundances but with no spike-in, 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 spike-in maintains a higher level of the least competitive strains, but since these abundances are still low, this difference in not reflected in the m e values.

Unequal enzyme budgets
While we have assumed exact trade-offs to achieve diversity within a resource-competition model, the trade-offs present among real microorganisms will not be exact. For the serial dilution protocol with spike-ins, diversity is maintained by migration and so it is possible to relax the constraint of exact trade-offs. How does diversity depend on the nutrient supply if we allow species to have different enzyme budgets? We implemented random differences in species enzyme budgets by setting " ¼ 0:1, that is, a standard deviation of 10%, and plotted effective number of species m e in Figure 6C. As in the " ¼ 0 limit ( Figure 6A), at sufficiently small c 0 the spike-in procedure dominates both sampling noise and differential growth rates due to unequal enzyme budgets. Raising c 0 leads to a drop in m e (albeit still above the competitive-exclusion limit). Examining the species abundances in 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 ) 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 " ¼ 0 case. However, increasing c 0 results in increased abundance for the species with the highest enzyme budget -which would ultimately lead to its domination for sufficiently large c 0 . This increasing dominance of the species with the highest enzyme budget is another manifestation of the early-bird effect: as the amount of growth in a batch increases, the advantage of a larger enzyme budget further compounds. In short, for spike-in serial dilutions the influence of unequal enzyme budgets depends on the nutrient supply, such that the species with the largest budgets dominate for large, unbiased supplies.

Discussion
Natural ecosystems experience variations in the timing and magnitude of nutrient supply, and the impact of these variations on species diversity is not fully understood (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 'early-bird' effect associated with supplying nutrients as large seasonal boluses instead of continuously. Some species can capitalize on rapid initial growth on an abundant nutrient to reach a large population size, which then allows them to deplete the remaining nutrients at the expense of their competitors. This 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 trade-offs are combined with stabilizing mechanisms, the impact of the early-bird effect remains. For example, in the case of crossfeeding, the community diversity falls as a function of c 0 =K due to the early-bird advantages gained by species at higher trophic levels.
While the idea of species gaining early advantages has been explored, such as in the literature on founder effects and speciation (Barton and Charlesworth, 1984;Brown,, 1957), to the best of our knowledge this is the first demonstration of the influence of the early-bird effect on the diversity of seasonal ecosystems. We believe that this effect will occur in a variety of such ecosystems, as its only fundamental requirement is competition for multiple nutrients that are supplied in a time-dependent manner. Interestingly, while the early-bird effect plays a large role in our model, it is not the only bolus-dependent effect that influences diversity. We also observe another effect that can be viewed as a 'single-nutrient' version of the early-bird effect. This effect arises from a well-studied property of competition: as growth proceeds, a superior competitor for a nutrient gains an exponential advantage over inferior competitors for that nutrient. Like the early-bird effect, this shifts the system's biomass towards species more specialized in initially abundant nutrients, particularly for large but non-saturating nutrient bolus sizes. This single-nutrient effect can co-occur with the early-bird effect, for example in competition for an abundant nutrient between two early-bird species.
The form of seasonality we explore in this manuscript, where mixed boluses are supplied periodically, is only one possible form of seasonal nutrient supply. The impact of the early-bird effect and single-nutrient competition will likely differ between different forms of seasonality. For example, we show in Appendix 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 single-nutrient competition effect initially observed in our mixed-bolus models. We expect the principles gleaned from our models to be of use in understanding diversity in a variety of seasonal ecosystems.
Finding a general relation between the amount of nutrient supplied to a community and its diversity is a long-standing goal of theoretical ecology (Tilman, 1982;Abrams, 1995;Leibold, 1996). We found that in our model the form of the nutrient-diversity relation (NDR) can change based on model details. The model has two regimes: a low diversity and a high diversity regime. The former satisfies competitive exclusion (no more species coexisting than resources), whereas the latter exceeds competitive exclusion and occurs when the nutrient supply lies within the convex hull of the remapped metabolic strategies present (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 early-bird effect and other mechanisms. Even with a single trophic layer, the NDR can be U-shaped, hump-shaped, monotonically decreasing, or have multiple peaks. These trends can then be further modified by the addition of more trophic layers, differences in enzyme budgets, etc.
Experimental studies that characterize the NDRs of microbial ecosystems have reached similarly variable conclusions. For example, one work studying bacterial communities in Arctic deep-sea sediments found an increasing trend between energy input and richness (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 meta-analysis of aquatic microbial ecosystems found examples of both monotonic and non-monotonic 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 one-size-fits-all theory, but to analyze the role of different factors such as cross-feeding, trade-offs, and immigration in determining that particular ecosystem's NDR. While we have focused on microbial systems, the absence of a universal NDR is consistent with results from recent work in plants (Adler et al., 2011).
We found that the stringency of metabolic trade-offs has a large impact on community diversity. We imposed a metabolic enzyme budget on each species to reflect the reality that microbial cells have a finite capacity to synthesize proteins and must carefully apportion their proteome (Basan et al., 2015). However, while it is true that microbes have limited biosynthetic capacity, it is unclear how strict are the resulting trade-offs. For this reason, we characterized versions of the model with both exact and inexact trade-offs. Our results show that the form of an ecosystem's NDR can depend on the stringency of metabolic trade-offs. This finding is not exclusive to the serial dilution model. The stringency of trade-offs was also important in the original chemostat setting: in a birth-death-immigration framework, small violations of the enzyme budget still allowed for high levels of coexistence, but large violations disrupted coexistence (Posfai et al., 2017). These results suggest that an experimental characterization of the stringency of metabolic trade-offs among microbes would provide a valuable ecological parameter. Note that metabolic trade-offs are only one of the many types of trade-offs microbes are subject to; other types of trade-offs, such as constraints between biomass yield and growth rate (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 8i; 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 ) 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 trade-offs, since some species will have a strict advantage over others. Thus, it is likely that these unequal yields between species will lead to a reduction in community diversity. However, varying the yield in this manner also allows for the inclusion of new trade-offs that may impact diversity, such as the aforementioned trade-off between yield and growth rate (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 low-abundance nutrient also has a low K i , the early-bird effect favors species that were disadvantaged in the chemostat limit, thus reversing the equal-K i NDR and leading to humpshaped NDR curves. Indeed, large differences in K i values can lead to a multi-peaked 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 early-bird effect, are in principle testable. To connect our modeling assumptions to real microbial systems, we compare our growth model of substitutable and simultaneous nutrient consumption to previously published experimental data from Escherichia coli growing in batch and chemostat conditions. We find that our modeling assumptions are consistent with both datasets and outline potential future experiments to test the model's multispecies predictions, detailed in Appendix 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, nutrient-diversity relationships could then be obtained by competing strains in serial dilution culture and measuring the community diversity (e.g. via fluorescent tags or by 16S rRNA sequencing) as a function of the total concentration of multiple, substitutable nutrients provided at the start of each batch. 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. Additional files

Supplementary files
. Transparent reporting form

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).
The following datasets were generated: Byproduct matrix converting nutrient i 0 to nutrient i j s;i Nutrient i consumption rate by species s

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 wherẽ a 1 is to the left ofã 2 on the 1-simplex. Species 2 can invade Species 1 if the nutrient supply is to the right ofã 1 . Species 1 can invade Species 2 if the nutrient supply is to the left ofã 2 . Therefore, the two species can mutually invade and coexist if and only if the nutrient supply lies betweenã 1 andã 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 P i a s;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.
Perturbation theory for c 0 =K ( 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 0 =K, respectively. Alternatively, one might choose to expand around a third ratio, c 0 = 0 , which would be useful to model extremely nutrient-dilute conditions as found in some marine microbial ecosystems.
We define a perturbation expansion with respect to the small parameter f ¼ c 0 =K, We note that at Oð1Þ we have s ðtÞ ¼ s ð0Þ and c i ðtÞ ¼ 0 as expected. We begin by expanding the Monod function, Accordingly, in the kinetic equation for c i , We next solve for ð1Þ s using c i , and then we will use ð1Þ s to solve for c ð2Þ i . It is possible but not necessary for our purposes to iterate further. The kinetic equation for the biomass density s is, Substituting the f expansion gives, to leading order, Taking the long-time limit, t ) 1 g i , we obtain, Focusing on the leading order, we conclude that, Substituting for g i and for c : Explicitly stating the batch number d, at the end of the batch, that is, at time t ¼ t f ) g À1 i , the biomass density is, In the serial dilution model with complete consumption of all nutrients c 0 and initial biomass 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 s ðd þ 1; 0Þ ¼ s ðd; 0Þ: Our calculation, to order f, gives, : Dividing both sides by t f and defining,d we finally reach the c 0 =K ( 1 steady-state condition for the serial dilution system: : Averaged over a batch, s i is the average rate that nutrient i is supplied, andd 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 parameterd, which would be the rate all nutrients are continuously supplied per unit biomass, would need to equal d, the dilution rate of the chemostat to maintain steady state. Indeed, Equation 32 is precisely the steady-state condition for the chemostat (Equation 4 from Posfai et al., 2017) with s i andd ¼ d interpreted as above.
Thus we complete the proof that in c 0 =K ( 1, the steady state of our serial dilution model is identical to the steady state of the equivalent chemostat model. Second-order corrections to remapping of the coexistence boundaries for c 0 =K ( 1 We have demonstrated above that the leading terms in an expansion for small nutrient supply retrieve the steady-state solution of the chemostat model. However, we know from numerical simulations, that as c 0 =K f 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 f as proved above, to capture the remapping we expand to order f 2 in c i ðtÞ. To this end, we return to the f expansion and extract the f 2 contribution, Kg j te Àg i t À e Àg i t À e Àðg i þg j Þt g j ! : Now we can solve for the growth-function integrals I i for the case of a single species growing in isolation. We expand the growth-function integral, i þ ::: Substituting the order f from Equation 22 and integrating, gives: For a single species growing in isolation, g i ¼ 1 K a s;i s ð0Þ. Thus, to order f, the chemostat limit, we obtain, (in the chemostat limit, for a single species), Thus, to satisfy the coexistence boundary conditions: 8i : I i ¼ const, to order f it must be that 8i : c i ð0Þ=a s;i ¼ 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 f 2 .
To order f 2 we substitute Equation 22 for c ; (38) which upon substituting g i for a single species simplifies to: Collecting terms to order f 2 gives As before, the coexistence boundaries are defined by I i ¼ const. To order f 2 , Equation 40 can be used to solve for this remapping analytically. Equation 40 also clarifies why perfect generalists (8i; j : a s;i ¼ a s;j ) do not get remapped, as stated in the main text. This is because for a generalist 8i; j : a 2 s;j =ða s;i þ a s;j Þ ¼ const., meaning that the second-order I i will all be equal if the first-order I i are all the same. Moreover, the term I ð1Þ j / À1 s ð0Þ multiplies the order f 2 correction to the coexistence boundary, and as a result, the larger s ð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 0 =K increases because at high 0 =K the remapping is small.
Second-order corrections to the steady-state abundance manifold for c 0 =K ( 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 d and nutrient supply rate s i (Posfai et al., 2017), We use the asterisk in Ã s 0 to make explicit that the abundances are the steady-state abundances, after the system has moved from its initial conditions. As a result, the steady state constrains the abundances Ã s , such that they must lie on a manifold of solutions that satisfy P s a s;i Ã s ¼ E d 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 identifyingd ¼ c0 0tf and s i ¼ cið0Þ tf , to first order in f c 0 =K ( 1, the steady-state condition for the serial dilution system is identical to the steady state of the equivalent chemostat. Thus, to order f, In the serial-dilution framework, to leading order in f, the chemostat limit of the serial-dilutions steady state is, We note that Kg i from Equation 43 is a solution of Equation 42, with P i a s;i ¼ E. Having established Equation 43 as the leading order (chemostat limit) term in an expansion in f, we proceed to calculate Kg i to order f 2 to obtain corrections to the chemostat limit.
From Equation 3, we obtain s ðtÞ ¼ 0 ð0Þe P i as;iIi , 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 steady-state species abundances Ã s ð0Þ obey Indeed, Ã s ð0Þ cancels, yielding the serial-dilutions steady-state to order f 2 , X j a s;j c j ð0Þ=c i ð0Þ c i ð0Þ þ c j ð0Þ : Comparing Equation 47 with Equation 43 we note the small, order c 0 = 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, 8i : c i ð0Þ ¼ c 0 =p, we have X i ¼ 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. 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 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 well-defined. 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 mg/L, respectively (Lendenmann and Egli, 1998). The three fitting parameters were the yield Y and the strategies a 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. own abundance, c i =ðK þ c i Þ » c i =K, and there is little relative change in biomass, ðtÞ » 0 . The more abundant nutrient is consumed faster (since the strategy is matched to nutrient proportions) and the majority of it is consumed quickly. The less abundant nutrient is consumed more slowly and a significant portion of it remains after the more abundant nutrient is almost completely depleted. In this way, the growth timecourse integrals are balanced to be equal. Once c 0 increases relative to 0 , but c 0 is not large compared to K this balance is broken. The more abundant nutrient will still be depleted quickly. However, now that c 0 is larger this initial consumption results in an increased abundance of the consumer. This means that the less abundant nutrient is depleted more quickly, leading a smaller growth timecourse integral for this nutrient relative to the more abundant nutrient. Thus, in this regime, the difference between the growth timecourse integrals increases with c 0 . Restoring them to equality requires more equal starting distributions of nutrients. Once c 0 increases such that c 0 ) K and c 0 ) 0 , the increase in c 0 now drives the growth timecourse integrals towards equality. The growth function is now almost always saturated and this neutralizes the effect of one nutrient starting at a much larger concentration. There will now be a significant buildup of biomass before the nutrients are exhausted, meaning that the growth timecourses will subsequently drop very quickly. As the growth function becomes more and more saturated, the nutrients will be consumed in proportion to their strategy. Thus, the growth timecourse integrals will once again be equal since the strategies match the nutrient proportions and the 'crash' times will therefore be similar for both nutrients. Appendix 7-figure 4. Dependence of coexistence boundary remapping on c 0 =K. As a further exposition to Figure 3A in the main text, shown here is the difference between the remapped coexistence boundaries and the corresponding metabolic strategies as a function of c 0 =K and metabolic strategy with 0 =K ¼ 10 À3 .
Appendix 7-figure 5. Serial dilution model with three nutrients. (A) Example of remapping on the three-nutrient simplex, similar to Figure 2 from the main text. Here we show how the remapping analysis presented in the main text for two nutrients can be extended to three nutrients. Remapping of three strategies for c 0 =K ¼ 1 and 0 ¼ 10 À2 . Outer circles: strategies fã s g; inner triangles: remapped nodes; lines connecting outer circles: supplies within this convex hull of strategies lead to coexistence of all species in the chemostat regime c 0 =K ( 1; dashes connecting inner circles: approximate remapped convex hull boundary defining region of supplies leading to coexistence for c 0 =K ¼ 1. Note that, as in the two-nutrient case, the strategies map inwards on the simplex for c 0 =K » 1. (B) Steady-state effective number of species m e versus c 0 =K for equal initial inocula of 64 species equally spaced throughout the triangular simplex competing for three nutrients. Effective number of species shows the same trend of loss in diversity when c 0 =K » 1 as in the two-nutrient case in Figure 2C.
that varying the values of K i and Y i can influence the direction and magnitude of the remapping, using a s1 ¼ 0:2 with 0 ¼ 1 as an example. (A) Remapping with different combinations of K i . The strategy we are examining devotes most of its enzyme budget to consuming nutrient 2. When there is a large different in K i that favors nutrient 2, inward remapping is enhanced. When there is a large difference that favors nutrient 1, outward remapping is enhanced. Eventually, the remapping returns to the chemostat limit. (B): Remapping with different combinations of Y i . Note that when yields are variable the condition for the remapped point is The remapped points shown are normalized to yield (if the remapped point is ðx; 1 À xÞ the normalized form is ðx Ã ; 1 À x Ã Þ where x Ã ¼ Y1x Y1xþð1ÀxÞY2 ). Similar to the unequal K i case, inward remapping is enhanced when there are large differences in Y i favoring nutrient 2. When the Y i favor nutrient 1, outward remapping is enhanced. Unlike in the unequal Y i case, the remapping does not return to the chemostat limit.
Appendix 7-figure 8. Batch timecourses of a bitrophic model with only two species. To further investigate the difference between the unitrophic and bitrophic scenarios, we consider a toy system with only two species, Species 1 with strategy ð0:05; 0:95Þ and Species 2 with strategy ð0:95; 0:05Þ. We set the byproduct matrix for perfect conversion, G 1;2 ¼ 1 and the nutrient bolus composition so that only Nutrient 2 is provided. By the end of each batch, the same amounts of Nutrient 1 and Nutrient 2 have been consumed. (A) Simulations with constant c 0 = 0 and variable c 0 =K. The 'earlybird' effect becomes stronger with decreasing c 0 =K, since this allows the dominant species to grow more before the byproduct can be readily consumed. (B) Simulations with variable c 0 = 0 and constant c 0 =K. The 'early-bird' effect is stronger at higher c 0 = 0 because the larger amount of supplied nutrient allows the dominant species to build a large population which can then outcompete other species. Conversely, if instead of cross-feeding we were to supply in the nutrient bolus equal quantities of Nutrients 1 and 2, the result would be equal abundance of both species.