Diversity waves in collapsedriven population dynamics
Abstract
Abstract
Populations of species in ecosystems are often constrained by availability of resources within their environment. In effect this means that a growth of one population, needs to be balanced by comparable reduction in populations of others. In neutral models of biodiversity all populations are assumed to change incrementally due to stochastic births and deaths of individuals. Here we propose and model another redistribution mechanism driven by abrupt and severe collapses of the entire population of a single species freeing up resources for the remaining ones. This mechanism may be relevant e.g. for communities of bacteria, with strainspecific collapses caused e.g. by invading bacteriophages, or for other ecosystems where infectious diseases play an important role.
The emergent dynamics of our system is cyclic “diversity waves” triggered by collapses of globally dominating populations. The population diversity peaks at the beginning of each wave and exponentially decreases afterwards. Species abundances are characterized by a bimodal timeaggregated distribution with the lower peak formed by populations of recently collapsed or newly introduced species, while the upper peak  species that has not yet collapsed in the current wave. In most waves both upper and lower peaks are composed of several smaller peaks. This selforganized hierarchical peak structure has a longterm memory transmitted across several waves. It gives rise to a scalefree tail of the timeaggregated population distribution with a universal exponent of 1.7. We show that diversity wave dynamics is robust with respect to variations in the rules of our model such as diffusion between multiple environments, speciesspecific growth and extinction rates, and bethedging strategies.
pacs:
Author Summary
The rate of unlimited exponential growth is traditionally used to quantify fitness of species or success of organizations in biological and economic context respectively. However, even modest population growth quickly saturates any environment. Subsequent resource redistribution between the surviving populations is assumed to be driven by incremental changes due to stochastic births and deaths of individuals. Here we propose and model another redistribution mechanism driven by sudden and severe collapses of entire populations freeing up resources for the growth of others. The emergent property of this type of dynamics are cyclic “diversity waves” each triggered by a collapse of the globally dominating population. Gradual extinctions of species within the current wave results in a scalefree timeaggregated distribution of populations of most abundant species. Our study offers insights to population dynamics of microbial communities with local collapses caused e.g. by invading bacteriophages. It also provides a simplified dynamical description of market shares of companies competing in an economic sector with frequent rate of bankruptcy.
Introduction
Mathematical description of many processes in biology and economics or finance assumes longterm exponential growth fisher1930 (); kelly () yet no reallife environment allows growth to continue forever malthus (); verhulst (). In biology any growing population eventually ends ups saturating the carrying capacity of its environment determined e.g. by nutrient availability. The same is true for economies where finite pool of new customers and/or natural resources inevitably sets a limit on growth of companies. Population dynamics in saturated environments is routinely described by neutral “community drift” models hubbell2001 (); woodcock2007 () sometimes with addition of deterministic differences in efficiency of utilizing resources sloan2005 ().
Here we introduce and model the saturatedstate dynamics of populations exposed to episodic random collapses. All species are assumed to share the same environment that ultimately sets the limit to their exponential growth. Examples of such systems include local populations of either a single or multiple biological species competing for the same nutrient, companies competing to increase their market shares among a limited set of customers, etc. Specific case studies can be drawn from microbial ecology where susceptible bacteria are routinely decimated by exposure to bacteriophages (see e.g. middleboe2001 (); haerter2014 () and references therein), or paleontological record, where entire taxonomic orders can be wiped out by sudden extinctions happening at a rate independent of order size bornholdt2009 ().
Model
Population growth of a single exponentially replicating species in a growthlimiting environment is traditionally described by Verhulst’s verhulst () or logistic equation , where the carrying capacity of the environment without loss of generality can be set to . In this paper we consider the collective dynamics of multiple populations competing for the same ratelimiting resource:

Local populations are subject to episodic random collapses or extinctions. The probability of an extinction is assumed to be independent of the population size. Immediately after each collapse the freedup resources lead to the transient exponential population growth of all remaining populations . The growth stops once the global population reaches the carrying capacity .

