Dynamics of bacterial populations under the feast-famine cycles

Bacterial populations in natural conditions are expected to experience stochastic environmental fluctuations, and in addition, environments are affected by bacterial activities since they consume substrates and excrete various chemicals. We here study possible outcomes of population dynamics and evolution under the repeated cycle of substrate-rich conditions and starvation, called the"feast-famine cycle", by a simple stochastic model with the trade-off relationship between the growth rate and the growth yield or the death rate. In the model, the feast (substrate-rich) period is led by a stochastic substrate addition event, while the famine (starvation) period is evoked because bacteria use the supplied substrate. Under the repeated feast-famine cycle, the bacterial population tends to increase the growth rate, even though that tends to decrease the total population size due to the trade-off. Analysis of the model shows that the ratio between the growth rate and the death rate becomes the effective fitness of the population. Hence, the functional form of the trade-off between the growth and death rate determines if the bacterial population eventually goes extinct as an evolutionary consequence. We then show that the increase of the added substrate in the feast period can drive the extinction faster. Overall, the model sheds light on non-trivial possible outcomes under repeated feast-famine cycles.


INTRODUCTION
As already pointed out in the 18th century, exponential growth is the most prominent feature of population dynamics [1], and bacterial systems are probably the best-studied model system about exponential growth. However, as pointed out by J.
Monod [2], the exponential growth is only one of the growth phases of bacteria, and the stationary phase, the death phase, and the lag phase are all as important for the bacterial population dynamics.
A variety of theoretical studies on the bacterial population dynamics tend to focus on the competition for nutrients under constant environment [3,4], where the competition takes place mainly in the form of exponential growth under a constant influx of substrate (combined with dilution/death to keep the environment constant).
While these models have provided fundamental insights into bacterial population dynamics, in natural environments like ponds, soils, and puddles, the nutrients may be supplied by rarely happening events rather than continuous influx. Under such natural environments, bacteria experience substrate rich conditions and poor conditions alternately.
This cycle between substrate rich and poor conditions is called the feast-famine cycle [5][6][7][8][9][10][11]. In contrast to the continuous nutrient supply (or the constant environment) condition, under the feast-famine cycle there is no steady-state for the amount of the substrate, and accordingly for the number of the cells. While the environment becomes substrate-rich for some time after the substrate addition event, once the cells in the environment run out all the substrates, they have to tolerate until next substrate addition event which is typically highly stochastic. The feast-famine cycle is more than just a fluctuating environment, in that the rate of the substrate consumption affects the feast and famine period. Cells starve until the substrate is supplied, and once the environment gets substrate-rich, the cells use it quickly.
During the feast (substrate-rich) period the growth of the cells changes the state of the environment. If cells use the substrate slowly, the feast period lasts longer and vice versa.
In this sense, the feast-famine cycle is one example that cells are not just affected by the environmental condition, but also changing the environment. While the population dynamics without such feedback are well studied [12,13], our understandings with the feedback between the environment and the population are still underdeveloped.
In the seminal Lenski and coworkers experiment of long-term bacterial evolutionary adaptation [14], bacteria have been diluted into fresh media every 24 hours while the growth reaches the saturation after less than 10 hours, hence the cells do experience repeated feast-famine cycles. However, the famine period is rather short that there is no visible death during the period, and the major observed adaptation was to increase the growth rate and decrease lag-time during the feast period, without a significant increase of the death rate during the famine period [15]. Interestingly, however, in a separate experiment, Vasi and Lenski had isolated mutants after the 30 to 49 days starvation [16], and most of them are found to be inferior in competing with their progenitors in fresh medium, but having better resistance to starvation.
One of the five strains clearly showed an inferior fitness to its original strain for a oneday growth competition but had superior survivability in 15-day starvation. These observations suggest that evolutionary adaptation under repeated feast-famine cycles with long-enough famine periods could result in fairly non-trivial results, especially if there is a trade-off between survivability in the long famine period and the growth rate in the feast period. It is also worth mentioning that, the trade-off of growth rate and survivability has been clearly documented for E. coli mutants with genetically manipulated rpoS activity [17]; the rpoS gene encodes for the stationary/starvation sigma factor σ S and known to cause trade-offs in multiple stressed conditions [18,19].
In the previous works, the trade-off between the growth rate and the growth yield are well-documented [20][21][22][23][24][25][26]. Fast-growth but low-yield results in an extended starvation time and a reduction of the population size at the end of the feast period.
It has been theoretically predicted that a fast-growth but low-yield species will wipe out a slow-growth but high-yield species in a well-mixed environment [21,25,26].
However, the effect of the trade-off between the growth rate and the survivability has not been theoretically analyzed yet.
In this paper, we construct a population dynamics model where the population growth is driven by a discrete and stochastic substrate addition events. In order to understand the simplest case, we consider a well-mixed system with a single niche, i.e, there is only one kind of nutrient/substrate that bacteria can consume to grow.
Bacteria cells divide by consuming substrates, and if there is no substrate the cells die at a constant rate. In addition to the growth/death dynamics, mutations take place to change the growth rate. We introduce a trade-off between the growth rate and the death rate or the growth yield: The fast-growing cells have a higher chance to win the competition for the nutrients in the feast period, but instead, are less tolerant in the famine period because of high death rates or a small number of the population due to the low yield.
Using the model, we show that the bacterial population faces the Tragedy of The Commons (TOC)" type phenomenon [27] under feast-famine cycles; the fast-grower tends to take over during the feast period, making the population less tolerant for the famine period. This evolutionary dilemma drives the population to faster growth and eventual extinction in the famine period for the growth-yield trade-off as previously predicted [21], but we find that in the growth-death rate trade-off case, the longterm outcome depends on the functional form of the trade-off. By analyzing the model, we show that this is because the effective fitness under repeated feast-famine cycle is determined by the ratio of the growth rate and the death rate. We then focus on the non-trivial case of the growth-death rate trade-off where functional form drives the population to extinction in repeated feast-famine cycles, and show that the TOC effect can be prompted by increasing the amount of substrate added to the environment. Finally, we discuss the model assumptions and possible extensions in comparison with experimental results in the literature.

