Diversity waves in collapse-driven population dynamics

Diversity waves in collapse-driven population dynamics



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 strain-specific 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 time-aggregated 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 self-organized hierarchical peak structure has a long-term memory transmitted across several waves. It gives rise to a scale-free tail of the time-aggregated 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, species-specific growth and extinction rates, and bet-hedging strategies.


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 scale-free time-aggregated 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.


Mathematical description of many processes in biology and economics or finance assumes long-term exponential growth fisher1930 (); kelly () yet no real-life 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 saturated-state 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 ().


Population growth of a single exponentially replicating species in a growth-limiting 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 rate-limiting 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 freed-up 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, S1-S7 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


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.


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 time-courses of populations in a system with species and .

Figure 1: Population dynamics. The simulated model has species and . (A) Time-courses of populations of all species. The color denotes population sizes (see the color scale on the right) with the dominating species visible as red horizontal bands. Note five diversity waves ending at purple dashed lines. Transitions between these waves were triggered by extinctions of the dominating species # 5, 15, 6, 19, 16 correspondingly. (B) The time-course of the species # 6 with the logarithmic y-axis. Note the pattern of intermittent periods of exponential growth fueled by local extinctions.

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 one-by-one 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 time-course 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 re-population 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 dot-dashed 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 one-by-one. Thus it is determined by or

Figure 2: Diversity dynamics. The grey shaded area shows the the time course of the population diversity in our model with and . Purple dashed lines mark the beginnings of diversity waves when a collapse of the dominant species with leads to an abrupt increase in population diversity from to . The diversity subsequently decreases (dash-dotted line)

Figure 3 shows the time-aggregated distribution of population sizes for and (the grey shaded area). This logarithmically-binned distribution defined by were collected over 20 million individual collapses (time-steps in our model). Thus, a time-aggregated 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 time-aggregated 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.

Figure 3: Time-aggregated population size distribution. The grey shaded area shows the time-aggregated distribution of population sizes in our model with and collected over 20 millions collapses. The green and red lines show the population size distributions collected, respectively, at the very end of each wave and at the very beginning of the next wave correspondingly as described in the text. Note that they roughly correspond to two peaks of the time-aggregated distribution.

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 time-aggregated 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 4: Time-aggregated distributions for different values of and . Time-aggregated distributions of population sizes in our model with A) and (blue), (red), and (green). B) and varying ranging between (green) to (red) in ten-fold decrements. Note the emergence of a nearly universal scale-free tail of the distribution fitted with (dashed line).

Figure 4A shows time-aggregated distributions of population sizes for and different values of ranging between ad while Figure 4B shows time-aggregated 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 scale-free 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 time-aggregated 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 sub-peaks in some waves (color bands in Figure 5B). Remarkably this multi-modal 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 time-aggregated population distribution in our model compared to its history-free version is a simple manifestation of a complex interplay between ”upstairs” and ”downstairs” subpopulations transmitting memory across waves.

Figure 5: Memory of population size distribution is preserved across several diversity waves. Time course of jumps in the logarithm of surviving populations following a collapse of a substantial population in A) the simplified model in which at the start of each wave all populations are set equal to each other; B) our basic model. Both were simulated at and . Note that our basic model, unlike its simplified counterpart, preserves memory of population sizes distribution across several subsequent diversity waves. This is manifested in similar fractal structure of jumps sizes in waves #2-6 shown in panel B). Colors of symbols represent the of the number of substantial populations during the the previous wave, when a given population originated at the small size . Thus red dots mark populations originated at the very end of the previous wave, while yellow dots - those originated when there were two large populations left in the previous wave. Finally, green, blue, and purple dots refer to older populations in the previous wave.

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.