The environment is periodically reseeded with new species starting from the same small population size (measured in units of environment’s carrying capacity).
We initially assume that the growth rates and collapse probabilities of all species are equal to each other. We also disregard the neutral drift hubbell2001 () in sizes of individual populations during the time between subsequent collapses. All these assumptions will be relaxed, simulated, and discussed later in the paper (see Supplementary Materials S1 Text, S1S7 Figures). The number of species in the steady state of the model is determined by the competition between the constant rate of introduction of new species and the overall rate of extinctions in the environment that is proportional to the number of species. To simplify our modeling we will consider a closely related variant of the model in which the number of species is kept strictly constant. In this case each extinction event is immediately followed by the introduction of a brand new species. We have verified that the two version of our model have very similar steady state properties. The dynamics of the fixed model is described by
(1) 
where is the random variable which is equal to zero for surviving populations and has a large positive value for populations undergoing an extinction/collapse.
To speed up our simulations we do not continuously calculate Eq. (1) since most of the time the carrying capacity of the environment is saturated when local populations do not change. Instead we simulate the model at discrete time steps marked by extinction events. At every time step a randomly selected local population goes extinct and a brand new species with population is added to the environment. We then instantaneously bring the system to its the carrying capacity by multiplying all populations by the same factor.
Results
In spite of its simplified description of the ecosystem disregarding pairwise interactions between species our model has a rich population dynamics. Figure 1A shows the timecourses of populations in a system with species and .
At certain times the carrying capacity of the environment is nearly exhausted by just one dominant species with visible as dark red stripes in Figure 1A. When such dominant species goes extinct a large fraction of the resources suddenly becomes available and consequently all other populations increase by a large ratio . This marks the end of one and the start of another diversity wave that initially is dominated by a large number of species with similar population sizes. In the course of this new wave these species are eliminated onebyone by random extinctions until only one dominant species is left standing. Its collapse ends the current and starts the new wave. In Figure 1A one can clearly distinguish about 5 such waves terminated by the extinctions of dominant species #5, 15, 6, 19, and 16 correspondingly.
Figure 1B shows the timecourse of just one local population of the species #6 on a logarithmic scale. Between time steps 100 and 150 the population of this species nearly exhausts the carrying capacity of the environment. Its local extinction at the time step 154 ended the third diversity wave and started the fourth one. Note somewhat erratic yet distinctly exponential growth of this species happening on the slow timescale set by the frequency of extinctions. This growth should not be confused with exponential repopulation of recently collapsed environments that happens much faster (a small fraction of one time step).
Figure 2 follows the population diversity (grey shaded area) defined as as a function of time in a system of size . In general can vary between when one abundant species dominates the environment and when all species are equally abundant. The diversity is inversely proportional to the largest population . The diversity waves (purple dashed arrows in Figure 2) are initiated when a dominating species collapses. As a consequence, at the start of each wave the diversity abruptly increases from to a substantial fraction of the maximal possible diversity . After this initial burst triggered by the global redistribution of biomass, the diversity exponentially declines as (the dotdashed line in Fig.2), driven by ongoing extinctions of individual populations. The typical duration, of a diversity wave is equal to the time required for all of the species present at the start of the wave to go extinct onebyone. Thus it is determined by or
(2) 
Figure 3 shows the timeaggregated distribution of population sizes for and (the grey shaded area). This logarithmicallybinned distribution defined by were collected over 20 million individual collapses (timesteps in our model). Thus, a timeaggregated distribution is rather different from a typical “snapshot” of the system at a particular moment in time characterized by between and of highly abundant species in the current diversity wave. The timeaggregated distribution is bimodal with clearly separable peaks corresponding to two population subgroups. The upper peak consists of the species that have not yet been eliminated in the current wave.
To better understand the dynamics of the system in Figure 3 we also show the distribution of populations sizes at the very end of each diversity wave (green line) and at the beginning of the next wave (red line). That is to say, for the green curve we take a snapshot of all populations immediately after the dominant species with was eliminated, but before the available biomass was redistributed among all species. At those special moments, happening only once every time steps, most population sizes are between and while a small fraction reaches all the way up to . During the rapid growth phase immediately after our snapshot was taken, all populations grow by the same factor thereby moving all of them to the upper peak of the timeaggregated distribution thereby starting the new wave. The red curve corresponds to the snapshot of all populations immediately after this rescaling took place. It shows that at the very beginning of the new wave local populations are broadly distributed between and with a peak around .
Figure 4A shows timeaggregated distributions of population sizes for and different values of ranging between ad while Figure 4B shows timeaggregated distributions with and for a wide range of . One can see that for , the tail of the distribution for most abundant populations between and is well fitted by a power law (dashed line in Figure 4B) corresponding to the power law distribution of population sizes on the linear scale . Overall Figures 4A,B demonstrate that the exponent for different values of and is remarkably universal. Indeed, for a range of parameters where the upper and the lower peaks of are clearly separated, approaches a universal value .
An insight into the origins of the scalefree tail of the distribution of population sizes is gained by considering a simplified version of our model in which at the start of each wave the populations of all species are artificially set to be equal to each other resulting in the maximal diversity. We further assume that . The passage of time elapsed since the beginning of the current wave, leads to a decrease in the number of surviving species , which all have the same population size jointly filling up the carrying capacity of the environment. Above we ignore a negligible fraction () of the total population of the lower peak species. The timeaggregated probability for a species to have a population size is naturally given by and thus
The exponent predicted by this simplified model is realized in our actual model for moderately high , whereas smaller values of give rise to a different universal exponent . The decrease of the exponent from 2 to 1.7 in our original model is the result of unequal population sizes at the beginning of a new wave. In fact, we verified numerically that is recovered if at the start of each wave one equilibrates all species abundances to . The first section of the S1 Text in supplementary materials provides additional details on how the reduced population diversity at the start of population waves affects the exponent .
Two panels in Figure 5 illustrate the difference between the simplified (panel A) and the real (panel B) models. In both versions of the model the average jump in the logarithm of surviving populations grows exponentially with time elapsed since the start of the current wave: . However, unlike the simplified model, the population distribution in our real model has a rich hierarchical structure with multiple subpeaks in some waves (color bands in Figure 5B). Remarkably this multimodal distribution reappears in subsequent waves, implying that the memory about the hierarchical structure in the upper part of the distribution is transmitted to emerging populations in the lower part with sizes starting at . At the start of the next wave these same populations will move to the upper part of the distribution thereby transmitting the history across waves. Colors of symbols in Figure 5 illustrate the origin of multiple peaks. Indeed, populations in each of these peaks were born during the previous wave under similar conditions (the number of substantial populations) as described in the caption. Thus, the broadening of timeaggregated population distribution in our model compared to its historyfree version is a simple manifestation of a complex interplay between ”upstairs” and ”downstairs” subpopulations transmitting memory across waves.
The population distributions in both upper and lower peaks in our model are described by the same exponent . This similarity reflects the fact that individual populations in the lower peak are exposed to the same multiplicative growth as the ones in the upper peak. Finally, the intermediate region of the distribution connecting two peaks is shaped by rescaling of all populations in the lower peak as they are moved up at the beginning of a new diversity wave. When peaks are well separated (as e.g. for low values of in Figure 4) the slope of the logarithmic distribution in this region has almost exactly the same value and the opposite sign to the slopes in both the upper and the lower peaks.
Discussion
In this paper we explore the population dynamics in saturated environments driven exclusively by nearcomplete collapses of subpopulations of competing species. This type of dynamics strongly contrasts the gradual changes implied in for example the “community drift” neutral models hubbell2001 () in ecology, or for the most part incremental random walks of stock values of individual companies. Conversely, collapsedriven dynamics represents a sudden and usually large change of system composition. In ecology such collapses may be caused e.g. by invading predators or diseases , whereas in the economy, companies of any size routinely go bankrupt e.g. through excessive debt amplifying the effects of external perturbations.
First, let us consider biological systems. One of the predictions of our model is a multimodal logarithmic distribution of population sizes. Indeed, while the timeaggregated distribution is bimodal with distinct upper and lower peaks, populations within any given diversity wave cluster together in several smaller peaks persisting over several waves (color stripes in Figure 5B). This overall finding is supported by a growing body of literature gray2005 (); dornelas2008 (); vergnon2012 (); matthews2014 () where multimodal Species Abundance Distributions (SAD) in real ecosystems were reported for plants, birds, arthropods matthews2014 (), marine organisms including single cells, corals dornelas2008 (), nematodes, fishes, entire seafloor communities gray2005 (), and even extinct brachiopods olszewski2004 (). Like in our model, the empirical SADs range over many orders of magnitude with a noticeable depletion (or several depletions) at intermediate scales. The magnitude of this dip is usually somewhat less than predicted by our basic model but is consistent with several of its variants described below. This includes the model variant #1 inspired by the neutral theory of biodiversity hubbell2001 () thought to apply to a variety of ecosystems including microbial communities sloan2005 (); woodcock2007 () (see S1 Figure in supplementary materials).
Needless to say, our model is not unique in generating multimodal distributions (see e.g. vergnon2012 () for other examples). Conversely, some of the variants of our model give rise to interesting population dynamics including diversity waves even without any depletion in the middle of the logbinned SAD. We argue that a more reliable characterization of underlying dynamical processes can be obtained from timeseries data. First, all systems capable of diversity waves are described by rapid large changes in populations of individual species. Such sudden, populationscale shifts can occur e.g. due to introduced diseases or the disappearance of keystone species paine1969 (); cohn1998 () thereby changing the composition of the entire foodweb. On the microbial scale, sudden invasion of a new bacteriophage may lead to multiple orders of magnitude reduction in the population of susceptible bacteria levin1977 (); middleboe2001 (), potentially leading to their complete local extinction haerter2014 (). Phagedriven collapses do not spare bacteria with large populations and may even be biased towards such as postulated in the KilltheWinner (KtW) hypothesis thingstad1997 (). The magnitude and characteristic timescales of population changes in microbial ecosystems is still being actively discussed in the literature. While Ref. campbell2011 () reports that over half of all bacterial species in marine environments varied between abundant and rare over a threeyear period, other studies rodriguezbrito2010 () found more modest variability at the level of species or genera over weeks to months period. However, everyone seem to agree on dramatic and rapid (often on the scale of days castberg2001 ()) population shifts at the level of individual bacterial strains middleboe2001 (); middelboe2009 (); rodriguezbrito2010 () caused by phage predation castberg2001 (). Except for interchangeable gene cassettes (metagenomic islands) responsible for either phage recognition cites or alternatively resistance to phages rodriguezvalera2009 () these strains routinely have very similar genomes and thus may have near identical growth rates. Hence, they are capable of coexistence in the saturated state implicitly assumed in our model.
Extinctions and collapses in our model are assumed to be caused exclusively by exogenous effects such as natural catastrophes or predation by external species not sharing the carrying capacity of our environment. Reallife ecosystems can also collapse due to endogenous effects, i. e. internal interactions between species. Such intrinsic collapse mechanisms were the sole focus of earlier models by us and others (see e.g. bak1993 (); sole1996 (); haerter2014 ()).
On vastly longer, geological timescales, the collapsedriven dynamics of our model resembles that of species extinctions and subsequent radiations in the paleontological record raup (); vanvalen (). One example is the recolonization by mammals of a number of ecological niches vacated by dinosaurs after the endCretaceous mass extinction thought to be preceded by a gradual depletion of diversity among dinosaurs who were finally wiped out by a singular catastrophic event sloan1986 (). Interestingly, the extinction rate of taxonomic orders appears to be independent of their size quantified by the number of genera they contains bornholdt2009 (), which is also one of the assumptions of our basic model of collapsedriven dynamics.
The second area of applications of our model is to describe company dynamics in economics. The size or the market share of a publicly traded company reflected in its stock price is well approximated by a random walk with (usually) small incremental steps bachelier1900 (). However, as in the case of ecosystems, this smooth and gradual drift does not account for dramatic rapid changes such as bankruptcies or market crashes. In case of companies the main driver of sudden changes is their debt fisher1933 (). When the intrinsic value of a company is much smaller than its debt, even small changes in its economical environment can make it insolvent not sparing even the biggest companies from bankruptcies (think of Enron and Lehman Brothers). Empirically, the yeartoyear volatility of company’s market share varies with its size , econo_hurst ().
Abundance distributions in our original model and many of its variants are characterized by a powerlaw tail with an exponent close to . This is in approximate agreement with the empirical data on abundance distributions of bacteria in soil samples bacterial_zipf (), marine viruses (phages) viral_zipf ().
In the economics literature, the distributions of company sizes filiasi2014 (), as well as those of wealth of individuals klass2004 () are known to have similar scalefree tails. Recent data for companies filiasi2014 () and personal wealth klass2004 () suggest tail of the former distribution and tails of the latter one. Traditionally, scalefree tails in these distributions were explained by either stochastic multiplicative processes pushed down against the lower wall (the minimal population or company size, or welfare support for low income individuals) levy1996 (); levy1997 (); sornette (), by variants of richgetricher dynamics simon1955 (), or in terms of SelfOrganized Criticality bak1987 (); bak1993 (). The emphasis of the latter type of models on large systemwide events (avalanches bak1987 (); bak1993 () or collapses newman1996 ()) and on separation of timescales is similar in spirit to the collapsedriven dynamics in our models. A potentially important socioeconomic implication of our model is that during each wave contingent subpeaks in the “upstairs” part of the distribution are imprinted on the “downstairs” part and thereby can be repeated in the new wave following the “revolution”.
Needless to say, our models were simplified in order to compare and contrast the collapsedriven dynamics to other mathematical descriptions of competition in saturated environments. The S1 Text in supplementary materials describes several variants of our basic model that in addition to population collapses include the following elements:

