Boom-bust population dynamics can increase diversity in evolving competitive
communities

Michael Doebeli∗, Eduardo Cancino Jaque∗∗, Iaroslav Ispolatov∗∗

∗ Department of Zoology and Department of Mathematics, University of British
Columbia, 6270 University Boulevard, Vancouver B.C. Canada, V6T 1Z4;
<EMAIL_ADDRESS>
∗∗ Departamento de Física, Universidad de Santiago de Chile, Santiago, Chile,
Casilla 302, Correo 2, Santiago, Chile

November 29, 2020

Running title: Boom-bust dynamics increases ecosystem diversity

Keywords: Competitive exclusion, limiting similarity, adaptive
diversification, population fluctuations, ecosystem diversity, niche packing,
temporal niche segregation, non-equilibrium population dynamics,
overcompensation

The processes and mechanisms underlying the origin and maintenance of
biological diversity have long been of central importance in ecology and
evolution. The competitive exclusion principle states that the number of
coexisting species is limited by the number of resources, or by the species’
similarity in resource use. Natural systems such as the extreme diversity of
unicellular life in the oceans provide counter examples. It is known that
mathematical models incorporating population fluctuations can lead to
violations of the exclusion principle. Here we use simple eco-evolutionary
models to show that a certain type of population dynamics, boom-bust dynamics,
can allow for the evolution of much larger amounts of diversity than would be
expected with stable equilibrium dynamics. Boom-bust dynamics are
characterized by long periods of almost exponential growth (boom) and a
subsequent population crash due to competition (bust). When such ecological
dynamics are incorporated into an evolutionary model that allows for adaptive
diversification in continuous phenotype spaces, desynchronization of the boom-
bust cycles of coexisting species can lead to the maintenance of high levels
of diversity.

The amazing diversity of life has sustained the debate about the origins and
limits of biodiversity. While random, selectively neutral processes are
thought by some to play an important role, e.g. in ecosystem dynamics 2 and in
molecular evolution, it seems that a majority of researchers would agree that
non-neutral ecological interactions – competition, predation, mutualism – are
central to understanding diversity, with competition having received the most
attention. Coexistence between competing species requires that intraspecific
competition is strong enough relative to interspecific competition. This is
captured by the concept of limiting similarity: to coexist, populations must
be sufficiently different in their resource use. If populations use the same
resource in the same way, they cannot coexist, a phenomenon known as the
competitive exclusion principle 3.

The exclusion principle has faced challenges from many empirical counter
examples, in which the number of coexisting and ecologically interacting
species was significantly higher than the number of limiting resources. The
best known such example is the Paradox of the Plankton 4, which is based on a
comparison between the relatively small number of biochemical resources
essential for plankton growth, and the number of known coexisting plankton
species, which is orders of magnitude larger. Different theoretical
explanations for conditions that circumvent the exclusion principle have been
proposed, and it is known that fluctuating population sizes can lead to
violations 3, 5, 6, 7. With fluctuating population sizes, the storage effect
8, 9, as well as relative non-linearities 3 can lead to coexistence of more
competitor species than essential resource species, e.g. because of cyclical
dominance between competitors 10. Most of these examples involve models with a
finite number of distinct resources, but coexistence due to fluctuating
population dynamics has also been shown in models with continuous niches. With
continuous, externally imposed (seasonal) periodic cycles in population sizes,
time essentially becomes an additional niche dimension along which populations
can segregate and coexist 11. In an evolutionary context, it has been shown
that limiting similarity in a continuous niche space used by an evolving
community whose member species are undergoing externally forced population
fluctuations, larger amplitude fluctuations lead to more diversity, and hence
effectively to smaller limiting similarity 12.

Most previous models used in his context have assumed that the population
fluctuations are externally imposed, and that there is a finite number of
distinct resources. Here we investigate the questions of evolving diversity
and limiting similarity in a setting where population fluctuations are not
externally imposed, but are instead due to competitive interactions within and
between the evolving species. In addition, we address the question of
diversity in continuous phenotypes spaces, corresponding to continuously
varying resource use. Rather than the 1-dimensional phenotype spaces that are
usually assumed with continuous resource distributions 12, 13, 14, our
phenotype spaces are potentially high-dimensional.

We use the models of 15, 16, 17, 18, 19, which are extensions of classical
competition models to high-dimensional continuous phenotype spaces. In
previous work, we assumed that the underlying ecological dynamics have a
stable equilibrium (the carrying capacity). However, by using difference
equations rather than differential equations to describe ecological dynamics,
it is straightforward to extend these models to allow for complicated
ecological dynamics, which are by definition endogenously generated (i.e., the
population fluctuations are a result of the competitive interactions). We show
that for certain types of endogenously generated fluctuations, which we term
“boom-bust” dynamics, the amount of diversity that evolves and is maintained
at evolutionary steady state can be much larger than the diversity maintained
without ecological fluctuations.

Boom-bust dynamics consist of long periods of (near-) exponential growth
followed by a deep crash, in such a way that the boom-bust cycles of different
species become spontaneously desynchronized. The amount of excess diversity
enabled by boom-bust dynamics increases with the dimension of phenotype space,
so that species can be much more tightly packed in high-dimensional spaces,
corresponding to a much smaller limiting similarity necessary for coexistence
in high-dimensional niche spaces.

Apart from asynchronous boom-bust dynamics, essentially all other types of
complex fluctuations, including asynchronous chaotic dynamics not exhibiting
the boom-bust fluctuations, do not increase the diversity at evolutionary
steady. Our models thus provide a specific and robust mechanism for the
evolutionary origin and maintenance of highly diverse competitive communities.

Model and simulation methods

To accommodate various types of ecological dynamics, we consider ecological
models given by difference equations, and hence set in discrete time. The
basic ecological model we use is a difference equation 20, 21, 22 that links
population densities $N$ of two consecutive generations $t$ and $t+1$,

$\displaystyle N(t+1)=F(N(t))=N(t)\frac{\lambda}{1+aN(t)^{\beta}}.$ (1)

where $\lambda>0$ is the per capita number of offspring, and $a>0$ and
$\beta>0$ are parameters describing the effect of competition. For
$\lambda<1$, $N(t)$ converges to 0 for any initial condition $N(0)>0$, and
hence extinction is the only possible outcome. We therefore assume $\lambda>1$
in what follows. In that case, model (1) has a non-zero equilibrium at
$K=((\lambda-1)/a)^{1/\beta}$, which is the carrying capacity of the
population. It is then convenient to write (1) as:

$\displaystyle N(t+1)=N(t)\frac{\lambda}{1+(\lambda-1)(N(t)/K)^{\beta}},$ (2)