MODEL
The model consists of a state vector X defined by (M + 1) integers that denote the population of M species (either genotypes or phenotypes) and the number of substrates in the environment. It is expressed as X = (N 0 , N 1 , · · · , N M −1 , S), where N i is the number of the ith bacteria, and S represents the number of the substrates.
All elements of the state vector X are non-negative integers. Each species can have a different growth rate, growth yield, and death rate.
A single individual of the ith species proliferates at a constant rate µ i being given by µ i = (i + 1)∆µ if S is larger than zero, while it dies at rate γ i under the starving condition. Each species has a different growth yield Y i , and ⌈1/Y i ⌉ substrates are consumed when a single bacterium of the ith species divides where ⌈·⌉ is the ceiling function. The number of the substrate in the environment is recovered to S = S m when the substrate addition event takes place at a constant rate of 1/λ.
We introduce the mutation among species to the model. It occurs with probability ρ when an individual divides, and then, the daughter cell of the ith species becomes either the (i − 1)th or (i + 1)th species in an equal probability.
By assuming that each event occurs as the Poisson process, the master equation for the simplest one species case is given by where P (N, S) is the probability of the state with N bacteria and S substrates.  (1) with the Gillespie algorithm [28]. The oscillation in the number of bacteria results from the feast-famine cycle.