“Neutral drift model” assumes changes of population sizes during time intervals between collapses as described in Ref. hubbell2001 (). In this model in addition to collapses a population of size randomly drifts up and down . The resulting diversity waves and timeaggregated distributions can be found in the supplementary S1 Figure.

“Exponential fluctuations model” is another variant of the neutral scenario where the population sizes between collapses undergo slow multiplicative adjustments restricted by the overall carrying capacity of the environment. Details and the resulting timeaggregated distribution can be found in the supplementary S2 Figure.

“Interconnected environments model” is another neutral variant of our basic rules in which spatially separated subpopulations of the same species are competing with each other for the same nutrient. Subpopulations are connected by the diffusion, that is much slower than the diffusion of shared nutrient. In this model a collapsed subpopulation is replenished by a small number of individuals diffusing from other environments, see the supplementary S3 Figure.

”KilltheWinner (KtW) model” where collapse probability systematically increase with the population size as suggested by the studies of phagebacteria ecosystems thingstad1997 (). In this particular case the diversity dynamics and the scalefree tail of the population distribution becomes sensitive to the extent that the large populations are disfavoured by collapse. When the collapse probability is proportional to population size, one obtain a flat distribution where numbers of species in each decade of population sizes are equal to each other, see the supplementary S4 Figure.