as this makes it easy to formulate the model in terms of continuous phenotypes
(see below). Model (2) was shown to fit well a wide range of data 21, and for
$\beta=1$ can be derived from the logistic differential equation by
integration over a finite time interval 23.

It is well known that that the basic quantity underlying the dynamic behaviour
of this model is the derivative $dF/dN$ evaluated at the equilibrium $K$:

$\displaystyle\left.\frac{dF}{dN}\right|_{N=K}=1-\frac{\lambda-1}{\lambda}\beta.$
(3)

For $\lambda>1$, the population dynamics converges to the steady state $K$ if
and only if $|1-\frac{\lambda-1}{\lambda}\beta|<1$. Thus, for a given
$\lambda$, for small values of the exponent $\beta$ the system exhibits stable
equilibrium dynamics, and increasing $\beta$ gives way to a period-doubling
route to chaos. Biologically, increasing $\beta$ can be viewed as reflecting a
gradual change from contest to scramble competition 24. This is reflected by
the shape of the per capita number of offspring as a function of population
size, given by the right hand side of (2) divided by $N(t)$, and viewed as a
function of $N(t)$: for any $\lambda$, and for high $\beta$, the per capita
number of offspring is approximately constant until the population size
reaches the vicinity of $K$, but as $N$ increases above $K$ the number of
offspring falls rapidly to very low values, essentially generating a
population crash as soon as the population size is above $K$.

Because the difference equation (2) has three rather than two parameters, it
can exhibit certain dynamical properties that better known difference
equations, such as the Ricker equation 21, 25, do not have. For example, for
small values of $\lambda$, model (2) can exhibit highly chaotic dynamics (as
e.g. measured by the Lyapunov exponent) despite the fluctuations in population
size being (arbitrarily) small (Figure 1). Importantly, for small values of
$\lambda$, and for large enough $\beta$, model (2) exhibits ”boom-bust”
dynamics (Figure 1), in which long periods of near-exponential growth (due to
small $\lambda$) are followed by deep crashes (due to high $\beta$) once the
population size is above $K$ for the first time after the exponential phase.
This cycle repeats itself qualitatively, but the dynamics is in fact chaotic
and exhibits sensitive dependence on initial conditions, because the
population size after the crash, and hence the length of the subsequent
exponential phase, is different in each cycle. We note that such boom-bust
dynamics cannot be observed in the Ricker model. Some of the possible dynamic
regimes of model (2) and transitions between them are shown in Figure 1.

We note that there are in principle many different models that can exhibit
boom-bust dynamics (including models set in continuous time, see Discussion
section). We chose model (2) as a generic model with boom-bust dynamics for
certain parameter regions, viz. for $\lambda$-values close to 1 and large
enough $\beta$. Rather than being interested in the likelihood of a particular
model exhibiting boom-bust dynamics, we are interested in the consequences of
such dynamics for the evolution of diversity. Therefore, while pointing out
the contrast to the consequences of other types of ecological dynamics, such
as cyclic or ”regular” chaotic dynamics, delineating the different regions in
parameter space generating the different types of dynamics is not relevant for
our purposes.

Figure 1: A): Examples of population dynamics; convergence to steady state for
$\beta=2,\lambda=1.2$ (red line), stable periodic oscillations for
$\beta=16,\lambda=1.2$ (blue line), chaotic dynamics when
$\beta=20,\lambda=1.2$ (magenta line) and $\beta=5,\lambda=5$ (green line),
and boom-bust chaotic dynamics when $\beta=45,\lambda=1.2$ (black line with
individual generations shown by circles). Note the almost exponential multi-
generation growth phases in the boom-bust regime, followed by a single-
generation bust. B): bifurcation diagram for Eq. (2) for $\lambda=1.2$. Both
panels computed with $K=1$.

We now consider a generalization of Eq. (2) that includes competition between
$S$ phenotypically monomorphic populations. Each population is characterized
by its phenotype $x_{s}=(x_{s}^{1},...,x_{s}^{d})\in\mathbb{R}^{d}$,
$s=1,...,S$, where $d$ is the dimension of phenotype space (which is assumed
to be Euclidean $d$-space). The population size of phenotype $x_{s}$ is
denoted by $N_{s}$. The ecological dynamics of all $S$ clusters are determined
by the competition kernel $\alpha(x_{s},x_{r})$, which measures the
competitive impact of phenotype $x_{r}$ on phenotype $x_{s}$, and the carrying
capacity $K(x_{s})$, which is the equilibrium population size of phenotype
$x_{s}$ in the absence of any other phenotypes. (The competition kernel and
the carrying capacity are functions
$\alpha:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ and
$K:\mathbb{R}^{d}\rightarrow\mathbb{R}$, respectively.)

The discrete time dynamics of each phenotype in the competitive community is
then given by

$\displaystyle
N_{s}(t+1)=N_{s}(t)\frac{\lambda}{1+(\lambda-1)\left[\sum_{p=1}^{S}N_{p}(t)\alpha(x_{s},x_{p})/K(x_{s})\right]^{\beta}},$
(4)

$s=1,...,S$. The sum in the denominator on the right hand side of (4) is the
effective population size experienced by phenotype $x_{s}$. Eq. (4) is a
discrete time analog of the continuous-time many-species logistic competition
model in multidimensional phenotype space that was used in several previous
articles 17, 18, 19. In contrast to the continuous time models used
previously, in the discrete time model populations can undergo ecological
fluctuations and sudden collapses not only after a population itself exceeds
the carrying capacity, but also when the cumulative competition from other
phenotypes is strong enough (i.e, when the effective population size is above
$K$).

For simplicity, and following 16, 17, 18, we used the following functions for
the competition kernel and the carrying capacity,

$\displaystyle\alpha(x,y)=\exp\left[-\sum_{i=1}^{d}\frac{(x^{i}-y^{i})^{2}}{2\sigma_{\alpha}^{2}}\right],\>K(x)=K_{0}\exp\left[-\frac{\sum_{i=1}^{d}(x^{i})^{4}}{4\sigma_{K}^{4}}\right].$
(5)

Thus, competition is symmetric and strongest between phenotypes that are
similar, and the carrying capacity has a unique maximum $K_{0}$ at 0. We set
the scaling parameters $K_{0}$ and $\sigma_{K}$ to $K_{0}=1$ and
$\sigma_{K}=1$. (We note that in general, using Gaussian forms for both the
competition kernel and the carrying capacity can result in structural
instabilities 26).

