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 reduction in size of the 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 characterized by 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 have 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.


Introduction
Mathematical description of many processes in biology and economics or finance assumes long-term exponential growth [1,2] yet no real-life environment allows growth to continue forever [3,4]. 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 often described by neutral "community drift" models [5,6] sometimes with addition of deterministic differences in efficiency of utilizing resources [7].
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 decimated by bacteriophages (see e.g. [8,9] 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 [10].

Model
Population growth P(t) of a single exponentially replicating species in a growth-limiting environment is traditionally described by Verhulst's [4] or logistic equation dP/dt = O Á P Á (1 − P/ P tot ), where the carrying capacity of the environment P tot without loss of generality can be set to 1. 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 P i . The growth stops once the global population ∑ j P j reaches the carrying capacity P tot = 1.
• The environment is periodically reseeded with new species starting from the same small population size γ ( 1 (measured in units of environment's carrying capacity).
We initially assume that growth rates and collapse probabilities of all species are equal to each other. We also disregard the neutral drift [5] 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 Figs). 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 N is kept strictly constant. In this case each extinction event is immediately followed by the introduction of a new species. We have verified that two versions of our model have very similar steady state properties. The dynamics of the fixed-N model is described by where η i (t) 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 and 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 new species with population γ ( 1 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. Fig 1A shows time-courses of populations in a system with N = 20 species and γ = 10 −4 . At certain times the carrying capacity of the environment is nearly exhausted by just one dominant species with P max ' 1 visible as dark red stripes in Fig 1A. When such dominant species goes extinct, significant fraction of resources suddenly becomes available and consequently all other populations increase by a large ratio 1/(1 − P max ). 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 Fig  1A one can clearly distinguish about 5 such waves terminated by extinctions of dominant species #5, 15, 6, 19, and 16 correspondingly. Fig 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 inverse 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). The diversity is inversely proportional to the largest population P max (t) = max i P i (t). The diversity waves (purple dashed arrows in Fig 2) are initiated when a dominating species collapses.
As a consequence, at the start of each wave the diversity abruptly increases from * 1 to a substantial fraction of the maximal possible diversity N. After this initial burst triggered by the global redistribution of biomass, the diversity exponentially declines as exp(−t/N) (the dotdashed line in Fig 2), driven by ongoing extinctions of individual populations. The typical duration, t wave 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 N Á exp(−t wave /N) * 1 or Fig 3 shows the time-aggregated distribution of population sizes for γ = 10 −9 and N = 1000 (the grey shaded area). This logarithmically-binned distribution defined by π(P) = dProb(P i (t) > P)/dlog 10 P 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 1 and N 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.
To better understand the dynamics of the system in Fig 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 P max > 1 − 1/N was eliminated, but before the available biomass was redistributed among all species. At those special moments, happening only once every t wave time steps, most population sizes are between γ and γ Á N while a small fraction reaches all the way up to * 1/N. During the rapid growth phase immediately after our snapshot was taken, all populations grow by the same factor 1/(1 − P max ) ' N 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 * N Á γ and 1 with a peak around 1/N. One can see that for γ < 0.01/N, the tail of the distribution for most abundant populations between 1/N and 1 is well fitted by a power law π(P) / 1/P τ−1 ' 1/P 0.7 (dashed line in Fig 4B) corresponding to the power law distribution of population sizes on the linear scale dProb(P i (t) > P)/dP * 1/P τ ' 1/P 1.7  demonstrate that the exponent τ for different values of γ and N is remarkably universal. Indeed, for a range of parameters where the upper and the lower peaks of π(P) are clearly separated, τ approaches a universal value τ = 1.7.

. Overall Fig 4A and 4B
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 γ ( 1/N. The passage of time t elapsed since the beginning of the current wave, leads to a decrease in the number of surviving species N surv (t) = Nexp(−t/ N), which all have the same population size P = 1/N surv (t) 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 P i > P = 1/N surv (t) is naturally given by N surv (t)/N / 1/P and thus The exponent τ = 2 predicted by this simplified model is realized in our actual model for moderately high γ * 0.1, whereas smaller values of γ give rise to a different universal exponent τ ' 1.7. 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 τ = 2 is recovered if at the start of each wave one equilibrates all species abundances to 1/N. The first section of the S1 Text in supplementary materials provides additional details on how the reduced population diversity D ¼ 1= P P 2 i < N at the start of population waves affects the exponent τ.
Two panels in Fig 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 t elapsed since the start of the current wave: −log(1 − P collapsed (t)) ' exp(t/N)/N. 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 Fig 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 Fig 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.
The population distribution in both upper and lower peaks is 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 Note that our basic model, unlike its simplified counterpart, preserves memory of population sizes distribution across several subsequent diversity waves. This is manifested e.g. in similar fractal structure of jumps sizes in waves #2-6 shown in panel B). Colors of symbols (see colorbar below) represent the log10 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. separated (as e.g. for low values of γ) the slope of the logarithmic distribution in this region has almost exactly the same value τ − 1 = 0.7 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 near-complete collapses of sub-populations of competing species. This type of dynamics strongly contrasts gradual changes implied in "community drift" neutral models [5] in ecology, or incremental random walks of stock valuations of individual companies in economics. Conversely, collapse-driven dynamics assumes sudden and usually large changes 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 can 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 Fig 5B). This overall finding is supported by a growing body of literature [11][12][13][14] where multi-modal Species Abundance Distributions (SAD) were reported for plants, birds, arthropods [14], marine organisms including single cells, corals [12], nematodes, fishes, entire seafloor communities [11], and even extinct brachiopods [15]. 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 [5] thought to apply to a variety of ecosystems including microbial communities [6,7] (see S1 Fig in supplementary materials).
Needless to say, our model is not unique in generating multimodal distributions (see e.g. [13] for other examples). Conversely, some of the variants of our model have 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 [16,17] 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 [8,18], potentially leading to their complete local extinction [9]. Phage-driven collapses do not spare bacteria with large populations and may even be biased towards such as postulated in the Killthe-Winner (KtW) hypothesis [19]. The magnitude and characteristic timescale of population changes in microbial ecosystems is still being actively discussed in the literature. While Ref. [20] reports that over half of all bacterial species in marine environments varied between being abundant and rare over a three-year period, other studies [21] 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 [22]) population shifts at the level of individual bacterial strains [8,21,23] caused by phage predation [22]. Except for interchangeable gene cassettes (metagenomic islands) responsible for either phage recognition cites or alternatively resistance to phages [24], 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 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. [9,25,26]).
On vastly longer, geological timescales, the collapse-driven dynamics of our model resembles that of species extinctions and subsequent radiations in the paleontological record [27,28]. 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 [29]. Interestingly, the extinction rate of taxonomic orders appears to be independent of their size quantified by the number of genera they contains [10], which is also one of the assumptions of our basic model.
The second application of our model is to company size 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 [30]. 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 [31]. 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 yearto-year volatility of company's market share varies with its size S, ΔS/S / S −0.2 [32].
Abundance distributions in our original model and many of its variants are characterized by a power-law tail with an exponent τ close to 2. This is in approximate agreement with the empirical data on abundance distributions of bacteria in soil samples [33], marine viruses (phages) [34]. In the economics literature, the distributions of company sizes [35], as well as those of wealth of individuals [36] are known to have similar scale-free tails. Recent data for company sizes [35] and personal wealth [36] suggest 1/P 1.8 tail of the former distribution and 1/P 2.3 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) [37][38][39], by variants of rich-get-richer dynamics [40], or in terms of Self-Organized Criticality [25,41]. The emphasis of the latter type of models on large system-wide events (avalanches [25,41] or collapses [42]) and on separation of timescales is similar in spirit to 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 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: 1. "Neutral drift model" assumes changes of population sizes during time intervals between collapses as described in Ref. [5]. In this model in addition to collapses a population of size P i randomly drifts up and down DP i / AE ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi P i ð1 À P i Þ p . The resulting diversity waves and time-aggregated distributions can be found in the supplementary S1 Fig.   2. "Exponential fluctuations model" is another variant of the neutral scenario where population sizes between collapses undergo slow multiplicative adjustments ΔP i / ±O i P i restricted by the overall carrying capacity of the environment. Details and the resulting time-aggregated distribution can be found in the supplementary S2 Fig.   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 the 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 Fig.   4. "Kill-the-Winner (KtW) model" where collapse probability c systematically increases with the population size as suggested by studies of phage-bacteria ecosystems [19]. 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 obtains a flat distribution where numbers of species in each decade of population sizes are equal to each other, see the supplementary S4 Fig.   5. "Kill-the-Loser (KtL) model", where collapse probability c systematically decreases with the population size P as c(P) * P −0.2 as suggested by the empirical studies of company dynamics [32]. As seen in the supplementary S5 Fig, 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 O i during rapid re-population phase and its own collapse probability c i . The supplementary S6 Fig 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" γ i defined by the population drop following a collapse: P i ! γ i P i . As in the previous variant each of the species is also characterized by its own growth rate O i . The supplementary S7 Fig 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. S1 Text and captions to supplementary S1-S7 Figs 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 waves, and a broad time-aggregated distribution of population sizes with scale-free tail for the most abundant species.
The classic definition of the fitness of a species in terms of its long-term exponential growth rate [43] 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 a variety of ways [44] reflecting (among other things) different timescales of growth and evolution [25], and relative emphasis on population dynamics vs. risk minimization [45]. An appropriate way to quantify species' success in a steady state system like ours is in terms of their time-averaged population size hP i (t)i t .
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 O i (model 6 and 7), its propensity to large population collapses quantified by their frequency c i (model 6), and the severity of collapses quantified by surviving fraction γ i of the population (model 7). Fig 6 shows the average population size hP i (t)i t as function of O i and γ i in the model 7. As expected, 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 c i and magnitude γ i of collapses are assumed to be independent of growth rate O i 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 fast growth even at the cost of a higher risk of collapse, whereas others 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 [2,[45][46][47] to optimize their long-term growth rate. This is obtained by splitting their populations into "growth-loving" and "risk-averse" phenotypes [45,47,48]. One example of this type of bet-hedging is provided by persistor sub-populations of some bacterial species consisting of γ * 10 −4 of the overall population [49,50] increasing to as much as γ = 10 −2 in saturated conditions (S. Semsey, private communications). A bet-hedging strategy with persistor sub-population 10 −2 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. 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 Diversity Waves in Collapse-Driven Population Dynamics 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 grows, while diversity exponentially 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. We believe that this model and its variants provide the foundation for future studies of collapse-driven dynamics in ecosystems, market economies, and social structures.