Supplementary materials for : Population cycles and species diversity in dynamic Kill-the-Winner model of microbial ecosystems

Determinants of species diversity in microbial ecosystems remain poorly understood. Bacteriophages are believed to increase the diversity by the virtue of Kill-the-Winner infection bias preventing the fastest growing organism from taking over the community. Phage-bacterial ecosystems are traditionally described in terms of the static equilibrium state of Lotka-Volterra equations in which bacterial growth is exactly balanced by losses due to phage predation. Here we consider a more dynamic scenario in which phage infections give rise to abrupt and severe collapses of bacterial populations whenever they become sufficiently large. As a consequence, each bacterial population in our model follows cyclic dynamics of exponential growth interrupted by sudden declines. The total population of all species fluctuates around the carrying capacity of the environment, making these cycles cryptic. While a subset of the slowest growing species in our model is always driven towards extinction, in general the overall ecosystem diversity remains high. The number of surviving species is inversely proportional to the variation in their growth rates but increases with the frequency and severity of phage-induced collapses. Thus counter-intuitively we predict that microbial communities exposed to more violent perturbations should have higher diversity.

When viewed over a long period of time any given species would repeatedly cycle between low and high population numbers. Such cyclic dynamics of populations of individual species masked by an approximately constant total population saturated at the carrying capacity of the environment is discussed in the ecological literature as "cryptic cycles" [22][23][24] .

Model
Consider a number of bacterial species/strains sharing the same environment and competing for the same rate-limiting nutrient defining its carrying capacity. Their populations sizes at time t are denoted as P i (t), where i = 1, 2, … , N. Each of these individual species is exposed to rare but severe collapse events in which its population is suddenly and drastically reduced by a constant factor γ ≪ 1. We assume that these collapses happen relatively rarely so that the total population of all bacterial species has sufficient time to reach the steady state value given by the overall carrying capacity of the environment. Without loss of generality carrying capacity can be set to 1, so that in between collapses one has ∑ j P j (t) = 1. In our model we assume that while the total population of all stays constant, relative population sizes of individual species continue to change exponentially in-between collapse due to differences in their fitness in the saturated environment.
In the spirit of Kill-the-Winner principle we assume that the rate of collapse of the species i is proportional to its population size P i . Due to a broad distribution of population sizes this rule strongly biases collapses towards one or few largest populations. We assume that collapse events are independent of each other, so that the time interval between consecutive collapses is exponentially distributed with mean τ.
One update cycle in our model consist of three steps: (1) Draw a time interval Δ t until the next collapse event from the exponential distribution P(Δ t) = exp (− Δ t/τ)/τ. (2) Calculate population sizes at the time of collapse. In between collapse events relative population sizes are assumed to change exponentially while the total population stays saturated at 1 (the carrying capacity of the environment): (3) Select one species to collapse with the probability equal to its relative population size P i (t + Δ t) and multiply its population by γ.
In our simulations each of N species is assigned its individual growth rate drawn from the Gaussian distribution with zero mean and standard deviation σ. The value of the mean is not important since normalization of the overall population to 1 ensures that only relative growth rates matter. Furthermore, the exact form of the distribution of growth rates is not particularly important. In our mathematical analysis we will use a more convenient exponential distribution of growth rates: P(g i ) = exp(− g i /σ)/σ, while delegating more cumbersome derivations for the Gaussian P(g i ) to supplementary materials.