To simulate the evolutionary process, we start with a number $S$ of phenotypes
(typically, $S=1$, and the phenotype is randomly chosen in the vicinity of the
maximum of the carrying capacity). We then simulate the ecological dynamics in
discrete time, using (4) for each of the phenotypes. In each generation, a new
phenotype is generated with a probability $\mu$ (typically $\mu=0.1$). The new
phenotype is a mutant of one of the existing phenotypes. Of those, a parental
phenotype is chosen with a probability proportional to its population size,
and the offspring phenotype is chosen randomly from a Gaussian distribution
with the average centered at the parental phenotype and a small standard
deviation $\Delta x$ (typically $\Delta x=10^{-2}$).

After addition of the new phenotype, the community now comprises $S+1$
phenotypes, and the process is repeated for many generations. What one wants
to know from this process is how the distribution of phenotypes changes over
time. To keep the number of phenotypes from increasing to very large numbers
that would render the simulation computationally impossible, we periodically
merge phenotypes that are very close together. Specifically, once every
$t_{merge}$ generations (typically $t_{merge}=1000$), phenotypes that are
within a distance $\Delta x$ of each other are merged (preserving their
phenotypic center of mass) and their population sizes added. In addition,
every generation all clusters with populations densities below a threshold
(typically $=10^{-12}$) are declared extinct and removed from the system.
Together, these procedures preserve the phenotypic variance necessary for
evolution, but prevent undesirable computational complexity.

To define and count the number of phenotypically distinct species in the
community at any given point in time, with each species possibly consisting of
a number of similar phenotypes, the phenotypes in the community are clustered
with a larger distance $\Delta x_{species}$ (typically $\Delta
x_{species}=10^{-1}$). This phenotypic distance is still significantly smaller
than the typical scales of ecological interactions as long as
$\sigma_{\alpha}$ and $\sigma_{K}$ in (5) are of order 1. This implies that
the phenotypes within a designated species experience very similar competitive
interactions and generally follow the same population dynamics. Note that
species designation is only used to gather statistical data from simulations,
but not in the actual computational steps of the simulations.

Results

For $\sigma_{\alpha}<\sigma_{K}$ in (5), the continuous time analog of the
model presented above undergoes adaptive diversification and radiation into a
steady state species distribution 14, 18, 27, see also 28, 29. For example, if
we set $\sigma_{\alpha}=0.5$ and $\sigma_{K}=1$ in (5), then for $\beta=1$
system (4) is equivalent to the corresponding continuous time system 23, and
in a 2-dimensional phenotype space undergoes diversification into a stable
community of 16 coexisting phenotypic clusters (species) with approximately
constant population sizes. This is illustrated in Figure 2a and the
corresponding video. As long as $\beta=1$, the observed diversification is
independent of $\lambda$ and $K_{0}$. The main purpose of this paper is to
explore the effect of increasing $\beta$ to values $>1$, which eventually
makes the local dynamics (2) unstable. As the exponent $\beta$ is increased,
stationary populations lose stability and, similarly to the single-species
model shown in Fig. 1, the ecological dynamics of populations in an evolving
community become first periodic and then chaotic (see below for specific
examples of non-equilibrium dynamics).

Figure 2:  Snapshots of cluster distributions for low (left column) and high
(right column) $\beta$-values after $10^{7}$ generations, for 2-dimensional
(top row) and 3-dimensional (bottom row) phenotype spaces. A): 16 species for
$\beta=1,\lambda=1.2$, $\sigma_{\alpha}=0.5$; B): 47 species for
$\beta=45,\lambda=1.2$, $\sigma_{\alpha}=0.5$; C): 8 species for
$\beta=1,\lambda=1.2$, $\sigma_{\alpha}=0.75$; D) 20 species for
$\beta=55,\lambda=1.2$, $\sigma_{\alpha}=0.75$. The corresponding videos of
the evolutionary process that led to these configurations from single
ancestors can be found at figshare.com/s/f2d8ecf480fa372319e1. In the top row
panels and corresponding videos (2-dimensional case), all phenotypes present
in the evolving community are shown (represented as dots). In the bottom row
panels and corresponding videos (3-dimensional case), species resulting from
clustering of populations of similar phenotypes using a merging distance of
$\Delta x_{species}=10^{-1}$ (see main text) are represented by circles whose
size is proportional to a species’ population size.

For low intrinsic growth rates $\lambda$ this has profound effects on the
amount of diversity in the system, as illustrated in Figure 2 (and
accompanying videos). For such $\lambda$-values, increasing $\beta$ in the
local dynamics (2) has the effect of eventually inducing pronounced boom-bust
population dynamics (cf. Figure 1). In an evolving community, increasing
$\beta$ induces boom-bust dynamics in each of the phenotypes present in the
community, with each phenotypic cluster (species) undergoing multiple
generations of (near-)exponential population growth punctuated by deep
crashes. These ecological dynamics unfold in such a way that the dynamics of
neighbouring clusters of phenotypes are desynchronized, i.e., such that
crashes and subsequent exponential growth phases occur at different time
points. With such desynchronized boom-bust ecological dynamics, evolution can
generate a drastic increase in diversity compared to that evolving in
ecologically stable communities (Figure 2). Increased diversity due to boom-
bust ecological dynamics typically occurs as long as the ecological conditions
for adaptive diversification due to frequency-dependent competition are met
(i.e., as long as $\sigma_{\alpha}<\sigma_{K}$ in (5)). While the amount of
diversity depends on the exact values of $\sigma_{\alpha}<\sigma_{K}$,
significantly more diverse communities tend to evolve with boom-bust dynamics
than with stable equilibrium ecological dynamics.

Figure 3: The number of coexisting species as a function of the exponent
$\beta$ for $\lambda=1.2$ (circles), $\lambda=1.6$ (squares), $\lambda=2$
(diamonds), $\lambda=5$ (triangles). Dimension of phenotype space is $d=2$,
and $\sigma_{\alpha}=0.5$. Species are counted after $10^{7}$ generation and
after clustering of populations of similar phenotypes using a merging distance
of $\Delta x_{species}=10^{-1}$ (see main text).

Figure 3 shows the number of species coexisting at the evolutionary saturation
state as a function of the parameter $\beta$ for different values of
$\lambda$. The figure illustrates that $\lambda$ has to be small enough for a
substantial increase in diversity to be observed for high $\beta$. Indeed, in
the local model (2) boom-bust dynamics can only be observed for small
$\lambda$, and it is exactly these kinds of population dynamics that allow for
increased diversity. For larger $\lambda$, increasing $\beta$ results in more
“traditional” forms of chaotic dynamics with irregular, high frequency
oscillations of increasing amplitude. Such local dynamics also lead to chaotic
ecological dynamics in populations comprising an evolving community, but they
do not generate an increase in the diversity that can evolve and be
maintained. In general, as $\beta$ is increased to very high values, the model
becomes less relevant biologically: the population crashes become very severe,
which results in extinctions that are frequent enough for diversity not to be
able to evolve anymore (see below).