RESULTS
The Tragedy of Commons dilemma in the bacterial evolution under the feast-famine cycle We study the effect of the feast-famine cycle in the multi-species system The master equation for M (M > 1) species system is obtained by just extending Eq.(1) for M species and introducing mutation among the species. There is no direct interaction among the species, but the species interact via the competition for the substrate. For the exact expression, see Eq. (A1) in Appendix. We compare the dynamics of the model with the growth rate-yield (µ − Y ) and the growth rate-death rate (µ − γ) trade-off separately.
We carried out stochastic simulations of the M-species model (Eq.(A1)) with a fixed λ and S m value using the Gillespie algorithm. Fig.1b  arguing that under the growth-yield trade-off, the species with a higher growth rate outcompetes others even though it leads to a reduction of the population size due to the small yield [21,25,26]. On the other hand, in the µ − γ trade-off case, the evolution of the growth rate can either leads to population collapse ( To understand the differences of the outcome, we constructed a simplified version of the model (Eq.(A1)) which allows us analytical calculations. We approximate that the population size is a continuous quantity and consider deterministic growth and death. We denote the number of the ith species right before the nth substrate addition event by a continuous variable N i (n). Also, we regard S as a continuous variable, and thus, remove the ceiling function from the yield.
After the n-th substrate addition, the species grow exponentially until all the substrate runs out. The length of the feast period after the nth addition event, τ (n), , because the increment of the total population divided by the growth yield should sum up with the added substrate S m . If the interval between the nth and the (n + 1)th addition events is longer than τ (n), the cells experience the famine period to die at the rate γ i . Thus, the number of cells right before the (n + 1)th substrate addition event is given by where ∆t(n) is the stochastic variable representing the interval between the nth and the (n + 1)th addition events, which follows the exponential distribution with average λ. By taking average of the effective growth rate, ln(N i (n + 1)/N i (n)), over the exponential distribution, we obtain a deterministic, discrete map system which describes the dynamics of the population growth as In Appendix, we show that the map dynamics has at least M fixed points that only one species exists and the number of cells is zero for the other species, which is given by The linear stability analysis for the fixed points showed that only one fixed point among the M fixed points is stable. The condition for the fixed point to be stable is to have the largest ratio of the growth rate to the death rate µ i /γ i (see Appendix).
In the µ − Y trade-off case, the death rate γ i has no index dependency, and thus, the stability or the fitness is simply determined by the growth rate. There, the growth yield Y i has no effect on the competition among the species, and thus, the bacterial population evolves to increase the growth rate without caring of how the growth is efficient in terms of the substrate consumption.
In contrast, the stability criterion tells us that the functional form of µ − γ tradeoff affects whether the evolution lasts until the whole population goes extinct or not.
When the trade-off is linear (γ = a+bµ), the ratio µ/γ has no upper bound, while the ratio is bounded in the square trade-off (γ = a + bµ 2 ) case. Therefore, in the square trade-off case, once the growth rate reaches the optimal point (the maximum µ/γ), the system stays at that state. On the contrary, due to the lack of the maximum, the growth rate will never stop increasing for the linear trade-off case and leads to the population collapse.
In both the growth-yield trade-off and the linear growth-death trade-off case, the growth rate increases over generation while the tolerance to the famine period gets worse due to the decrease or increase of the yield or the death rate. The evolution then leads to the reduction of total population size. Eventually, the population size becomes too small to tolerate fluctuations of the substrate addition time, and it results in the extinction of the whole population when the substrate addition times are longer than usual by chance. This can be seen as one of the typical consequences of "The Tragedy of The Commons" Dilemma.
Impact of the feast-famine cycle for the survival of the bacterial population In the previous section, we have seen that under a repeated feast-famine cycle, the TOC dilemma is evoked if the trade-off relationship is such that the death rate increase linearly or slower than linear with the growth rate, or when there is a growth rate-yield trade off. In the following sections, we study how the "degree" of the feastfamine cycle affects the survival of the bacterial population and the evolutionary dynamics when the TOC dilemma is evoked. We chose to focus on the less trivial case of linear growth-death trade-off model γ(µ) = a + bµ in the following. Firstly, we introduce the "degree" of the feast-famine cycle. In the following, we change the value of the average famine period λ while keeping the time-averaged substrate supplyS = S m /λ constant. With this constraint, the change in λ, and accordingly in S m , controls the severeness of the feast-famine cycle. A large λ (and S m ) value indicates that a large amount of the substrate is supplied to the environment less often, corresponding to the severe feast-famine cycle. On the contrary, the limit of λ → 0 with keepingS constant would correspond to the continuous substratesupply limit, though strictly speaking in the present model this limit cannot be taken due to the discreteness of S m .
We compared the dynamics under the different degrees of the feast-famine cycle.

Crossover from the directed evolution to neutral evolution
In order to have a better understanding of the observed behavior, we now focus on the cross-over of the survival time T s shown in Fig.2b from being approximately proportional to λ −2 to λ −1 .
In the moderate feast-famine cycle (λ ≪ 10) depicted in the top panel of Fig.2a, it appears that the evolution speed of the growth rate has two regimes: The evolution speed significantly slows down after the average growth rate reaches ≈ 1, which happens around time 1 × 10 7 in this example. Since the mutation probability is constant, we hypothesized that the dynamics of how a new species takes over the majority of the population changes with the average growth rate.
To quantify this change of the dynamics, we studied the dynamics of takingover among species by setting N ini cells of the ith phenotype at t = 0, and ran the population dynamics until the dominant species becomes another.
For the simulations performed in this section, we add the spontaneous migration term to the model. The spontaneous migration makes it possible to increase the number of the cells by one without consuming any nutrient even under the starved condition (for the detailed expression of the term, see Eq.(A1)). An introduction of the term is a mathematical treatment for gathering many trajectories efficiently for better statistics. Since the spontaneous migration rate ǫ is set to be sufficiently smaller than any other parameters in the model, it is effectively the same as introducing one cell into the system when the population extinct.
From this computation, we obtained the transition probabilities from the ith species to another. With our default parameter set, it never happened that the species other than the nearest neighbors of i becomes dominant before i−1th or i+1th dominates the system. Thus, the obtained transition probabilities were for increasing the growth rate by ∆µ (probability p(µ)) or decreasing it by ∆µ (probability 1 − p(µ)). Fig.3a shows the asymmetry of the probability for increasing/decreasing the growth rate, defined by the difference between the two probabilities (2p(µ) − 1).
As clearly seen, the evolution of the growth rate takes place in a directed manner up to µ ≈ 1, whereas the dynamics of the evolution resembles the random walk when µ ≫ 1. This qualitative difference in the evolution dynamics above and below µ ≈ 1 may consistently describe the crossover of the survival time T s , which is happening at around µ c ≈ 1. Namely, when the critical growth rate µ c is below one for long enough λ, the extinction happens relatively quickly since the growth rate systematically increase through the evolution, but when µ c is above one, the evolution takes a lot longer time due to the diffusive behavior, hence survival time T s grows faster as decreasing λ.
Where does this transition from the directed to neutral evolution come from?
The simple map dynamics (Eq.(4)) just tells us that the species with the largest µ/γ dominates the population which does not explain the random walk-like behavior of the averaged growth rate. In order to gain more insights, we studied the detailed dynamics of the take over events of the population by dominant species in stochastic simulation.
Since the main purpose of this simulation was to ask how the dominant species changes one to another, we simulated the model with only two species which have a slightly different growth rate to each other. One has a growth rate µ l and the other has a higher growth rate, µ h given as µ h = 1.05 · µ l . Fig.3b shows the timeaveraged populations of the species are plotted against the growth rate of the slowlygrowing species (µ l ). The fast-growing species dominates the whole population and the number of individuals is much larger than that of the slow grower in the small µ l region, whereas the difference of the population sizes of the two species shrinks as µ l increases and it gets indistinguishably small at µ l ≈ 1. Examples of time courses are plotted in Fig.3c. The dynamics with large µ l shows rather stochastic changes between the fast-grower dominating and slow-grower dominating states, while with small µ l the fast grower is stably dominates the system. This shrinkage of the gap in the two populations explains the transition from the directed to the neutral evolution of the growth rate. Intuitively, the mechanism of this shrinkage can be understood by considering the effective fitness µ/γ. Since we use the linear trade-off γ(µ) = a + bµ, the difference of the effective fitness between the fast species with the growth rate µ h = (1 + δ)µ l and the slow species with the growth rate µ l is given by which approaches zero as µ l increases. In other words, the larger the value of µ l is, the harder it becomes for the fast species to take over the population. As a result, the growth rate performs almost a random walk through evolution for the large value of µ l .
The more quantitative understanding can be obtained by applying the Wright-Fisher (WF) model [29]. The WF model is a stochastic model describing temporal changes of the population structure such as the fixation probability and the fixation time. While the set-up of the present model does not fully fit the WF framework, we can apply the framework to the model with some assumptions which are described in the Appendix. The WF framework enables us to calculate the probability of the fast grower to be fixed in the population under no-mutation no-migration condition.
The fixation probability is given by where N t , y and u represents the total number of the cells, the initial fraction of the fast growers, and the relative fitness of the fast grower defined as  5)), we can see that the decrease of the fixation probability is led by two effects, namely, the decrease of the relative fitness advantage and the population size-effect.
One effect is the form of u being a decreasing function of µ, hence as discussed before the advantage of the fast growth is reduced even the population size stays constant.
In addition, the population size shrinks as the growth rate increases, and it makes the population dynamics noisier, making the small fitness difference no longer be the determinant of the dynamics. The two effects similarly contribute to the change of the fixation probability as shown by the dashed lines in Fig.3d, where the fixation probability Eq. (5) with a constant relative fitness u or a constant total population N t are also plotted.
The WF model also in principle shows the parameter dependence of the crossover point, which should correspond to the point where fixation probability is sufficiently close to 1/2 when starting from the equal population (y = 1/2). The closed-form is difficult to obtain because the complex dependence of N t on µ, but the form indicates that the crossover growth rate depends on the trade-off function parameter values (a and b).
1 To compute the fixation probability of the present model in the stochastic simulation, we set the mutation rate and the migration rate to zero, and run the dynamics from the fixed initial value of N l , N h , and S. A single run finishes if one of them extincts, and it is repeated with different random number seeds to compute the fixation probability.  The fast grower is expected to fail the fixation (1 − p fix )/p fix times on average, and new mutant need to appear at every fixation failure. Therefore, the time for the dominant species to become n to n + 1 is given by Note that p fix , τ h , and τ l are the functions of µ n and µ n+1 . In addition, at the every change of the dominant species that increase the average growth rate, the average population size decreases. This population size was already calculated in Appendix D from the master equation, and we assume that the extinction occurs when this population size becomes smaller than unity. Then, by summing up T n,n+1 over n, until the population size reaches to unity, we obtain the estimate of the survival time.
With the present parameter values shown in Fig.4, the dominant term of n T n,n+1 is the time for the fast grower appearance by mutation, i.e., Σ(1 − p fix )/p fix · τ m . The comparison between the simulation and this expression is compared in Fig. 4 shows a reasonable agreement. Also, the agreement indicates that the reduction of the survival time T s is mainly due to the increased rate of getting a fast-growing mutant with increasing S m because it increases the typical population size. The more detailed comparison with the rest of the terms is given in Appendix F.