Results
Collapses supports Diversity. Figure 1A shows a typical outcome of a simulation of our model with γ = 10 −3 and σ = 4 over around 100 population collapses after which only D = 3 out of N = 8 species survive. The relative growth rate g i of the species is the main predictor on whether it will survive or not. Indeed, as shown by the rainbow coloring of curves in Fig. 1A ranging from dark red (the slowest growing) to purple (the fastest growing) the 3 surviving species have the largest values of g i .
A natural question to ask is what determines the number of surviving species/strains in the steady state of the model? In the limit of very rare collapses the fastest growing species would diverge from the rest of the population so much that it will be the only one to survive, as indeed expected from the competitive exclusion principle 1 .
The situation is more complex for intermediate rate of collapses where more than one of the fastest growing species can coexist with each other but some of the slowest growers become extinct. In the steady state each of these surviving species repeatedly cycles between low and high populations. Faster growing species reach large population sizes more often which makes them to collapse more frequently thus eliminating their growth advantage. As we show below this balance can be sustained within a finite range of growth rates.
For each of the species its individual growth rate g i is reduced by the same negative number − g cc due to the overall resource competition quantified by the denominator in Eq. 1. In the steady state the excess growth rate of each of the surviving species (g i − g cc ) must be exactly compensated by the logarithmic population losses |log γ| due to collapses happening at the species-specific probability c i : Note that as the probability of collapse (per update) c i = (g i − g cc )Δ t|log γ| needs to be positive and normalized. Positive values of c i means that only the fastest growing species with g i > g cc would survive in the long run. The collapse rates of these D surviving species are further constrained by normalization ∑ = = c 1 j D j 1 , reflecting the requirement of one collapse per update. Using eq. 2 the threshold g cc is then determined by: For a given set of species, this allows us to self-consistently calculate g cc and D.
For g i selected from the exponential distribution with standard deviation σ the diversity D is given by (see Supplementary Information for This expression holds for average diversity provided that it is larger than 1. This is because a single fastest growing species would always survive. Clearly D is also capped from above by N. Similar relation holds for the Gaussian distribution of growth rates and is in agreement with our numerical simulations of the model shown in Fig. 1B. For the exponential distribution the growth rate threshold above which a species survives is given by g cc = σ log(N/D) = σ log(NσΔ t/|log γ|). Note that while threshold explicitly depends on the starting number of species, the final diversity given by Eq. 4 is independent of N. This particular property of the exponential distribution would be modified for other distributions resulting in a mild dependence of D on N.
A Gaussian distribution of growth rates would slightly increase the diversity compared to the exponential distribution with the same spread, while a more fat-tailed *(say, power law) distribution would decrease it.
Our basic model can be generalized to the case where different species have different collapse ratios γ i . This may for example reflect their different degrees of vulnerability to phages, or different ways to partition their population in physical space. The only consequence of this modification is that log γ in the equations above needs to be replaced by its average value across species (see supplementary materials for simulation results).
In our model the collapse probability of a given species is proportional to its population size. Thus time-averaged relative population size of each of the species species is equal to its overall collapse frequency 〈 P i (t)〉 t = c i . This is consistent with "Kill-the-Winner" principles according to which species with larger populations collapse more often. Rainbow colors correspond to growth rates ranging from the slowest (dark red) to the fastest (purple). For each of the species, the likelihood of collapse is proportional to its population sizes ("Kill-the-Winner" rule) and the collapse ratio γ = 10 −3 is the same for all species. Only 3 fastest growing species survive in the long term. (B) The final diversity D counted as number of surviving species as function of σΔ t -the spread of growth rates integrated over the average time between collapses. Each black dot represents the outcome of one simulation started with N = 500 species exposed to collapse ratio γ = 10 −6 . The dashed line is the analytical fit similar to Eq. 4 but here done for the Gaussian distribution of growth rates used in the simulation (see Supplementary Materials for details). Figure 2B illustrates this cyclic dynamics in a system containing a mixture of slow and fast growing species. Surviving populations mostly grow, but do so at different rates. Their coexistence is possible only because of the negative feedback via "Kill-the-Winner" rule where populations of an individual species get severely reduced once it starts to dominate the overall biomass. The population of each of the species goes through approximately periodic cycles of growth and collapses with the period T i = 1/c i = |log γ i |/(g i − g cc )Δ t (in units of collapse events). Thus the slowest surviving species (marked blue in Fig. 2B) nearly never collapse, whereas the fastest growing species (marked red in Fig. 2B) obtain dominance and expose themselves to a collapse on a much shorter timescale. Individual collapse events of these species are marked in Fig. 2B with red and blue arrows correspondingly. Note that the population of the slowest growing species often decreases not due to a phage-mediated collapse but simply because it gets temporarily outgrown by other species with a faster growth rate.
For comparison in Fig. 2A we show a system of the same size (D = 10) but where all species have exactly the same growth rate. In that case the system has a very long memory of the initially imposed order of species populations, because even after a long time each of the species would collapse the same number of times. That is, if one species have experienced one more collapse than the others, it would be smaller by a factor γ and thus be much less exposed to subsequent collapses until it would regrow to the size where it again may collapse with a non-negligible probability. Indeed, populations shown in Fig. 2A follow much more regular oscillatory dynamics than those with unequal growth rates shown in Fig. 2B. Model with collapses to a fixed threshold. In our standard version of the KtW model the collapsing population is reduced by a constant factor: P i → P i ⋅ γ. An alternative possibility is that following a collapse the population starts at a fixed small threshold value γ irrespective of its earlier population size. This would be the case e.g. when following a collapse the local population is completely eliminated and is reintroduced by one individual from a neighboring region. It can also happen when a collapse drives one species extinct only to be quickly replaced by a single bacterium of a new species. In the thus defined fixed threshold kill-the winner model (KtWT) the diversity remains close to what was reported in Fig. 1B (data now shown). The dynamics is also characterized by individual population undergoing cycles of duration defined by their relative growth rates much similar to what is shown in Fig. 2 for our original model. However the long term memory of cycle order is reduced compared to the constant factor model discussed above, simply because every collapse completely erases the population history of the collapsed species. In what follows we explore the dynamical properties of the fixed threshold model and its variants.
"Kill-the-King" Model. To better understand the cyclic dynamics in the KtWT model we first consider its extreme and deterministic version in which the next collapse always happens at the largest population. We will refer to this version as Kill-the-King (KtK) model. Furthermore, we assume that the growth rates g i of all species ) each with growth rates selected from the Gaussian distribution with σ = 3, and identical collapse ratios γ i = 10 −6 . The blue arrow and the red arrows mark times for collapse events of these two species. Note how the fastest growing species (red line) collapses much more often than the slowest growing species (blue line) which only collapsed once during the time shown. The growth rate difference between these two species is g max − g min = 2.44.
Scientific RepoRts | 7:39642 | DOI: 10.1038/srep39642 are equal to each other. Thereby the asymptotic dynamics becomes periodic with period N when time is measured by collapse events.
To concentrate on slow trends in population size dynamics we only measure them between intervals where each population collapsed exactly once, which in KtK secure that the order of populations is exactly preserved. We relabel species in the order of decreasing population sizes and calculate the ratios δ i = P i+1 /P i < 1 between successive population sizes (The ratio δn for the currently smallest population Pn is defined by its value after the next collapse when it becomes the second smallest). As shown in the supplementary materials, in KtK model these ratios evolve according to the following discrete equation describing changes acquired after a full round of N collapses so that each member of the population collapsed exactly once: The steady state of the equation is reached when all ratios are equal to each other, i.e. δ i+1 − δ i = 0. In this case the logarithms of population sizes are equidistantly distributed in the interval of length |log γ| so that δ i (∞ ) = δ * = γ 1/N . Figure 3 shows a simulation of KtK model with N = 10 and γ = 10 −6 One can see how it asymptotically approaches this steady state.
The asymptotic dynamics of KtK is described by the discrete Eq. 5 which for large N can be approximated by a continuous PDE (see SI for more details) in which the continuous coordinate x replaces the species rank i: Here δ(x) has periodic boundary conditions over x-interval [0, N]. As its discrete counterpart this equation describes the state of our system every N'th time-step. This equation is closely related to the Burgers equation 25,26 , although it differs in terms of the diffusion coefficient that instead of being constant as in refs 25,26 and is proportional to δ(t).
Having finished with Kill-the-King model we return to Kill-the-Winner fixed Threshold (KtWT) model. In the KtWT model population collapses do not always happen in the order dictated by their relative sizes. This results in a somewhat chaotic dynamics illustrated in Fig. 4A. When a smaller population collapses out of turn it causes only a very small rescaling of other population sizes. The (very likely) subsequent collapse of the largest population leads to a situation where these two just collapsed populations become nearly equal in size (δ  1). This dramatically increases the likelihood for further re-orderings between these two species, resulting in an extended period where these two species fight for dominance. This intermittent dynamics switching the order of populations is clearly visible in Fig. 4B with D = 3. The nearly vertical lines clustered around collapse events 5100, 5200, and 5400 correspond to frequent shifts in the population order of three species within the cycle.
An intermittent region ultimately ends with a particular order winning over. After this all populations slowly relax back to the steady state with equal ratios δ * (curved lines ending in horizontal plateaus in Fig. 4B). The exact form of the relaxation to the steady state is derived in supplementary materials. While δ(t) ≫ δ * the relaxation is proportional to 1/(t/N). The expected number of collapse events for iδ(t) to go from ~1 to ~δ * is ~N/δ * or about 300 for the parameters used in Fig. 4. When δ δ  ⁎ t ( ) the relaxation crosses over to δ(t) − δ * ~ exp(− δ * t).