Figure 4: Examples of population dynamics of 5 neighbouring phenotypes in a
community consisting of the 25 phenotypes indicated by color in the inset.
Parameter values are $\beta=30,\lambda=1.6$, and $\sigma_{\alpha}=0.5$, which
corresponds to the right-most square with more than 1 species in Figure 3.
Dimension of phenotype space is $d=2$. With equilibrium population dynamics
($\beta=1$) the 25 phenotypes cannot coexist. Note that the population
dynamics of immediate neighbours are desynchronized.

Figure 4 illustrates the desynchronized boom-bust dynamics in an artificial
community of 25 species, with each species represented by a single phenotype,
and such that the phenotypes are arranged on a regular grid in phenotype space
(see inset in Figure 4). This community corresponds to the case indicated by
the right-most square with more than 1 species in Figure 3, in which $\beta$
is large enough for the diversity to increase to 25 coexisting species, rather
than the 16 that would evolve for lower $\beta$. Figure 4 shows the population
dynamics of a subset of 5 phenotypes arranged on a line in the grid, resulting
from simulating the ecological dynamics (4) of the whole community of 25
phenotypes. For each phenotype, the dynamics exhibits boom-bust cycles, and
neighbouring cycles are desynchronized.

It is worth noting that the increased diversity seen for higher $\beta$-values
occupies approximately the same phenotypic range as the lower diversity for
lower $\beta$-values (Figure 2). This implies that with higher diversity, the
different species are more closely packed in phenotype space, and hence that,
generally speaking, conditions of limiting similarity are relaxed in the boom-
bust dynamic regime. Because of lower thresholds for limiting similarity,
i.e., denser packing, the increase in diversity at evolutionary stationary
state due to boom-bust ecological dynamics becomes more pronounced with higher
dimensions of phenotype space, as illustrated in Figure 5.

Figure 5: The number of species vs. $\beta$ for $\lambda=1.2$,
$\sigma_{\alpha}=0.5$ and three different dimensions of phenotype space, $d=1$
(squares), $d=2$ (circles), and $d=3$ (diamonds). Species are counted after
$10^{7}$ generation and after merging the populations of similar phenotypes
using a merging distance of $\Delta x_{species}=10^{-1}$ (see description in
main text). Note that the level of diversity shown for small $\beta$-values
corresponds to the diversity evolving with stable equilibrium ecological
dynamics.

Relaxed limiting similarity conditions require a decrease in competitive
pressure that species in neighbouring regions of phenotype space exert on each
other. Such a decrease can be achieved if neighbouring populations are
fluctuating in opposite phases, as shown for an artificial example in Figure
4. Figure 6 illustrates that in full simulations of the evolutionary process,
neighbouring species indeed generally exhibit such an anti-correlation for
high $\beta$-values. Essentially, the anti-correlation between populations of
neighbouring species stems from the asynchrony of their boom-bust cycles. In
the Supplementary Material, we show that such desynchronization is expected to
emerge spontaneously from an arbitrary small initial difference between
populations: in a simple idealized configuration of two competing species with
boom-bust dynamics, their population sizes converge to a state of complete
anti-synchronization.

Figure 6: The correlation between population sizes of phenotypes $x$ and $y$,
$C_{xy}\equiv\langle(N_{x}(t)-\langle N_{x}\rangle)(N_{y}(t)-\langle
N_{y}\rangle)\rangle$ as a function of phenotypic distance $|x-y|$ for
$\beta=50$, $\lambda=1.2$ (black line), $\beta=20$, $\lambda=1.2$ (red line)
and $\beta=5,\lambda=5$ (blue line). Dimension of phenotype space is $d=2$,
and $\sigma_{\alpha}=0.5$. The correlation was calculated by taking into
account all possible pairs of phenotypes over $5\times 10^{6}$ generations
after the steady state level of diversity was reached. Anti-synchronization is
only seen for boom-bust dynamics (black) allowing for increased diversity (cf.
Figure 3), but not for higher frequency chaotic dynamics with small (red) or
large (blue) amplitudes.

For smaller values of $\beta$ or larger values of $\lambda$, for which
populations do not undergo boom-bust dynamics this anti-correlation is not
seen (Figure 6).

The explanation for higher diversity based on anti-correlated boom-bust cycles
of phenotypically close species suggests that to make this mechanism work,
these cycles should be of sufficient length. This means that the population
crashes should be sufficiently severe (large $\beta$), and the intrinsic
growth rate $\lambda$ should be sufficiently small. Essentially, the
exponential phase should be long enough for robust desynchronization. This
effect cannot be achieved with high intrinsic growth parameters $\lambda$
(Figure 3): the increase in diversity is noticeably diminished for
$\lambda=1.6$, and is absent for larger $\lambda$. In particular, the type of
chaotic population fluctuations induced by high $\beta$ for $\lambda$-values
that are significantly larger than $1$ do not lead to increased diversity,
because for larger $\lambda$, complex dynamics are not of the boom-bust type.

To confirm that the boom-bust cycles, rather than chaoticity or other features
of the population dynamics defined by Eq. (4), are the essential mechanism for
the observed increase in diversity, we stripped the model (4) from all other
features except its ability to run boom-bust cycles: We assumed that each
phenotypic population grows exponentially with an exponent $\lambda$ until the
effective density experienced by a given phenotype, i.e., the cumulative
competitive effect of all phenotypes, given by the term
$\sum_{p=1}^{S}N_{p}(t)\alpha(x_{s},x_{p})$ in denominator of (4), becomes
greater than the carrying capacity of that phenotype. When that happened, the
population of that phenotype was reduced to a small fraction of its population
size, simulating a severe crash. The mutation and merging procedures were
implemented as in the original model. This modified model shows qualitatively
very similar results (not shown) and exhibits significant increases in
diversity at a level very similar to the original model, as long as
$\lambda$-values are close to 1, so that the exponential phase starting from
low densities is long and slow, and as long as the population crashes are
severe enough. This confirms that the key for the evolution of higher
diversity is the existence of pronounced boom-bust dynamics for all phenotypic
populations.