DISCUSSION
Here, we developed a stochastic population dynamics model in which the vital substrates were supplied to the environment by discrete, stochastic events rather than continuously. This stochastic substrate addition separates the dynamics into two phases, namely, the feast and the famine phase. During the feast period with plenty of substrates, the cells with a higher growth rate increase their population more quickly than the others. On the other hand, during the famine period, the cells could not grow but die due to the lack of substrates. The survival time, or the average time to extinction, was shown to have a powerlaw dependency to the average waiting time λ with the constant time-averaged substrate supplyS = S m /λ, with a cross-over from close to λ −2 to λ −1 . The cross over stemmed from the transition from the directed evolution to the undirected, random walk-like evolution dynamics. As the average growth rate increases, the difference in the fitness between two species with similar growth rates reduces. In addition, the population dynamics become noisier due to the decrease in the average population, and thus, the difference in fitness becomes less influential to the dynamics. Finally, it was shown that a pure increase of supplied substrate per nutrition addition event enhanced the extinction, and the reduction of the mutant appearance time due to increased population size was turned out to be the main part of the decrease of For the cross-over from the directed to the random evolution to occur, the form of the trade-off is probably essential. With the square trade-off, one can choose the parameter values so that the evolution stops without enough shrinkage of the fitness gap and the population size, and then, only the directed evolution is expected.
On the other hand, with the trade-off without the upper bound of the fitness, the population size becomes considerably small before the extinction. Also, the fitness gap always shrinks as the growth rate increases. By this argument, a cross-over being similar to what we presented above is expected to occur.
It is worth noting that the trade-off between the growth rate and death rate (killing rate) are well reported under a variety of stress conditions for bacteria [17-19, 22, 30-33], and the trade-offs between the growth and death rate are linear in some cases [34,35]. Also, the linear relationship between the growth-and the death rate is documented in the continuous culture of the fission yeast [36]. With the linear trade-off, the evolution could not find an optimal point in the present model and the whole population faced extinction.
Undoubtedly, bacteria do not go extinct so easily as the present model has shown.
Indeed, it is reported that the long-term evolution experiment of E. coli leads to the coexistence of the variety of bacterial genotypes [14,37,38].Since the present model is a simple toy model, there is a large gap between the model and the real experiment. For the further investigations of bacterial population dynamics under the feast-famine cycle, it might be fruitful to argue the differences between them.
First of all, the present model has no spatial degree of freedom. It is well-known fact that the introduction of the spatial structure will allow replicator models to have coexistence solutions [39][40][41][42][43][44]. Introduction of the spacial structure was also proposed to avoid the TOC caused by the growth-yield trade-off [21]. Furthermore, it was experimentally shown that E. coli cells with the attenuated rpoS outcompeted the wild-type cells in the well-mixed culture but coexisted with the wild-type in the spatially-structured culture [45]. With the spatial structure introduced, it is possible that the present model shows the coexistence of multiple species that work to prevent extinction due to the TOC scenario. It is also worth noting that, in reality, a small fraction of the bacterial population can be non-growing persisters [46][47][48]. The introduction of the persister phenotype makes the big jump between the high-growth and high-death and low-growth and low-death state possible, and it might change the evolutionary strategies. Also, the lag and stationary phases were not implemented in the model. The period of the lag phase is considered to increase with the starvation time (or the time in the stationary phase) [11,49,50], and the prolonged lag time is clearly disadvantageous for the competition for substrates. The real bacteria might design their growth strategy also to cope with the length of lag time. Furthermore, it has been shown that during the death phase, the substrate provided by dead cells are utilized by alive cells to survive longer [51], and this feedback from cell death to the environment can be another factor to be considered to compare with the reality.
Another factor to consider is the stochasticity in the fixation process. It is theoretically shown by using a Wright-Fisher-type model that the trade-off between the fitness advantage in WF scheme (e.g., the growth rate) and the carrying capacity (the maximum population of given species) makes a fixation probability of the species with smaller fitness advantage higher than that of the species with bigger fitness advantage under a certain condition [52]. When a single cell of fast-growth and small carrying capacity is introduced into a community formed only of the other type, the chance of the former species to be selected is small due to a small frequency in the community. By contrast, in the opposite case, the population size of the community is small, and thus, the frequency and the fitness advantage are less influential to the selection. There, the randomness dominates the selection process rather than other factors.
The effect of the feast-famine cycle to the bacterial evolution has been addressed experimentally, though typically the famine period is the order of a day or shorter that the bacteria enters the stationary phase but not the death phase. In Lenski and coworkers evolutionary adaptation experiment [14,38,[53][54][55][56] has been continued for more than 30 years and a variety of mutations have been confirmed. For instance, some of the new genotypes can use citrate as a sole carbon source which the ancestral strain was unable to utilize [54,56] and cross-feeding polymorphism to appear [55].
Those observations mean that new niches were created during the evolution time series. As only a single nutrient is considered in the present model, there is only one niche and the possibility of the cross-feeding was not taken into account. An extension of the model to have more than two types of the nutrient and secretion of the chemicals from the cells will allow the model to have more niches and hence polymorphism to appear.
In a recent experiment of bacteria under repeated feast-famine cycle with the exchange of medium [11], a large amount of the bacteria cells were also flushed out when the fresh media is added. In this case, the cells that form aggregates are selected, even though the growth was in the liquid culture. Though this particular experiment imposed the selection for aggregates, in general one should have in mind that the aggregation and hence spatial heterogeneity can arise even in a liquid culture [57], which can make a deviation from the prediction of a "well-mixed" model.
Having these possible deviations from the present model in mind, it will still be interesting to perform an evolution experiment with longer famine period, to see the role of the trade-off between the growth rate and the death rate in the repeated feast-famine cycle.
As discussed, the present model has plenty of choices to be extended for emulating the strategy of the real bacterial population. But in spite of its simplicity, it provides several insights into how the feast-famine environment could affect the bacterial population dynamics, which hopefully helps the future development of our understandings of bacterial population dynamics and evolution.

ACKNOWLEDGMENTS
The authors thank Nen Saito for fruitful discussions and Hiroshi Kori for suggesting the reference [17]. This work was funded by the Danish National Research Foundation (BASP: DNRF120). • ǫ : migration (invasion) rate (normally it is set to zero) • λ : the average waiting time for the substrate addition , where IE i is the step operator of species i which acts to an arbitrary function f (· · · , N i , · · · ) as IE i f (· · · , N i , · · · ) = f (· · · , N i + 1, · · · ). Here, we allow to P (N, S) with negative S to have non-zero value because instead of setting the boundary at S = 0 because it needs an exception handling, for instance, assuming that all the transition from such states go to the state with S = 0. While the state S < 0 is reachable, the right hand side of the equation is the same among S ≤ 0. Therefore, by assuming the sum of the probabilities for S ≤ 0 represents the probability of the starving state, we can avoid the introduction of the exception handling.

Appendix B: Derivation of the map system
Here, we derive the map system, first for the one-species case. To derive it, we ignore the interactions among species (i.e., deal with the growth/death dynamics of only one species.) and stochasticity in the bacterial growth/death dynamics. First, we calculate the duration of τ in which the bacterial population use up the newly added S m substrates. Since the bacterial population grows exponentially as long as the substrate remains, the number of the bacteria at time t,Ñ(t) is given as where we set t = 0 as the time the substrates are newly added. The integral of this equation from t = 0 to t = τ gives us the cumulative consumption of the substrate given by (e µτ −1)Ñ(0)/Y . By the definition of τ , the cumulative consumption equals to S m , and thus, we get Next, we introduce the waiting time between nth and (n + 1)th substrate addition periods, ∆t(n). Note that if τ (n) < ∆t(n), the bacteria run out all the substrates and start to die, otherwise, the population grows toÑ (τ ). Thus, the population at the (n + 1)th substrate addition event is determined by the population at the nth addition event as follows; where N(n) indicates the number of the bacteria at the nth substrate addition. Since τ depends N(n), we put (n) to explicitly show that τ can have different values for different n's.
Note that the average of ln(N(n + 1)/N(n)) over ∆t(n) with the exponential distribution gives the effective growth rate and leads to the deterministic time evolution of the bacterial population. The average results in For the multi-species case, the map system is derived in a similar way. It is given as where τ (n) is given as the solution of the following equation; Although we performed many numerical computation of the map dynamics with a variety of parameter choice, any attractor at which the species coexist was not found. Thus, we concentrate to the attractors at which only one species has nonzero population. By settingμ i (n) = 0 and N j (n) = 0 for any j = i, we get the steady solution given by us suppose that there are M species. The Jacobi matrix of the system at the nth substrate addition event (N 0 (n), N 1 (n) · · · , N M −1 (n)) is given as By renumbering the species index, we can assume that only the species with index 0 survive at the attractor without losing generality. Since N 0 (n → ∞) = N st 0 and N 1 (n → ∞) = · · · N M −1 (n → ∞) = 0 hold, the elements of the Jacobi matrix get simplified at the fixed point as where ·| st indicates the value of the functions at the steady state which we are interested in. This equation shows that the Jacobi matrix is triangle, and thus, the eigenvalues are The first eigenvalue 1 + N st 0 λ ∂μ 0 ∂N 0 st. is zero, andμ i | st is given asμ i | st = µ i (γ 0 /µ 0 − γ i /µ i )/(1 + γ 0 /µ 0 ) because τ is given by λ ln(1 + γ 0 /µ 0 ) holds at the fixed point.
Since λ is positive constant, γ 0 /µ 0 < γ i /µ i , for all i = 1, 2, · · · , M − 1 is the stability condition for this fixed point. Therefore, the only one fixed point with the largest µ/γ is stable among all the fixed points at which only one dominant species exists. The full-model is approximated by assuming that λ ≫ ln(1 + S m /N st )/µ holds to ignore the terms P S (N)/λ for S > 0. Here, we also introduced a spontaneous migration with rate ǫ, because the master equation has an absorbing state at (N, S) = (0, S m ) and without the migration term, only the steady state is P (0, S m ) = 1 and 0 for the others. The approximated master equations are given as Here, we introduce the moment-generating function for each S defined as We get following equations at the steady state where G ′ S represents the first-order derivative of G S 's respect to z. By setting z = 1, we get Thus an equality holds.
Also, by differentiating Eq.(D1) respect to z again, we get By summing up these S m equations and setting z = 1, we get From Eq.(D3) and (D5), the first moment is given by Since G 0 (1) is indeterministic from the equations (D2)-(D5), we make an anzats that G 0 (1) has the form with Here, we use the steady state solution the map dynamics in the main text Eq.(4) aŝ N and ⌈·⌉ represents the ceiling function.
At the end, the steady state solution of the bacterial population is given by We made the anzats that G 0 (1) has the form shown in Eq.(D7) because G 0 (1) represents the fraction of time that the system stays at the states with no substrate.
Apparently, the timescale to escape from the no-substrate states is λ. On the other hand, the length of time that the system stays at the with-substrate states is given by H 0 /µ. Suppose that the system is at the state with (N, S) = (n, S m ). The spontaneous immigration rate is ignorable if there are the substrates, and thus, the bacteria proliferates without experiencing any other elemental processes until all the substrates are run out. Since the number of bacteria increases strictly one by one as a single substrate is used, the timescale for reaching to the non-substrate states is obtained by summing up the timescales of the proliferation with given number of the substrate. It is given as (1/n + 1/(n + 1) + · · · + 1/(n + S m − 1))/µ. The number of the bacteria right after the substrate addition is typically given by ⌈N⌉, and thus, the sum is written as H 0 /µ.
Lastly, there is finite possibility that the system reaches the state (N, S) = (0, S m ) which becomes the absorbing state under ǫ → 0 limit. Once the system falls into the state, the only one elemental process could take place is the spontaneous immigration occurring at rate ǫ. It might be natural to estimate the probability of the system to fall into this state from the ratio of two timescales, namely, the timescale of the substrate addition and cell death. A single bacterium in N bacteria dies within 1/γN on average. Thus, the whole population is expected to extinct within H 1 /γ. Therefore, the probability that the total population becomes zero before the substrates being supplied can be approximated by (exp[−H 1 /γλ]).
The sum of these three parts gives the typical length of time consumed in between two substrate addition events, and thus, the ratio between λ and the sum corresponds to what we desired. WF model, while detailed descriptions can be found in standard textbooks [29].
The WF model is a Markov-chain model which predicts the temporal evolution of relative abundances of the species. Suppose that there are two species, species h and l (high and low), which has the relative fitness 1 + u (here we assume u > 0) and 1, respectively. The total population is kept to the constant number N t , and the generational change is assumed to take place simultaneously. The probability that one individual at the next generation to be the species h is given by in one generation is given by If the mutation is not included in the model, the model has two absorbing states at which only one species dominates the whole population corresponding to the fixation of the species. By using the diffusion approximation of the WF model, the fixation probability of species h is analytically given as p fix = (1 − exp(−N t uy))/(1 − exp(−N t u)), where y is the initial abundance of the species h relative to the total population.
For fitting the present model into the framework of the Wright-Fisher (WF) model, we have to make several assumptions, while they are not always fulfilled by the present model exactly. Here, we first list the assumptions which we need to make with comments.
(i). The total number of the population is constant : Clearly this requirement is violated, while we assume that the total number of the population is fixed at the steady-state value (Eq.(D8)) with the growth rate of the slow grower.
(ii). Proliferation/Cell death takes place simultaneously : Since cell division and cell death occur as a Poisson random event in the present model. We can neither fulfill the two conditions. We assume, however, that it is effectively satisfied by coarse-graining the time. We suppose that on average every Nλ/S m time interval, N individuals die and N bacteria are newly born, and accordingly, one generation in the WT framework corresponds to Nλ/S m in the timescale of the present model.
(iii).The fitness function : Due to the differences between the WF framework and the present model, it is unclear which parameter in the present model should correspond to the relative fitness advantage. So, here we adopt the ratio between the growth rate and death rate as the fitness function because it is already shown to determine the stable fixed point under the deterministic approximation of the model.
Here, we calculate the fixation probability and the fixation time of the two-species WF model. We assume that there are two species with slightly different growth rates and death rates. We do not include the mutation process in the analysis. Instead, we set a very small portion of the total population (y) as the bacteria with higher fitness. The fixation probability of the bacteria with the higher fitness is given as under the continuous (large N) limit, where α is defined as α = uN. Since one of the two species is eventually fixed, the fixation probability of the other species is given as 1 − p fix .
Next, the fixation time τ i (i = l, h) is given by solving the backward Kolmogorov equation below, The solutions are given by the following equations With the sufficiently small mutation rate ρ, the survival time can be estimated from the recursive taking-over process between two species with different fitness values because there are only two neighboring species in the system most of the time. So, we calculate the average length of time needed for the (i + 1)th species dominates the whole population which was initially dominated by the ith species. We denote this as T i,i+1 , and consecutively sum up them to get the total length of time.
The time for the (i + 1)th species to take over the population consists of three A single mutant with the higher fitness appears within τ m = (1−(1−ρ/2) N ) ≈ 2/ρN generation on average. Note that we can ignore the appearance of the mutants with lower fitness because they typically fail to increase their portion in the whole population due to the low fitness and the small population. Thus, the result does not alter even if the (i − 1)th species is dealt as the ith species. In the following argument, we assume that the initial portion of the bacteria with higher fitness, y, equals to the inverse of the population size, N, because it is less likely to happen that more than two mutant appears at the same time due to low mutation rate.
The expected number of failures of the fixation (the count of the event that a mutant appeared but eventually eliminated) is given by R(u(µ l , µ h ), N(µ l , S m )) = R(µ l , µ h , S m ) (F1) ≡ Π l (u, N, 1/N)/Π h (u, N, 1/N) (F2) = coth(u/2) sinh(uN/2) − cosh(uN/2) e −uN/2 (F3) At every failure, the population has to wait for a next mutant to appear and one failure takes time with length τ i,i+1 l . The fixation of the species with higher fitness has to take place only once, and thus, the population is expected to evolve their dominant species from ith to (i + 1)th species in Lastly, we need to convert the timescale because the unit of time is the number of generations in the WT model. Since S m bacteria divide after a single substrate addition event, N/S m substrate addition events are needed for a single generation to pass (N bacteria divide). The substrate addition events take place every λ interval on average. Thus, the timescale is converted from WT to the present model by multiplying N(µ)λ/S m .
Note that here we assumed that the (i+1)th species always has a fitness value higher than that of the ith species (u > 0). This assumption makes it unlikely to happen that the (i − 1)th species appear by mutation and eliminate the ith species. This directedness is vital to make an estimate of the survival time by simply summing up the length of time for each taking-over event T i,i+1 . The assumption is violated if the fitness function µ/γ(µ) has a peak, for example, the fitness function with the square trade-off. [2] Jacques Monod. The growth of bacterial cultures. Annual review of microbiology, 3,371, (1949).