”KilltheLooser (KtL) model”, where collapse probability systematically decreases with the population size as as suggested by the empirical studies of company dynamics econo_hurst (). As seen in the supplementary S5 Figure the diversity dynamics and the scalefree tail of the population distribution are both remarkably robust with respect to introduction of sizedependent collapse rate.

”Fitness model” in which each of the species has its own growth rate during rapid repopulation phase and its own collapse probability . The supplementary S6 Figure show that the overall shape of the timeaggregated distribution is similar to that in our basic model, whereas its lower panel illustrate the interplay between population size and the the two variables that define the species’ fitness.

“Resilience model” as a variant of the above fitness scenario, in which collapsing species do not necessarily go into extinction. Instead, each species is assigned its own “survivor ratio” defined by the population drop following a collapse: . As in the previous variant each of the species is also characterized by its own growth rate . The supplementary S7 Figure shows that for intermediate populations the timeaggregated distribution is described by a power law scaling. Compared to the basic model it has a broader scaling regime and larger likelihood to have most of the “biomass” collected in one species.
Captions to supplementary S1S7 Figures provide more detailed description of model dynamics in each of these cases. Overall, the ’ simulations of the variants of our basic model described above preserve the general patterns of collapsedriven dynamics such as diversity wave dynamics, and a broad timeaggregated distribution of population sizes with scalefree tail for the most abundant species.
The classic definition of the fitness of a species in terms of its longterm exponential growth rate fisher1958 () is clearly inappropriate for our model. Indeed, the longterm growth rate of each of our species is exactly zero. One must keep in mind though that fitness is a very flexible term and has been defined in in a variety of ways orr2009 () reflecting (among other things) different timescales of growth and evolution bak1993 (), and relative emphasis on population dynamics vs. risk minimization bergstrom2004 (). An appropriate way to quantify species’ success in a steady state system like ours is in terms of their timeaveraged population size .
In the last two variants of our basic model we add fitnessrelated parameters to each of the species: its shortterm exponential growth rate (model 6 and 7), its propensity to large population collapses quantified by their frequency (model 6), and the severity of collapses quantified by surviving fraction of the population (model 7). Fig. 6 shows the average population size as function of and in the model 6. As one can naively expects species with larger shortterm growth rates and larger surviving ratios also tend to have substantially larger populations (the red area in the upper right corner of Fig. 6).
While in our models the probability and severity of collapses are assumed to be independent of growth rate in reality they are often oppositely correlated. Indeed, in biology much as in economics/finance fast growth usually comes at the cost of fragility and exposure to downturns forcing species to optimize tradeoffs.
Some environmental conditions favor the fast growth even at the cost of a higher risk of collapse, whereas other could call for sacrificing growth to minimize probability or severity of collapses. Species in frequently collapsing environments considered in our study are known to employ bethedging strategies kelly (); zhang (); bergstrom2004 (); ms_kelly2015 () to optimize their longterm growth rate. This is obtained by splitting their populations into “growthloving” and “riskaverse” phenotypes bergstrom2004 (); kussell2005 (); ms_kelly2015 (). One example of this type of bethedging is provided by persistor subpopulations of some bacterial species consisting of of the overall population baleban2004 (); gerdes2013 () increasing to as much as in saturated conditions (S. Semsey, private communications). A bethedging strategy with persistor subpopulation somewhat reduces the overall growth rate (only by 1%) while dramatically reducing the severity of collapses caused by antibiotics. From Fig. 6 one infers that this is indeed a good tradeoff.
In this study we presented a general modeling framework for systems driven by redistribution of ratelimiting resources freed up by episodic catastrophes. In spite of their simplicity the population dynamics in such systems happens on at least four hierarchical timescales. At the shortest timescale the populations grow exponentially repopulating resources vacated during a catastrophic extinction event. This exponential growth results in a steady state at which the system is poised exactly at the carrying capacity of the environment. At even longer timescales the system is described in terms of diversity waves that are the main focus of this study. These waves are an emergent dynamical property of the system in which population of surviving species grow exponentially while species diversity decays. Remarkably the information about the “upstairs”and “downstairs” population peaks survives the “revolution” at the end of each wave. This memory gives rise to the final, longest timescale in our system correlating several consecutive waves. All of this complexity is already present in our basic oneparameter model. In spite of its simplicity this model and its variants provide the foundation for future studies of collapsedriven dynamics in ecosystems, market economies, and social structures.
References:
Supplementary Materials
FokkerPlanck equation for the basic model
Let’s consider a version of our model with a very low value of to ensure the complete separation between the upper and lower peaks. The multiplicative dynamics of surviving populations in the upper peak is described by the following elementary step:
. Here is the population at time step of the species that went extinct at this time step. It can be also easily integrated for all times since the beginning of the current wave at time step :
(S1) 
Indeed, is the total initial (at time step ) populations of all species that went extinct by the time step . Hence  is the total initial population of all surviving species used to normalize their initial populations to give their populations at the time step (population ratios of surviving species are preserved in our basic model). Taking the logarithm of both sides of Eq. S1 and approximating , which holds as long as the system is still far away from the end of the wave (, one gets:
(S2) 
The stochastic dynamics within a single wave can thus be described by the following equation:
(S3) 
The exponent of the population size distribution in our model is determined by the balance between the noisy multiplicative population dynamics and the exponential loss of surviving species due to collapses. It can be approximated by the following FokkerPlanklike equation:
(S4) 
Here is the timedependent population abundance distribution in the upper peak,  the loss term due to population collapses,  the logarithmic drift velocity and  is the logarithmic dispersion, which is totally absent in the simplified model where all populations start at the same size.
The Eqs.S2S3 allow us to derive the parameters of the FokkerPlank equation in terms of the distribution of population sizes at the start of the wave. Indeed, early in the wave one has , and . Note an unusual connection between the population diversity at the start of a wave , and the diffusion coefficient in the FokkerPlank equation.
The stationary solution for the timeaggregate distribution has an exponential tail . It corresponds to the power law tail of the species population distribution . The exponent is defined by one of the two solutions to the quadratic equation
(S5) 
In the version of the model, where all populations at the start of the wave are equal to each other, sizes of surviving populations increase deterministically as (see main text for the derivation) and thus have zero dispersion: . Hence in this simplified version the exponent is determined by balancing only the first two terms of this equation.
We have numerically verified that the decrease of the exponent from in the simplified model down to in our original model is driven entirely by noise (unequal population sizes) resulting in a finite value of . A nonzero value of in the Eq. S5 results in . For example, if the populations at the start of each wave had a Poisson distribution so that , the exponent would have been defined by the solution of the golden mean equation . While currently we have no firstprinciples argument allowing us to derive the value of in our basic model, the result from a Poisson distribution is not too far from the empirically observed exponent
Model variants
To test the robustness of our basic model with respect to rule changes we considered the following seven variants:

“Neutral drift model”.
This variant extends our basic model by adding to our standard model the random neutral drift of population sizes (Hubbell SP (2001) The unified neutral theory of biodiversity and biogeography (MPB32), Princeton University Press) between subsequent collapse events. To simulate this random drift, at every time step the population of each species changes up or down as prescribed by dispersion of binominal distribution: , where is the parameter quantifying the magnitude of fluctuations proportional to the inverse of the total population size and the square of the birth/death rate. After drift changes were applied to all populations we rescale them back to their carrying capacity . This is followed by a collapse event as in the standard model. Fig. S1 illustrates typical time courses of the diversity and timeaggregated species abundance distributions in this model variant for three values of and compares them to our basic model.

“Exponential fluctuations model”, where populations are exposed to random exponential shifts between successive collapse events. In this version of the model between successive collapse events all populations exponentially shift thus adjusting the way carrying capacity is divided between the. This model is similar to the ”Neutral random drift” model #1 above except that changes are proportional to and not to . The population of each species is characterized by its own exponential rate given by , where quantifies the overall rate of the redistribution, while  a random number uniformly distributed between 0 and 1  represents species and time specific shifts. We assume that differences in growth rates between species are not constant but instead fluctuate on the timescale when a single species collapses. Hence, after each collapse event we reset the exponential growth we randomly reset the rates for all species (and not just of the collapsed one). In this version of the model we also take into account the stochastic nature of the time interval between two successive collapse events, which is randomly chosen from the exponential distribution with mean value equal to 1. Thus between two successive collapse events each species population changes as and subsequently rescaled to the carrying capacity of the environment: . As in our basic model, at every time step one randomly selected population collapses and is reset to while all other populations are rescaled to fill up the carrying capacity: . Fig. S2 shows the simulations of this model for several value of compared to our basic () model.

“Interconnected environments model” with diffusion. In this version of the model there is a single species distributed between local environments. As in our basic model at every time step one local population is selected for collapse and reset to after which the populations are normalized back to the carrying capacity. The diffusion takes a fraction of the total population of and distributes it equally between all local environments: . Note, that here we implicitly assume that populations in all of these environments share the same carrying capacity. This is the case when diffusion rate of the ratelimiting nutrient is much faster then that of populations themselves. The main difference of this model from earlier variants is that populations in the lower peak with grow approximately linearly in time (as opposed to exponentially in other versions of the model). The rate of this linear growth is the same for all species and is equal to . The exponential growth is restored for populations that are larger than average. Fig. S3 shows the timeaggregated distribution of local population sizes in this model variant.

“KilltheWinner (KtW) model”. For bacterial populations the direction of the trend (if any) of collapse probability with population size is currently unknown. In fact one can plausible make a case for increasing of the probability of collapse with population size due to larger populations making easier to find and overall and more attractive targets for phages. In microbiology preferential targeting of large bacterial populations by virulent phages is known as ”KilltheWinner” (KtW) hypothesis (Thingstad TF and Lignell R (1997) Aquatic Microbial Ecology 13:1927). Here we simulate the version of our basic model where the collapse probability systematically increases with population size. At each time step we select a random population to collapse with probability . As before the collapsing population is reset to and all populations are subsequently grown with equal exponential rates to complete saturation: . Fig. S4 examines timeaggregated population distributions in KtV model variant for different values of . Whereas small and moderate preserve diversity wave dynamics, the version does not exhibit diversity waves and predicts a population size distribution distributions, , or (equal number of species in each decade of population sizes).

“Killthelooser (KtL) model” in which a collapse probability declines with population size in a powerlaw fashion. At each time step one select a random population to collapse with probability where is the current population size of species . In economics this corresponds to an intuitively plausible notion that larger companies are less likely to go bankrupt than smaller ones. Empirically, this trend is described by a power law with the exponent 0.2 (Nunes Amaral LA, Buldyrev SV, Havlin S, Leschhorn H, Maass P, Salinger MA, et al. (1997) Journal de Physique I. 7: 621633.) Notice the emergence of the lower peak distribution distribution above that is much more narrow than in our standard model. This makes sense as in the course of each wave small populations tend to collapse over and over. These repeated collapses don’t drive other populations up by much and thus their only consequence is clustering of small populations close to the very bottom of the lower peak distribution at and above . When the dominant upper peak population finally collapses all small populations are rescaled up to form a narrow distribution around . This is very similar to our simplified memoryfree model described in the main text and shown in Fig. 4A. The new wave starts with very high diversity which is subsequently reduced with time as . Here we ignore a relatively small fold decline in collapse frequency over the range of population sizes between and . Within the same approximation each surviving population grows as . The timeaveraged distribution of populations thereby approaches the scaling regime described by:
(S6) In reality the scaling exponent of the tail is around . It is the same as in our standard model but for a different reason. Indeed, taking into account that lifetime of a population before collapse scales as one gets .

“Fitness model” with heterogeneous, speciesspecific growth rates and extinction probabilities. Each species is assigned a growth rate used when it repopulates the freedup carrying capacity of the environment. It also has its own extinction probability . Both and are logarithmically distributed in the interval between 0.1 and 1. That is to say their are uniformly distributed between and . At each time step we select one of populations, with probability , this species goes extinct. It is immediately replaced by a new species with the population , new growth rate , and extinction probability . Subsequently all of the populations are rescaled proportional to their growth rates
to fill up the carrying capacity of the environment . The upper panel in Fig. S6 shows that the timeaggregated population distribution in this model preserves its powerlaw tail, whereas lower panel illustrates that in order for a species to reach substantial population size its fitness parameters need to be particularly favorable. Indeed, populations larger than tend to have smaller than average extinction probabilities , and larger than average growth rates .

“Resilience model” where heterogeneous, speciesspecific growth rates and survival ratios (population reduction following a collapse) are competing with each other. Each species is assigned a growth rate and collapse size , both logarithmically distributed (uniform distribution of the logarithm of the variable). At each time step we select one of populations, and collapse its population . Note that unlike in previous versions we scale down the population proportional to its size and not proportional to the carrying capacity of the environment. This reflects a different interpretation away from our basic model, where a collapse represents not an extinction of the species followed by the appearance of a new species at a fixed (very small) population size. In the new version of the model a collapse represents a sudden but proportionate reduction of a population e.g. due to species’ phenotypic or genotypic bethedging. In this version of the model an extinction happens only if a very low population is reached, i.e. when . If this lower bound is reached the old species goes extinct and a new species with the initial population is introduced. The new species is assigned new random values of and . As in the previous model following each collapse all populations are rescaled proportional to their growth rates to fill up the carrying capacity of the environment: .
References
 Fisher RA (1930) The Genetical Theory of Natural selection. Clarendon, Oxford
 Kelly JL (1956) A new interpretation of information rate. Bell System Technical Journal 35:917926
 Malthus TR (1798) An Essay on the Principle of Population. J. Johnson, in St. Paul’s Churchyard, London
 Verhulst PF (1838) Notice sur la loi la population poursuit dans son accroissement. Correspondance Mathematique er Physique, 10:113121
 Hubbell SP (2001) The unified neutral theory of biodiversity and biogeography (MPB32), Princeton University Press
 Woodcock S, Gast VC (2007) Neutral assembly of bacterial communities. FEMS Microbiol Ecol 62: 171180
 Sloan WT, Lunn M, Woodcock S (2006) Quantifying the roles of immigration and chance in shaping prokaryote community structure. Environmental Microbiology, 8: 732740
 Middelboe M, Hagstrom A, Blackburn N, Sinn B, Fischer U, et al. (2001) Effects of Bacteriophages on the Population Dynamics of Four Strains of Pelagic Marine Bacteria. Microb Ecol 42: 395406;
 Haerter JO, Mitarai N, Sneppen K (2014) Phage and bacteria support mutual diversity in a narrowing staircase of coexistence. The ISME journal. The ISME Journal 8: 23172326.
 Bornholdt S, Sneppen K, and Westphal H (2009) Longevity of orders is related to the longevity of their constituent genera rather than genus richness. Theory in Biosciences 128: 7583
 Gray JS, Bjoergesaeter A, Ugland KI (2005) The impact of rare species on natural assemblages. Journal of Animal Ecology 74: 11311139.
 Dornelas M, Connolly SR (2008) Multiple modes in a coral species abundance distribution. Ecology Letters 11: 10081016.
 Vergnon R, van Nes EH, Scheffer M (2012) Emergent neutrality leads to multimodal species abundance distributions. Nature Communications 3: 663.
 Matthews TJ, Borges PA, Whittaker RJ (2014) Multimodal species abundance distributions: a deconstruction approach reveals the processes behind the pattern. Oikos 123: 533544.
 Olszewski TD, Erwin DH (2004) Dynamic response of Permian brachiopod communities to longterm environmental change. Nature 428: 738741.
 Paine RT (1969) A note on trophic complexity and community stability. American naturalist 103:9193
 Cohn JP (1998) Understanding Sea Otters. BioScience 48: 151155
 Levin BR, Stewart FM, and Chao L (1977) ResourceLimited Growth, Competition, and Predation: A Model and Experimental Studies with Bacteria and Bacteriophage., The American Naturalist 111: 324
 Thingstad TF and Lignell R (1997) Theoretical models for the control of bacterial growth rate, abundance, diversity and carbon demand. Aquatic Microbial Ecology 13:1927
 Campbell BJ, Yu L, Heidelberg JF, Kirchman DL (2011) Activity of abundant and rare bacteria in a coastal ocean. Proceedings of the National Academy of Sciences 108: 1277612781.
 RodriguezBrito B, Li L, Wegley L, Furlan M, Angly F, et al. (2010) Viral and microbial community dynamics in four aquatic environments. The ISME journal 4: 739751.
 Castberg T, Larsen A, Sandaa RA, Brussaard C, Egge JK, et al. (2001) Microbial population dynamics and diversity during a bloom of the marine coccolithophorid Emiliania huxleyi (Haptophyta). Marine Ecology Progress Series 221. Available: http://depot.knaw.nl/10987/
 Middelboe M, Holmfeldt K, Riemann L, Nybroe O, Haaber J (2009) Bacteriophages drive strain diversification in a marine Flavobacterium: implications for phage resistance and physiological properties. Environmental microbiology 11: 19711982.
 RodriguezValera F, MartinCuadrado AB, RodriguezBrito B, Pai L, Thingstad TF, et al. (2009) Explaining microbial population genomics through phage predation. Nature Reviews Microbiology 7: 828836.
 Bak P and Sneppen K (1993) Punctuated Equilibrium and Criticality in a simple model of Evolution. Phys. Rev. Lett. 71: 40834086.
 Sole RV and Manrubia SC (1996) Extinction and selforganized criticality in a model of largescale evolution. Phys. Rev. E 54, R4245
 Raup DM. (1992) Extinction: bad genes or bad luck? WW Norton and Company.
 L. Van Valen. A new evolutionary law. Evolutionary Theory, 1:1, (1973).
 Sloan RE, Rigby JK, Van Valen LM, Gabriel D (1986) Gradual dinosaur extinction and simultaneous ungulate radiation in the Hell Creek Formation. Science 232: 629633.
 Bachelier L. (1900) Theorie de la speculation. GauthierVillars, 70 pp.
 Fisher I (1933) The debtdeflation theory of great depressions. Econometrica 1:337357.
 Nunes Amaral LA, Buldyrev SV, Havlin S, Leschhorn H, Maass P, Salinger MA, et al. (1997) Scaling Behavior in Economics: I. Empirical Results for Company Growth. Journal de Physique I. 7: 621633.
 Gans J, Wolinsky M and Dunbar J (2005) Computational Improvements Reveal Great Bacterial Diversity and High Metal Toxicity in Soil. Science 309: 13871390
 Breitbart M et al. (2002) Genomic analysis of uncultured marine viral communities. PNAS 99: 1425014255
 Filiasi M, Livan G, Marsili M, Peressi M, Vesselli E and Zarinelli E (2014) On the concentration of large deviations for fat tailed distributions, with application to financial data. Journal of Statistical Mechanics: Theory and Experiment P09030.
 Klass OS, Biham O, Levy M, Malcai O and Solomon S (2006) The Forbes 400 and the Pareto wealth distribution. Economics Letters 90: 290295.
 Levy M and Solomon S (1996) Power laws are logarithmic Boltzmann laws. International Journal of Modern Physics C 7:595601
 Levy M and Solomon S (1997) New evidence for the powerlaw distribution of wealth. Physica A 242:9094
 Sornette D and Cont R (1997) Convergent multiplicative processes repelled from zero: power laws and truncated power laws. Journal de Physique I 7:431444
 Simon HA (1955) On a class of skew distribution functions. Biometrika 42:425440
 Bak P, Tang C, Wiesenfeld K (1987) Selforganized criticality: An explanation of the 1/f noise. Physical Review letters, 59(4): 381384.
 Newman MEJ and Sneppen K (1996) Avalanches, scaling, and coherent noise. Physical Review E 54: 62266231.
 Fisher RA (1958) The Genetical Theory of Natural Selection, 2nd rev. ed., Dover, New York.
 Orr HA (2009) Fitness and its role in evolutionary genetics. Nature Reviews Genetics 10: 531539.
 Bergstrom CT and Lanhman M (2004) Shannon Information and Biological Fitness. Proceedings of the Information Theory Workshop, IEEE 0780387201: 5054; available at doi:10.1109/ITW.2004.1405273
 Maslov S, Zhang YC (1998) Optimal investment strategy for risky assets. International Journal of Theoretical and Applied Finance 1: 377387
 Maslov S and Sneppen K (2015) Welltemperate phage: optimal bethedging against local environmental collapses. Scientific Reports 5: 10523; http://arxiv.org/abs/1308.1646
 Kussell E, Kishony R, Balaban NQ and Leibler S (2005) Bacterial persistence a model of survival in changing environments. Genetics 169: 18071814
 Balaban NQ, Merrin J, Chait R, Kowalik L and Leibler S (2004) Bacterial persistence as a phenotypic switch. Science 305: 16221625
 Maisonneuve E, CastroCamargo M and Gerdes K (2013) (p)ppGpp controls bacterial persistence by stochastic induction of toxinantitoxin activity. Cell 154:11401150.