An interesting question concerns the effect of the frequency and size of
mutations on diversity. These were assumed to be $\mu=0.1$ and $\Delta
x=10^{-2}$ for the results presented so far (note that it is really only the
product of these two parameters that matters). A reduction in $\mu$ and/or
$\Delta x$ slows down evolution in general and diversification in particular.
This effect is illustrated in Figure S.3A, where we show the number of species
vs. time for 4 different mutation frequencies. Smaller mutation rates result
in longer times to reach the equilibrium level of diversity. There is,
however, another, less direct effect of mutation rate on diversity. For any
non-zero extinction threshold and even moderate $\beta$, there is a small but
finite probability that all phenotypes of a well-developed cluster, and hence
the corresponding species, go extinct during a particularly severe bust. The
extinct cluster can eventually get replaced by newly arising mutants, but the
time it takes mutations to undergo a sufficient number of phenotypic steps to
reach the vacated spot in phenotype space depends on $\mu$ and $\Delta x$. For
any given mutation rate and size, these processes may equilibrate at different
levels of diversity. In particular, lower extinction thresholds (making
extinction less likely) lead to higher levels of diversity at saturation. This
is illustrated in Figure S.3B.

In general, diversity decreases drastically for very high $\beta$-values and
eventually the system is reduced to just a single phenotypic cluster. This
occurs because with large $\beta$-values, the crashes due to the effective
density experienced being higher than the carrying capacity become
progressively more severe and can bring all phenotypes comprising one species
below the minimum population threshold, thus rendering the species extinct.
Even though diversification is still favoured by selection, the rate of
species extinction for high $\beta$-values is too high for diversity to
evolve. The very dynamic regime of this “competition” between extinction and
diversification is illustrated in Fig. S.2.

Discussion

We propose a possible explanation for the emergence and persistence of large
amounts of diversity based on competition models for evolving communities with
fluctuating population dynamics. When these fluctuations are in the boom-bust
regime, in which long periods of exponential growth are followed by deep
population crashes, diversity in continuous phenotype spaces evolves well
beyond what is expected based on limiting similarity with stationary
ecological dynamics. The key mechanism that results in higher diversity is the
spontaneous de-synchronization of boom-bust cycles between phenotypically
similar species, which essentially reduces interspecific competitive impacts
and allows for much denser packing of species in niche space.

Population fluctuations have long been considered as a potential mechanism
leading to violations of the competitive exclusion principle. For the most
part, past studies have either assumed a fixed set of resources 3, 10, 30, or
they have assumed externally imposed fluctuations 11, 12. In such models, the
mechanism of ecological fluctuations causing an increase in diversity can be
viewed as a form of the temporal storage effect 8, 9, which intuitively
corresponds to temporal segregation in niche space 11. In fact, there have
also been models showing that population fluctuations can decrease diversity
in an evolutionary context 31, but these models appear to allow for jack-of-
all-trades mutations on a finite set of resources, which can increase rather
than decrease the amount of interspecific competition in the system.

Our models extend previous models for the emergence and maintenance of
diversity under stable equilibrium ecological dynamics 14, 18, 28, 29, 27.
They differ from earlier models such as 28, 29 in key aspects: they consider
evolution in high-dimensional phenotypes that characterize continuously
variable and multivariate niche use, and persistent ecological fluctuations
are intrinsically generated by overcompensating competition. Desynchronized
boom-bust cycles provide a robust mechanism for a substantial increase in the
diversity that can evolve and be maintained in such models, an effect that
increases with increasing dimension of phenotype space. We note that this
latter result is not obvious, as with higher phenotypic dimensions the number
of phenotypically similar species (nearest neighbours in phenotype space)
increases linearly with the dimension, which may be expected to make de-
synchronization of neighbouring boom-bust cycles more difficult due to denser
phenotype packing.

This mechanism of “diversification in time” is similar to those previously
reported 11, 30: time acts as additional niche space, and separation along
this niche space can alleviate interspecific competition. In the language of
32, boom-bust desynchronization effectively increases the ”environmental
dimension”, which is a determinant of the amount of diversity that can be
sustained. Again, this is akin to the temporal storage effect 8, 9, although
the latter is mostly invoked for externally generated population fluctuations.
The longer the boom-bust cycles, the more temporal separation between similar
species is possible. If the population crashes in the boom-bust regime become
too severe, they produce frequent extinctions, which eventually leads to a net
negative effect on diversity.

In our models, higher diversity can only be observed in the presence of
pronounced boom-bust cycles, but not with other types of population
fluctuations, such as periodic or chaotic dynamics with high-frequency
oscillations. From a modeling perspective, it is worth noting that more
standard and more widely used discrete-time models, such as the Ricker
equation or the discrete-time logistic model, even in their chaotic regimes
cannot exhibit the type of chaotic boom-bust dynamics that model (2) exhibits
for low intrinsic growth rates $\lambda$ and large (overcompensating) $\beta$.
This reiterates old cautionary notes about the judicious use of discrete maps
for modeling ecological dynamics 22.

Discrete-time models have proved to be very useful for many different purposes
in ecology and evolution at least since Ricker’s s famous stock and
recruitment paper 25. However, we note that boom-bust dynamics can also be
generated using continuous-time models. To illustrate this, consider a
continuous-time analogue of the modified model introduced at the end of the
Results section. This continuous-time model has two phases, representing slow
exponential growth and fast exponential decline in continuous time. In the
first phase, long and slow exponential growth occurs from low densities for
each phenotypic population, while keeping track of the effective density
experienced by each phenotype, i.e., of the weighted sum over all phenotypic
population sizes, with weights given by the competition kernel (this
corresponds to the sum in the denominator of eq. (4)). Once the effective
density of a given phenotype reaches the carrying capacity of that phenotype,
there is a very fast exponential decline until the phenotype reaches a small
fraction of the population size it had before the decline, which corresponds
to a severe population crash. Simulations of this simple boom-bust model in
continuous time (results not shown) produce qualitatively identical results:
the amount of diversity that emerges and is maintained evolutionarily is much
larger than the diversity that would evolve with stable equilibrium dynamics
(as e.g. reported using a continuous-time logistic model in 18). This again
underscores the generality of the effects of boom-bust ecological dynamics on
diversity.

We speculate that the mechanisms and results reported here are not limited to
competition models, but could also be manifest in communities with other
ecological interactions, e.g. in communities with crowding effects 33, or in
communities containing both predators and prey 34. Whenever the population
dynamics exhibit patterns of rapid growth interspersed by crashes (as may e.g.
be expected in many predator-prey systems), temporal desynchronization can
occur spontaneously and thus lead to increased diversity. Such effects were
shown in 34, who reported that ”kill-the-winner” mechanism, in which predation
generates crashed in the most abundant consumer species, can generate
increased levels of diversity. These mechanism differ from the ones reported
here in that they are extrinsic to the crashing consumer species (and it is
difficult to compare those system to baseline systems with stable ecological
dynamics).