In this paper we explore the population dynamics in saturated environments driven exclusively by near-complete collapses of sub-populations 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, collapse-driven 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 time-aggregated 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 multi-modal 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 log-binned SAD. We argue that a more reliable characterization of underlying dynamical processes can be obtained from time-series data. First, all systems capable of diversity waves are described by rapid large changes in populations of individual species. Such sudden, population-scale shifts can occur e.g. due to introduced diseases or the disappearance of keystone species paine1969 (); cohn1998 () thereby changing the composition of the entire food-web. 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 (). Phage-driven collapses do not spare bacteria with large populations and may even be biased towards such as postulated in the Kill-the-Winner (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 three-year period, other studies rodriguez-brito2010 () 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 (); rodriguez-brito2010 () caused by phage predation castberg2001 (). Except for interchangeable gene cassettes (metagenomic islands) responsible for either phage recognition cites or alternatively resistance to phages rodriguez-valera2009 () 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. Real-life 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 collapse-driven 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 end-Cretaceous 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 collapse-driven 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 year-to-year 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 power-law 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 scale-free tails. Recent data for companies filiasi2014 () and personal wealth klass2004 () suggest tail of the former distribution and tails of the latter one. Traditionally, scale-free 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 rich-get-richer dynamics simon1955 (), or in terms of Self-Organized Criticality bak1987 (); bak1993 (). The emphasis of the latter type of models on large system-wide events (avalanches bak1987 (); bak1993 () or collapses newman1996 ()) and on separation of timescales is similar in spirit to the collapse-driven dynamics in our models. A potentially important socio-economic implication of our model is that during each wave contingent sub-peaks 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 collapse-driven 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:

  1. “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 time-aggregated distributions can be found in the supplementary S1 Figure.

  2. “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 time-aggregated distribution can be found in the supplementary S2 Figure.

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

  4. ”Kill-the-Winner (KtW) model” where collapse probability systematically increase with the population size as suggested by the studies of phage-bacteria ecosystems thingstad1997 (). In this particular case the diversity dynamics and the scale-free 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.

  5. ”Kill-the-Looser (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 scale-free tail of the population distribution are both remarkably robust with respect to introduction of size-dependent collapse rate.

  6. ”Fitness model” in which each of the species has its own growth rate during rapid re-population phase and its own collapse probability . The supplementary S6 Figure show that the overall shape of the time-aggregated 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.

  7. “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 time-aggregated 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 S1-S7 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 collapse-driven dynamics such as diversity wave dynamics, and a broad time-aggregated distribution of population sizes with scale-free tail for the most abundant species.

Figure 6: Average population vs species’ properties in the “Fitness model” variant #6. Time-averaged population of a species (see color scale on the right) plotted as a function of its re-population growth rate (x-axis) and population drop after collapses (y-axis). in a variant of our model with fitness differences between species. Note that the population increase with both and . Populations and fitness parameters of species were taken from 50 million snapshots of the model.

The classic definition of the fitness of a species in terms of its long-term exponential growth rate fisher1958 () is clearly inappropriate for our model. Indeed, the long-term 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 time-averaged population size .

In the last two variants of our basic model we add fitness-related parameters to each of the species: its short-term 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 short-term 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 trade-offs.

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 bet-hedging strategies kelly (); zhang (); bergstrom2004 (); ms_kelly2015 () to optimize their long-term growth rate. This is obtained by splitting their populations into “growth-loving” and “risk-averse” phenotypes bergstrom2004 (); kussell2005 (); ms_kelly2015 (). One example of this type of bet-hedging is provided by persistor sub-populations 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 bet-hedging strategy with persistor sub-population 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 trade-off.

In this study we presented a general modeling framework for systems driven by redistribution of rate-limiting 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 one-parameter model. In spite of its simplicity this model and its variants provide the foundation for future studies of collapse-driven dynamics in ecosystems, market economies, and social structures.


Supplementary Materials

Fokker-Planck 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 :


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:


The stochastic dynamics within a single wave can thus be described by the following equation:


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 Fokker-Plank-like equation:


Here is the time-dependent 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.S2-S3 allow us to derive the parameters of the Fokker-Plank 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 Fokker-Plank equation.

The stationary solution for the time-aggregate 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


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 non-zero 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 first-principles 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:

  1. “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 (MPB-32), 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 time-aggregated species abundance distributions in this model variant for three values of and compares them to our basic model.

    Figure S1: “Neutral drift model”. This variant extends our basic model with and collapse ratio by adding the neutral drift at rate taking place between subsequent collapse events in our standard model: . The lower panel shows the time-aggregated distributions in our system simulated for collapse events. The grey shaded area refers to our basic, unmodified model, i.e. to the case, while three color lines correspond to (blue), (green), and (red). The upper four panels illustrate typical time courses of the diversity in our basic model and for three values of the color-coded as in the lower panel.
  2. “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.

    Figure S2: “Exponential fluctuations model” with , , and (green), (blue), (red) system simulated for collapse events. The grey shaded area shows the time-aggregated population distribution in our basic model, corresponding to the limit.
  3. “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 rate-limiting 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 time-aggregated distribution of local population sizes in this model variant.

    Figure S3: “ Interconnected environments model”. Figure show time-aggregated species abundance distributions in the model with environments connected by diffusion of strength simulated over collapse events. The standard model with the same parameters is shown as the grey shaded area.
  4. “Kill-the-Winner (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 ”Kill-the-Winner” (KtW) hypothesis (Thingstad TF and Lignell R (1997) Aquatic Microbial Ecology 13:19-27). 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 time-aggregated 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).

    Figure S4: “Kill-the-Winner (KtW) model” in which larger populations are preferentially targeted for collapse: . Different colors correspond to time-aggregated SADs in the model with , , and (green), (blue), and (red) simulated over collapse events. The grey shaded area refers to time-aggregated population distribution in our basic, unmodified model with the same and .
  5. “Kill-the-looser (KtL) model” in which a collapse probability declines with population size in a power-law 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: 621-633.) 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 memory-free 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 time-averaged distribution of populations thereby approaches the scaling regime described by:


    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 .

    Figure S5: “Kill-the-loser (KtL) model”in which the collapse probability declines with population size. The figure shows an , system simulated for collapse events. The upper panel illustrates the recurrent diversity waves, whereas the lower panel shows time-aggregated distributions, with the grey shaded area referring to our standard model.
  6. “Fitness model” with heterogeneous, species-specific growth rates and extinction probabilities. Each species is assigned a growth rate used when it repopulates the freed-up 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 time-aggregated population distribution in this model preserves its power-law 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 .

    Figure S6: “Fitness-model” with heterogeneous, species-specific growth rates and extinction probabilities. Each species is assigned a growth rate used when it repopulates the freed-up carrying capacity of the environment. It also has its own collapse probability . Both and are logarithmically distributed in the interval between 0.1 and 1. The purple curve in the upper panel shows the time-aggregated population distribution whereas the grey shaded area is that for the standard model where species’ growth and collapse rates are all equal to each other. The lower panel shows the average collapse probability and growth rate binned by the size of the population. he average growth rate (blue) and the average collapse probability (red shaded area) of species binned by their collected at every time step. Both curves represent time-aggregated averages as individual populations change with time.
  7. “Resilience model” where heterogeneous, species-specific 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 bet-hedging. 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: .

    Figure S7: “Resilience model” with heterogeneous, species-specific growth rates and survival ratios following a collapse. Each of the species is assigned a growth rate and collapse size , both logarithmically distributed. A) The blue curve shows the time-aggregated population distribution, whereas the grey area refers to that in our standard model from the main text. B) The average growth rate (blue) and the average survival ratio multiplied by 200 to have the same range in the plot (red shaded area) as a binned by the population size collected at every time step. Both curves represent time-aggregated averages as individual populations change with time. C) The average (arithmetic) population size as a function of species’ survivor ratio . D) The average (arithmetic) population size as a function of species’ growth rate .


  1. Fisher RA (1930) The Genetical Theory of Natural selection. Clarendon, Oxford
  2. Kelly JL (1956) A new interpretation of information rate. Bell System Technical Journal 35:917-926
  3. Malthus TR (1798) An Essay on the Principle of Population. J. Johnson, in St. Paul’s Church-yard, London
  4. Verhulst PF (1838) Notice sur la loi la population poursuit dans son accroissement. Correspondance Mathematique er Physique, 10:113-121
  5. Hubbell SP (2001) The unified neutral theory of biodiversity and biogeography (MPB-32), Princeton University Press
  6. Woodcock S, Gast VC (2007) Neutral assembly of bacterial communities. FEMS Microbiol Ecol 62: 171-180
  7. Sloan WT, Lunn M, Woodcock S (2006) Quantifying the roles of immigration and chance in shaping prokaryote community structure. Environmental Microbiology, 8: 732-740
  8. 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: 395-406;
  9. 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: 2317-2326.
  10. 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: 75-83
  11. Gray JS, Bjoergesaeter A, Ugland KI (2005) The impact of rare species on natural assemblages. Journal of Animal Ecology 74: 1131-1139.
  12. Dornelas M, Connolly SR (2008) Multiple modes in a coral species abundance distribution. Ecology Letters 11: 1008-1016.
  13. Vergnon R, van Nes EH, Scheffer M (2012) Emergent neutrality leads to multimodal species abundance distributions. Nature Communications 3: 663.
  14. Matthews TJ, Borges PA, Whittaker RJ (2014) Multimodal species abundance distributions: a deconstruction approach reveals the processes behind the pattern. Oikos 123: 533-544.
  15. Olszewski TD, Erwin DH (2004) Dynamic response of Permian brachiopod communities to long-term environmental change. Nature 428: 738-741.
  16. Paine RT (1969) A note on trophic complexity and community stability. American naturalist 103:91-93
  17. Cohn JP (1998) Understanding Sea Otters. BioScience 48: 151-155
  18. Levin BR, Stewart FM, and Chao L (1977) Resource-Limited Growth, Competition, and Predation: A Model and Experimental Studies with Bacteria and Bacteriophage., The American Naturalist 111: 3-24
  19. Thingstad TF and Lignell R (1997) Theoretical models for the control of bacterial growth rate, abundance, diversity and carbon demand. Aquatic Microbial Ecology 13:19-27
  20. 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: 12776-12781.
  21. Rodriguez-Brito 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: 739-751.
  22. 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/
  23. 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: 1971-1982.
  24. Rodriguez-Valera F, Martin-Cuadrado A-B, Rodriguez-Brito B, Pai L, Thingstad TF, et al. (2009) Explaining microbial population genomics through phage predation. Nature Reviews Microbiology 7: 828-836.
  25. Bak P and Sneppen K (1993) Punctuated Equilibrium and Criticality in a simple model of Evolution. Phys. Rev. Lett. 71: 4083-4086.
  26. Sole RV and Manrubia SC (1996) Extinction and self-organized criticality in a model of large-scale evolution. Phys. Rev. E 54, R42-45
  27. Raup DM. (1992) Extinction: bad genes or bad luck? WW Norton and Company.
  28. L. Van Valen. A new evolutionary law. Evolutionary Theory, 1:1, (1973).
  29. Sloan RE, Rigby JK, Van Valen LM, Gabriel D (1986) Gradual dinosaur extinction and simultaneous ungulate radiation in the Hell Creek Formation. Science 232: 629-633.
  30. Bachelier L. (1900) Theorie de la speculation. Gauthier-Villars, 70 pp.
  31. Fisher I (1933) The debt-deflation theory of great depressions. Econometrica 1:337-357.
  32. 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: 621-633.
  33. Gans J, Wolinsky M and Dunbar J (2005) Computational Improvements Reveal Great Bacterial Diversity and High Metal Toxicity in Soil. Science 309: 1387-1390
  34. Breitbart M et al. (2002) Genomic analysis of uncultured marine viral communities. PNAS 99: 14250-14255
  35. 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.
  36. Klass OS, Biham O, Levy M, Malcai O and Solomon S (2006) The Forbes 400 and the Pareto wealth distribution. Economics Letters 90: 290-295.
  37. Levy M and Solomon S (1996) Power laws are logarithmic Boltzmann laws. International Journal of Modern Physics C 7:595-601
  38. Levy M and Solomon S (1997) New evidence for the power-law distribution of wealth. Physica A 242:90-94
  39. Sornette D and Cont R (1997) Convergent multiplicative processes repelled from zero: power laws and truncated power laws. Journal de Physique I 7:431-444
  40. Simon HA (1955) On a class of skew distribution functions. Biometrika 42:425-440
  41. Bak P, Tang C, Wiesenfeld K (1987) Self-organized criticality: An explanation of the 1/f noise. Physical Review letters, 59(4): 381-384.
  42. Newman MEJ and Sneppen K (1996) Avalanches, scaling, and coherent noise. Physical Review E 54: 6226-6231.
  43. Fisher RA (1958) The Genetical Theory of Natural Selection, 2nd rev. ed., Dover, New York.
  44. Orr HA (2009) Fitness and its role in evolutionary genetics. Nature Reviews Genetics 10: 531-539.
  45. Bergstrom CT and Lanhman M (2004) Shannon Information and Biological Fitness. Proceedings of the Information Theory Workshop, IEEE 0-7803-8720-1: 50-54; available at doi:10.1109/ITW.2004.1405273
  46. Maslov S, Zhang Y-C (1998) Optimal investment strategy for risky assets. International Journal of Theoretical and Applied Finance 1: 377-387
  47. Maslov S and Sneppen K (2015) Well-temperate phage: optimal bet-hedging against local environmental collapses. Scientific Reports 5: 10523; http://arxiv.org/abs/1308.1646
  48. Kussell E, Kishony R, Balaban N-Q and Leibler S (2005) Bacterial persistence a model of survival in changing environments. Genetics 169: 1807-1814
  49. Balaban NQ, Merrin J, Chait R, Kowalik L and Leibler S (2004) Bacterial persistence as a phenotypic switch. Science 305: 1622-1625
  50. Maisonneuve E, Castro-Camargo M and Gerdes K (2013) (p)ppGpp controls bacterial persistence by stochastic induction of toxin-antitoxin activity. Cell 154:1140-1150.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description