Discussion
Here and before 27 we investigated the impact of severe and sudden population collapses on ecosystem composition and diversity. This approach is complementary to a more traditional description of ecosystem dynamics at or around the steady state solution 2,5,9,28 . The emergent cyclic dynamic in our model is entirely collapse-driven and For clarity we show the state of the model only every 10 collapses, that is to say, after each species collapsed exactly once so that the population order is maintained. The steady state of KtK model where all ratios δ i between rank-ordered populations are equal to each other and to γ 1/N is approximately reached already after 300 collapses. The relaxation to this steady state is described by the discrete anisotropic Burgers equation 5 or its continuous counterpart Eq. 6.
The key assumption used in our study is that larger populations are more exposed to sudden collapses than the smaller populations. This is the foundation of "Kill-the-Winner" (KtW) principle proposed in ref. 5. The resulting negative (or stabilizing) frequency-dependent selection promotes the ecosystem diversity even in the simplest case considered above, where species interact with each other only via competition for a single rate-limiting resource. This KtW bias is very important as it shifts the collapse-driven dynamics away from "diversity waves" we reported before 27 towards population cycles investigated in this study. Indeed, as demonstrated in ref. 27 a version of collapse-driven dynamics in which the likelihood of a collapse is uncorrelated with population size or even biased towards smaller populations (Kill-the-Looser model) results in ebbing and flowing species diversity and bi-modal distribution of species abundances. This should be contrasted with constant diversity predicted in the KtW, KtWT, and KtK models studied here. In our scenario diversity is maintained as population reaches dominance and collapses in a particular order. The species abundance distribution in the here presented models is not bi-modal but uniform on the logarithmic scale (data not shown).
The model in this paper focus on one particularly lethal aspect of density dependent selection, and analyze it in detail. A key result is eq. 4, that quantify a dynamical maintenance of diversity, through frequent collapses of the largest populations. The falsifiable prediction is that higher flux of new phages makes well mixed microbial ecosystem more diverse. The obtained diversity is obtained by preferential attack on the largest population, whereas the traditional steady state KtW emphasize the coexistence of slow and fast growing bacteria in presence of phage 5 . Our argument for targeting the largest population preferentially is that new virulent phages tend to induce larger collapses for host population that facilitate more effective spreading of the phage. In effect our scenario is similar to the assumption that larger population densities have larger effective R0-factors for an "epidemic".
Our approach is a simplified approach to analyze an open, yet well defined system. Apart from assuming the preferential targeting of large populations, it is formulated in the limit where the total population reach carrying capacity before new invasions take place. This last assumption can easily be relaxed by assigning a collapse rate/ collapse size that depend on population sizes also before the sum of all populations reaches carrying capacity. We tested that such increasingly frequent collapses also allowed for higher microbial diversity, and found that it open for even higher diversity than predicted by eq. 4. This is because the bacteria are then so violently exposed to phages that they never sense their mutual competition at saturation. To test how sensitive are our results with respect to introduction of other types of interactions between bacterial species as well as to a more branched topology of Phage-Bacterial Infection Networks (PBIN) we simulated a variant of our model where in addition to abundant (KtW) species the infecting phage results in collapse of a constant number K of other bacterial species. This version of the model is reminiscent of the Bak-Sneppen model of species co-evolution 31 . We tested this model for K = 1 and K = 25 (out of N = 500). In the first case we observed no impact on diversity, while in the second case the diversity saturated at lower values of σΔ t. All together we can conclude that the diversity profile shown in Fig. 1B remains qualitatively (and sometimes even quantitatively) unaffected by additional interactions between microbial species or more interconnected PBINs.
According to our results the principal determinant of the ecosystem diversity D is the width σ of the distribution of logarithmic growth rates of individual bacterial strains or species. This difference is amplified during the average time Δ t between population collapses. Thus the overall frequency of collapses is a very important parameter with more frequent collapses counter-intuitively resulting in more diverse ecosystems. That is because in our scenario frequent collapses weaken the effect of competitive exclusion ultimately driving the diversity down to no more than single species per rate-limiting nutrient. Larger magnitude of collapses also promotes higher diversity but its impact increases only weakly (logarithmically) with the collapse ratio γ.
It is instructive to compare the determinants of microbial diversity in the static, steady state KtW model and in our more dynamic, collapse-driven variant. In the static KtW model 2,5 the steady state population size of each of the bacteria B * = δ/βη is determined exclusively by parameters of the phage to which it is susceptible: its burst size (β), death (or dispersal and dilution) rate (δ), and its infection rate (η) at a density equal the bacterial carrying capacity. This steady state population of a phage-controlled bacterium is usually much lower than the carrying capacity of the environment: B * ≪ 1. Thus a large number of bacteria each susceptible to its unique phage predator can coexists with each other 5 . Higher diversity can subsequently be achieved by carefully adding pairs of bacteria and phages, latter possibly supplemented by their epigenetic variants 32 , each consuming a small fraction of the carrying capacity 5,9 . Substantial diversity is found to be fragile to new invaders, in the form of bacteria that grow faster than resident ones or phages that prey on several bacteria at once 9 .
In contrast to this the diversity in our model is determined by both statistics of collapses as well as the spread of growth rates of resident bacterial species. In case of mild or infrequent collapses and large disparity in bacterial growth rates competitive exclusion principle is restored within our model as it then predicts an ecosystem dominated by just one fastest growing bacterium. When collapses are frequent (short Δ t) and severe (large |log γ|), while growth rates of individual bacterial strains or species are close to each other (small σ), Eq. 4 predicts high diversity of co-existing bacterial species. This prediction is robust with respect to exact causes of collapses, including the relatively frequent 16 invasion of phages that are capable of infecting several bacterial species.
Overall, the falsifiable (and counter-intuitive) prediction of the collapse-driven "Kill-the-Winner" model differentiating it from its stationary counterpart, is that by increasing frequency (and to a smaller extent severity) of collapses one could support higher diversity of microorganisms. In real world ecosystems, this particular aspect of density dependent selection is to be supplemented by more classical engines of microbial diversity 2,5,9,32 and in addition be modulated in case spatial heterogeneity 33,34 reduce the assumed collapse sizes.