There is some empirical support for the effect of boom-bust cycles on
ecosystem diversity. For example, such patterns were observed in carefully
staged long-running experiments with several plankton species 7, and the
experimental data showed that out-phase oscillations in predator-prey cycles
of zooplankton and phytoplankton were important for the maintenance of
diversity in this system 35. Predation from pathogens have also been reported
to induce algal boom-bust cycles 36. Generally, boom-bust cycles appear to be
common in many marine ecosystems, which are known to be very diverse. For
example, it has been suggested to call echinoderms a “boom-bust” phylum 37,
and recent work shows that in polar plankton communities, which constitute an
important ecosystem in the global ocean, phytoplankton dynamics are often
categorized by “boom-bust” cycles 38.

It is interesting to put our results in the context of observations of
“neutral evolution”. For example, 39 report neutral taxonomic distributions
during early metazoan diversification into relatively empty niche space. In
our models, such expansions could be classified as the boom stage, and
according to our model assumption would then indeed occur essentially unabated
and in the absence of competitive effects. The actual selection only occurs
during the bust stage with populations of less adapted species crashing
earlier and deeper. Such an application of the model (4) would be rather
speculative however, as the unrestricted exponential growth phase would have
to last for a very long time and would in any case represent a simplified and
unrealistic assumption for such scenarios. We also note that 39 consider
neutrality based on taxonomic data, not on functional data, whereas our model
only considers functional phenotypic data. It has been noted that the
distinction between taxonomic and functional data is very important in many
microbial ecosystems 40, and in particular that functional data can be
decidedly non-random even when taxonomic distributions look random.

Overall, we think that our results provide a useful evolutionary perspective
for thinking about diversity in natural ecosystems. Boom-bust population
fluctuations are a robust, intuitively appealing and probably under-
appreciated potential cause of significantly increased diversity in evolving
ecosystems.

Acknowledgement

MD was supported by NSERC Discovery Grant 219930. ECJ was supported funded by
the National Agency for Research and Development (ANID) Scholarship 21190785.
II acknowledges support from FONDECYT project 1200708.

Supplementary Material

Spontaneous desynchronization of coupled boom-bust populations

Here we present analytical arguments to illustrate that for large $\beta$ and
populations with sufficiently distinct phenotypes, and thus only moderate
mutual competitive effect, the boom-bust cycles inevitably desynchronize.

Consider two phenotypes $x_{1}$ and $x_{2}$ with population sizes $N_{1}$ and
$N_{2}$ whose phenotypes symmetrically diverged from the maximum of the
carrying capacity, so that $x_{1}=-x_{2}\equiv x$. Thus the carrying capacity
for each species is equal, $K(x_{1})=K(x_{2})\equiv K$. The population
dynamics of the first species is defined as

$\displaystyle N_{1}(t+1)=\frac{\lambda
N_{1}(t)}{1+(\lambda-1)\left\\{[N_{1}(t)+\alpha
N_{2}(t)]/K\right\\}^{\beta}},$ (S1)

where $\alpha\equiv\alpha(x_{1},x_{2})$ is the competitive effect that the two
species have on each other (recall that $\alpha$, given by (4) in the main
text, is assumed to be symmetric). The population dynamics of the second
species is defined analogously.

For argument’s sake, we assume that initially $N_{1}$ and $N_{2}$ are small,
and the second population is slightly larger than the first one,
$N_{1}(0)\equiv N$, $N_{2}(0)=N+\Delta N$. During the boom phase, the
effective population size experienced by phenotype $x_{1}$, $N_{1}(t)+\alpha
N_{2}(t)$, is (much) less than one, so that when elevated to the large power
$\beta\gg 1$, this term becomes negligible compared to 1. Then the denominator
of (S1) is very close to one. The same is true for $N_{2}$, so that both
populations grow approximately exponentially,

$\displaystyle N_{1}(t)\approx\lambda^{t}N$ (S2) $\displaystyle
N_{2}(t)\approx\lambda^{t}(N+\Delta N).$

We assume that the initial difference between the two population sizes is
small enough for both of them to crash in the same time step. The populations
of two species immediately after the crash step $T$ are

$\displaystyle
N_{1}(T+1)=\frac{\lambda^{T+1}N}{1+(\lambda-1)\left(\frac{\lambda^{T}[N+\alpha(N+\Delta
N)]}{K}\right)^{\beta}},$ (S3) $\displaystyle
N_{2}(T+1)=\frac{\lambda^{T+1}(N+\Delta
N)}{1+(\lambda-1)\left(\frac{\lambda^{T}[N+\Delta N+\alpha
N]}{K}\right)^{\beta}}.$

We introduce the notation

$\displaystyle C(t)\equiv\frac{N_{1}(t)+\alpha N_{2}(t)}{K},$ (S4)

and ignore the summand 1 in the denominators of (S3), as at the crash step
$C(T)$ exceeds 1 (for the first time in the cycle), so that for $\beta\gg 1$,
the second term in the denominators becomes very large, $C^{\beta}(T)\gg 1$.
Thus:

$\displaystyle N_{1}(T+1)$
$\displaystyle\approx\frac{\lambda^{T+1}N}{(\lambda-1)C^{\beta}(T)},$ (S5)
$\displaystyle N_{2}(T+1)$ $\displaystyle\approx\frac{\lambda^{T+1}(N+\Delta
N)}{(\lambda-1)\left[C(T)+\frac{(1-\alpha)\lambda^{T}\Delta
N}{K}\right]^{\beta}}=$
$\displaystyle=\frac{\lambda^{T+1}N}{(\lambda-1)C^{\beta}(T)}\times\frac{1+\frac{\Delta
N}{N}}{\left[1+\frac{(1-\alpha)\lambda^{T}\Delta N}{KC(T)}\right]^{\beta}}.$

Expanding expression (S5) for $N_{2}(T+1)$ for small $\Delta N$ and observing
that the first fraction in this expression is $N_{1}(T+1)$ and
$C(T)\approx\lambda^{T}N(1+\alpha)/K$, we get

$\displaystyle N_{2}(T+1)\approx N_{1}(T+1)\left\\{1-\frac{\Delta
N}{N}\left[\beta\frac{1-\alpha}{1+\alpha}-1\right]\right\\}.$ (S6)

For $\beta\gg 1$ and $\alpha$ noticeably less than 1, the factor in square
brackets in (S6) is larger than 1, so that the difference between the two
populations $N_{1}$ and $N_{2}$ after the bust will change sign,
$N_{2}(T+1)<N_{1}(T+1)$, and will grow in absolute value. Thus, after several
boom-bust cycles, the difference $|N_{1}-N_{2}|$ becomes large enough for the
minority population to avoid crashing at the same time step as the majority
population, and hence to delay its crash by at least one step, as it is
illustrated in Fig. S.1.

Figure S.1: Illustration of the desynchronization of two populations described
by (S5) with initially very similar population sizes. Parameter values
$\beta=40$ and $\alpha=0.7$.

In fact, the lag between the two populations increases until reaching
approximately half the boom-bust period, thus rendering the two populations
maximally anti-correlated. This can be seen by assuming that in each crash the
population is reduced to the initial population size $N$. This is an
approximation that ignores the chaotic nature of the population dynamics, but
serves to demonstrate the general trend. We denote the initial lag between the
first and the second population crashes as $\Delta t$ steps. Since the first
population crashed and resumed its growth earlier, its next crash will come no
later than that of the second population. It will occur at step $T$, where $T$
is the smallest integer that satisfies

$\displaystyle N\lambda^{T}\left(1+\frac{\alpha}{\lambda^{\Delta
t}}\right)\geq K.$ (S7)

Here we have taken into account that the second population at time $T$ is
$N\lambda^{T-\Delta t}$. After the first population crashes down to $N$, the
second population will continue growing for $\Delta s$ more steps and will
crash at the smallest $\Delta s$ satisfying the condition

$\displaystyle N\lambda^{T-\Delta t+\Delta s}\left(1+\alpha\lambda^{\Delta
s}\right)=N\lambda^{T-\Delta t+\Delta
s}\left(1+\frac{\alpha}{\lambda^{T-\Delta t}}\right)\geq K.$ (S8)

If $T-\Delta t>\Delta t$ (which means $\Delta t<T/2$), that is, if the initial
lag between crashes is smaller than half the cycle length, the two conditions
(S7, S8) can be simultaneously satisfied only when $\Delta s>\Delta t$, i.e.
when the lag between crashes increases. The lag stops increasing when $\Delta
t=T/2$, i.e., when the boom-bust cycles of the two populations become
maximally desynchronized.

In a general sense, the anti-synchronization of boom-bust cycles can be
thought of as the opposite of the well-analyzed phenomenon of synchronization
of biological oscillators 41. The simplest explanation for this disparity are
the opposite convexities of growth functions used in these models.

Effect of the extinction threshold on diversity

To show that the reduction in diversity for very large $\beta>60$ is indeed
caused by prevalence of cluster extinction over diversification, we performed
simulations varying the extinction threshold. Fig. S.2 (and accompanying
videos) shows the evolving community for extinction thresholds $10^{-12}$
(used for figures in the main text), $10^{-14}$ and $10^{-10}$. Three
snapshots of the evolving system at their steady state diversity are shown,
illustrating that changing the extinction threshold shifts the balance between
extinction and diversification: lower thresholds increase the level of
diversity, and with high enough thresholds (and high enough $\beta$) diversity
is reduced to a single species.

Figure S.2: Snapshots of phenotype distributions after $10^{7}$ generations
for three values of the extinction threshold: A) extinction threshold
$10^{-12}$ used for figures in main text; B) lower extinction threshold
$10^{-14}$; C) higher extinction threshold $10^{-10}$. All three scenarios are
simulated with very strong overcompensation, $\beta=110$, and $\lambda=1.2$.
Dimension of phenotype space is $d=2$, and $\sigma_{\alpha}=0.5$. The videos
of the evolutionary diversification processes that generated these
configurations can be found at figshare.com/s/f2d8ecf480fa372319e1.

Effect of the mutation rate on diversity

Figure S.3A shows the number of species as a function of time (in generations)
for different mutation rates (other parameters are as for Figure 2B in the
main text). The lower the mutation rate, the longer it takes for diversity to
reach saturation levels. Figure S.3B illustrates that the extinction threshold
has an effect on the saturation levels of diversity: lower thresholds (fewer
extinctions) lead to elevated levels of diversity. For Figure S.3B, the
mutation rate was $\mu=0.005$, so that one mutation of typical size less than
1% of the parental phenotype occurs every 200 generations in the entire
evolving community.

Figure S.3: A: Number of species in the community as a function of time for
different mutation rates $\mu$. Species were counted using the clustering
algorithm described in the Model and simulation methods section. Other
parameter values were as for Fig. 2B in the main text. B: Number of species in
the community as a function of time for two different extinction thresholds
and for a mutation rate $\mu=0.005$. Other parameter values were as for Fig.
2B in the main text. For the two scenarios shown in B, videos of the
evolutionary process can be found at figshare.com/s/f2d8ecf480fa372319e1.

## References

  * 1
  * Rosindell _et al._ 2012 Rosindell, J., Hubbell, S.P., He, F., Harmon, L.J. & Etienne, R.S. (2012).  The case for ecological neutral theory.  _Trends in Ecology & Evolution_, 27, 203–208.
  * Armstrong & McGehee 1980 Armstrong, R.A. & McGehee, R. (1980).  Competitive exclusion.  _The American Naturalist_ 115, 151–170.
  * Hutchinson 1961 Hutchinson, G.E. (1961).  The paradox of the plankton.  _The American Naturalist_ , 95, 137–145.
  * Koch 1974a Koch, A.L. (1974a).  Coexistence resulting from an alternation of density dependent and density independent growth.  _Journal of Theoretical Biology_ , 44, 373–386.
  * Koch 1974b Koch, A.L. (1974b).  Competitive coexistence of two predators utilizing the same prey under constant environmental conditions.  _Journal of Theoretical Biology_ , 44, 387–395.
  * Benincà _et al._ 2008 Benincà, E., Huisman, J., Heerkloss, R., Jöhnk, K.D., Branco, P., Van Nes, E.H. _et al._ (2008).  Chaos in a long-term experiment with a plankton community.  _Nature_ , 451, 822–825.
  * Chesson 1994 Chesson, P. (1994).  Multispecies competition in variable environments.  _Theoretical Population Biology_ , 45, 227–276.
  * Chesson 2000 Chesson, P. (2000).  Mechanisms of maintenance of species diversity.  _Annual Review of Ecology and Systematics_ , 31, 343–366.
  * Huisman & Weissing 1999 Huisman, J. & Weissing, F.J. (1999).  Biodiversity of plankton by species oscillations and chaos.  _Nature_ , 402, 407–410.
  * Barabás _et al._ 2012 Barabás, G., Meszéna, G. & Ostling, A. (2012).  Community robustness and limiting similarity in periodic environments.  _Theoretical Ecology_ , 5, 265–282.
  * Kremer & Klausmeier 2017 Kremer, C.T. & Klausmeier, C.A. (2017).  Species packing in eco-evolutionary models of seasonally fluctuating environments.  _Ecology Letters_ , 20, 1158–1168.
  * Roughgarden 1979 Roughgarden, J. (1979).  _Theory of Population Genetics and Evolutionary Ecology: An Introduction_.  MacMillan, New York.
  * Doebeli 2011 Doebeli, M. (2011).  _Adaptive diversification (MPB-48)_.  vol. 48.  Princeton University Press.
  * Doebeli & Ispolatov 2010 Doebeli, M. & Ispolatov, I. (2010).  Complexity and diversity.  _Science_ , 328, 494–497.
  * Doebeli & Ispolatov 2014 Doebeli, M. & Ispolatov, I. (2014).  Chaos and unpredictability in evolution.  _Evolution_ , 68, 1365–1373.
  * Ispolatov _et al._ 2016 Ispolatov, I., Madhok, V. & Doebeli, M. (2016).  Individual-based models for adaptive diversification in high-dimensional phenotype spaces.  _Journal of Theoretical Biology_ , 390, 97–105.
  * Doebeli & Ispolatov 2017 Doebeli, M. & Ispolatov, I. (2017).  Diversity and coevolutionary dynamics in high-dimensional phenotype spaces.  _The American Naturalist_ , 189, 105–120.
  * Ispolatov _et al._ 2019 Ispolatov, I., Alekseeva, E. & Doebeli, M. (2019).  Competition-driven evolution of organismal complexity.  _PLoS Computational Biology_ , 15, e1007388.
  * Maynard Smith & Slatkin 1973 Maynard Smith, J. & Slatkin, M. (1973).  The stability of predator-prey systems.  _Ecology_ , 54, 384–391.
  * Bellows 1981 Bellows, T. (1981).  The descriptive properties of some models for density dependence.  _Journal of Animal Ecology_ , pp. 139–156.
  * Doebeli 1995 Doebeli, M. (1995).  Dispersal and dynamics.  _Theoretical Population Biology_ , 47, 82–106.
  * Yodzis 1989 Yodzis, P. (1989).  _Introduction to Theoretical Ecology_.  Harper & Row, New York.
  * Brännström & Sumpter 2005 Brännström, A. & Sumpter, D.J. (2005).  The role of competition and clustering in population dynamics.  _Proceedings of the Royal Society of London B_ , 272, 2065–2072.
  * Ricker 1954 Ricker, W. (1954).  Stock and recruitment.  _Journal of the Fisheries Board of Canada_ 11: 559–623.
  * Gyllenberg & Meszéna 2005 Gyllenberg, M. & Meszéna, G. (2005).  On the impossibility of coexistence of infinitely many strategies.  _Journal of Mathematical Biology_ , 50, 133–160.
  * Dieckmann & Doebeli 1999 Dieckmann, U. & Doebeli, M. (1999).  On the origin of species by sympatric speciation.  _Nature_ , 400, 354–357.
  * Hernández-García _et al_ 2009 Hernández-García, E., López, C., Pigolotti, S. & Andersen, K. H. (2009)  Species competition: coexistence, exclusion and clustering.  _Phil. Trans. Roy. Soc. A_ 367, 3183–3195.
  * Scheffer & van Nes 2006 Scheffer, M. & van Nes, E. H. (2006).  Self-organized similarity, the evolutionary emergence of groups of similar species.  _Proc. Nat. Acad. Sci. USA_ 103, 6230–6235.
  * Vandermeer 2006 Vandermeer, J. (2006).  Oscillating populations and biodiversity maintenance.  _Bioscience_ , 56, 967–975.
  * Shoresh _et al._ 2008 Shoresh, N., Hegreness, M. & Kishony, R. (2008).  Evolution exacerbates the paradox of the plankton.  _Proceedings of the National Academy of Sciences USA_ , 105, 12365–12369.
  * Parvinen _et al._ 2020 Parvinen, K., Metz, J.A. & Dieckmann, U. (2020).  Environmental dimensionality determines species coexistence.  _Journal of Theoretical Biology_ , 110280.
  * Gavina _et al._ 2018 Gavina, M. K., Tahara, T., Tainaka, K., Ito, H., Morita, S., Ichinose, G., Okabe, T., Togashi, T., Nagatani, T., % Yoshimura, J. (2018).  Multi-species coexistence in Lotka-Volterra competitive systems with crowding effects.  _Scientific Reports_ 8, 1–8.
  * Xue & Goldenfeld 2017 Xue, C. & Goldenfeld, N. (2017).  Coevolution maintains diversity in the Stochastic ”Kill the Winner” Model.  _Phys. Rev. Lett._ 119, 268101.
  * Benincà _et al._ 2009 Benincà, E., Jöhnk, K. D., Heerkloss, R., Huisman, J.(2009).  Coupled predator–prey oscillations in a chaotic food web.  _Ecol. Lett._ , 12, 1367–1378.
  * Lehahn _et al._ 2014 Lehahn, Y., Koren, I., Schatz, D., Frada, M., Sheyn, U., Boss, E. _et al._ (2014).  Decoupling physical from biological processes to assess the impact of viruses on a mesoscale algal bloom.  _Current Biology_ , 24, 2041–2046.
  * Uthicke _et al._ 2009 Uthicke, S., Schaffelke, B. & Byrne, M. (2009).  A boom–bust phylum? Ecological and evolutionary consequences of density variations in echinoderms.  _Ecological Monographs_ , 79, 3–24.
  * Behrenfeld _et al._ 2017 Behrenfeld, M.J., Hu, Y., O’Malley, R.T., Boss, E.S., Hostetler, C.A., Siegel, D.A. _et al._ (2017).  Annual boom–bust cycles of polar phytoplankton biomass revealed by space-based lidar.  _Nature Geoscience_ , 10, 118–122.
  * Mitchell _et al._ 2019 Mitchell, E. G., Harris, S., Kenchington, C. G., Vixseboxse, P., Roberts, L., Clark, C., Dennis, A., Liu, A. G., & Wilby, P. R. (2019).  The importance of neutral over niche processes in structuring Ediacaran early animal communities.  _Ecol. Lett._ 22, 2028–2038.
  * Louca _et al._ 2016 Louca, S. Parfrey, L. W & Doebeli, M. (2016).  Decoupling function and taxonomy in the global ocean microbiome.  _Science_ , 353, 1272–1277.
  * Mirollo & Strogatz 1990 Mirollo, R.E. & Strogatz, S.H. (1990).  Synchronization of pulse-coupled biological oscillators.  _SIAM Journal on Applied Mathematics_ , 50, 1645–1662.