Spatiotemporal ecological chaos enables gradual evolutionary diversification without niches or tradeoffs

  1. Aditya Mahadevan
  2. Michael T Pearce
  3. Daniel S Fisher  Is a corresponding author
  1. Department of Physics, Stanford University, United States
  2. Department of Applied Physics, Stanford University, United States


Ecological and evolutionary dynamics are intrinsically entwined. On short timescales, ecological interactions determine the fate and impact of new mutants, while on longer timescales evolution shapes the entire community. Here, we study the evolution of large numbers of closely related strains with generalized Lotka Volterra interactions but no niche structure. Host-pathogen-like interactions drive the community into a spatiotemporally chaotic state characterized by continual, spatially-local, blooms and busts. Upon the slow serial introduction of new strains, the community diversifies indefinitely, accommodating an arbitrarily large number of strains in spite of the absence of stabilizing niche interactions. The diversifying phase persists — albeit with gradually slowing diversification — in the presence of general, nonspecific, fitness differences between strains, which break the assumption of tradeoffs inherent in much previous work. Building on a dynamical-mean field-theory analysis of the ecological dynamics, an approximate effective model captures the evolution of the diversity and distributions of key properties. This work establishes a potential scenario for understanding how the interplay between evolution and ecology — in particular coevolution of a bacterial and a generalist phage species — could give rise to the extensive fine-scale diversity that is ubiquitous in the microbial world.

Editor's evaluation

This important study explores the question of “what gives rise to diversity in ecological settings?”. By considering the interplay between ecology and evolution, this study proposes a scenario of spatiotemporal chaos, in which interactions between strains drive large changes in the relative abundances of strains. The presented theoretical approach is compelling and goes beyond the current state of the art. This innovative theoretical work is of broad interest to the field of ecology and evolution.


A remarkable discovery of the DNA sequencing revolution is the vast diversity of microbes (Biller et al., 2015; Acinas et al., 2004; Kashtan et al., 2017; Rosen et al., 2018). Increasingly it has become clear that this diversity extends far below the level of conventionally defined species to finer and finer genetic scales (Bonilla-Rosso et al., 2020; Kashtan et al., 2014; Rosen et al., 2015; Tikhonov et al., 2015), and in some cases, a great multitude of strains coexist and compete in the same spatial location. Why doesn’t “survival of the fittest” drive almost all strains extinct, at least locally? Traditional explanations invoke the existence of a great many spatial or functional niches which limit competition between strains, down to “micro-niches” involving finer differences. However, especially for bacteria in relatively simple environments such as the marine cyanobacterium Prochlorococcus (Kashtan et al., 2014), does it make sense to postulate nano- or pico-niches, ad absurdum? Or is a statistical description of the small subtle differences more appropriate? Community ecology models with many similar strains competing for a mixture of resources have been much studied, but in their simplest manifestations the maximum number of coexisting strains is limited by the number of chemicals via which they interact, which in effect create a series of niches, each of which can be occupied by at most one strain (Chesson, 1990). Perfect “tradeoffs” are sometimes invoked to enable higher diversity (Posfai et al., 2017; Beardmore et al., 2011; Erez et al., 2020) but even tiny differences will destroy this coexistence (Caetano et al., 2021).

An alternative to the multi-niche scenario is the neutral theory of ecology which postulates that species are similar enough that they are somehow ecologically equivalent, with their population dynamics dominated by stochastic births, deaths, and migration. The predictions of this theory for abundance and spatial distributions are intriguingly similar to some data (Volkov et al., 2003; Volkov et al., 2007). However for microbes with short generation times and huge populations without tight bottlenecks, the neutral scenario is not viable: Even if the differences between strains could be neglected over the long times for which they have coexisted, the dynamics from stochastic fluctuations are far too slow. Instead, rapid population dynamics with large changes of relative abundance are often observed (Ignacio-Espinoza et al., 2020; Martin-Platero et al., 2018). “Selection”, in the broad sense of differential population growth rates, is clearly involved. Thus if a highly diverse population appears “neutral” in some respects (including close-to-perfect tradeoffs) this must emerge from the complex ecological and evolutionary dynamics; it should not be assumed.

It is often said that pathogens promote diversity (Bever et al., 2015; Rodriguez-Valera et al., 2009; Thingstad, 2000; Thingstad et al., 2014). However, there is thus far little understanding of how or under what circumstances ongoing coevolution of hosts and pathogens could cause and sustain extensive coexisting within-species diversity. Understanding this process theoretically is a long-term goal, towards which the present work is a step. To make progress, we need to distill the general phenomenon of fine-scale diversity to its most basic, and endeavor to develop potential scenarios in which evolution, coupled with ecology, might play out. For closely related strains, there is no compelling reason why interactions with siblings should be much stronger than those between distant cousins. Thus, we ask: Without assuming niche-like interactions, perfect tradeoffs, or spatial gradients, can a highly diverse collection of closely related strains stably and robustly coexist? If so, can such a highly diverse “phase” evolve and continue to evolve and diversify? If the evolution is fast, some amount of diversity will always exist (although the common ancestors of the population at any time may be recent and few). Thus we consider the most difficult regime for diversity: when the evolutionary dynamics are much slower than the ecological and spatial population dynamics.

In recent work (Pearce et al., 2020), referred to henceforth as PAF, we developed a new scenario for the coexistence of multiple closely related strains that are assembled all together into a community, leaving aside the question of their past or future evolution (or even how the community is assembled). In this scenario, we explored a particular key feature of models of many similar strains: the nature of interactions between pairs of strains. It is known that competition for resources in a well-mixed environment leads to positive correlations: if more A individuals are worse for B, then more B are worse for A. We consider the opposite case where the interactions are anticorrelated. This can arise if the competition is one-on-one: if A beats B, then B loses to A. A compelling biological motivation for anticorrelated interactions arises from a different scenario: a spectrum of generalist phage strains that prey, with varying efficacies, on a spectrum of bacterial strains. If a particular phage strain, a, does better than average against a particular bacterial strain, b, then more b individuals are better for a, and more a are worse for b, leading to anticorrelated interactions. While we are particularly interested in coevolving bacteria-phage diversity, to build up an understanding of the complex eco-evolutionary dynamics, we focus in this paper on simpler models that — as we have shown in PAF — capture many of the key features.

Host-pathogen, and other anticorrelated interactions, give rise to “kill the winner” ecological dynamics (Thingstad, 2000). If a strain rises to high abundance, other strains that do well against it will bloom and drive down the abundance of the first, and the process repeats. With many strains that do not have their own niches, this leads to wilder and wilder chaotic variations of abundances, soon driving most types extinct. In PAF we showed that rudimentary spatial structure — a large set of I islands with a low migration rate between all pairs of islands — can maintain much diversity without a commonly-invoked mainland (MacArthur and Wilson, 1967). In this spatial model, many strains go globally extinct, but a large fraction persists indefinitely in a spatiotemporally chaotic phase (hereafter STC). Crucially, the chaotic dynamics desynchronize across the islands allowing strains that go extinct locally to be repopulated from other islands. This mechanism is a manifestation of the “spatial storage effect” (Chesson, 2000). On each island, each persistent strain occasionally blooms up to high abundance and subsequently crashes (Figure 1). While it is at low abundance, dispersal from blooms on other islands rescues the strain from local extinction until conditions are favorable and its population blooms again, sends out migrants, and crashes. This STC is very robust: strains either go extinct rapidly, or persist globally for times that are exponentially long in the number of islands.

Dynamics of strain abundances on a single island in the spatiotemporally chaotic state (STC).

A subset of strains is plotted. Each persistent strain occasionally blooms up to high abundance and between blooms its abundance is sustained above a migration floor (here 10-7) set by migration from other islands, although a few marginal strains fluctuate below this threshold. (A) An example of the evolutionary process: at the beginning of an epoch (vertical dashed line), a new (here unrelated) strain, (black), is introduced at intermediate abundance. This new strain establishes and persists, causing two (red) strains, which persisted in the previous epoch, to go globally extinct by the end of this epoch. (B) Strains that would go extinct on a single island, can persist, and invade from low abundance, due to migration. The purple strain successfully invades. But at the vertical dashed line, migration is turned off for the purple strain only, and it proceeds to go extinct with average exponential decay rate given by its negative bias, schematically indicated by the dashed black line (with an extended range of log-abundance shown).

Complementary work (Roy et al., 2020; Roy et al., 2019) suggests the generality of the STC beyond anticorrelated pairwise interactions, although in the Lotka-Volterra models studied in these works of Roy et al., the diversity is limited by the strength of self-interactions, which also limits the diversity of stable communities. Indeed, much previous work has focused on ecological dynamics that reach a stable state, where diversity is limited by strength of niche interactions compared to inter-species interactions (Bunin, 2017). Here we approach evolutionary dynamics in similar generalized Lotka-Volterra models, but from the opposite starting point: all interactions are of comparable magnitude which makes the effects of self-interactions negligible compared to the effects of the total interactions from all other strains. Then there is no large stable community, and the diversity is maintained by spatial structure and chaotic dynamics.

With anticorrelated interactions, arbitrarily large numbers of strains can coexist in the STC even when spatial mixing — and hence competition — occur on timescales comparable to those of the local ecological dynamics. However, if the strains differ somewhat in their overall growth rate, or other ways that make some generally better, these advantages can limit the diversity of the community. A natural assumption is that, having all survived on evolutionary timescales, the persistent strains will be similar enough that such differences are very small. But this assumption — and even more so assumptions of close-to-rigid tradeoffs, (Posfai et al., 2017; Amicone and Gordo, 2021; Farahpour et al., 2018) — should surely be questioned. Such features must emerge from the evolution rather than being assumed.

Many theoretical (and some experimental) analyses, have, like our prior work, focused on ecological communities that are assembled without conditioning on their evolutionary histories: a number of species (or strains) is brought together, and the resulting community consists of the species that do not go extinct (Bunin, 2016; Bunin, 2017; Serván and Allesina, 2021; Friedman et al., 2017; Goldford et al., 2018; Hu et al., 2021). Although this is an important starting point, it is essential to incorporate evolution to understand how the processes of mutation, inheritance, selection, and extinction could give rise to highly diverse communities.

Previous theoretical work has shown that diverse communities in certain consumer-resource models are destabilized by evolution (Shoresh et al., 2008), at odds with the highly diverse continuously evolving microbial populations in nature. Others have focused on eco-evolutionary dynamics when the mutation rate is high enough to sustain diversity: in this case the common ancestor of coexisting strains is recent and extensive diversity over a wide range of genetic divergence does not have time to evolve (Xue and Goldenfeld, 2017). Yet others have shown that when niches in phenotype space are assumed, boom-bust dynamics can result in the evolution of higher diversity than occurs in stable equilibrium (Doebeli et al., 2021). But overall there is no clear consensus on whether evolution tends to destabilize or to increase diversity in ecologically interacting communities — indeed the answer to this question is likely context-dependent — though observations of the natural world suggest that evolution often results in increased diversity. Here, we investigate evolution starting from a state with spatiotemporally chaotic ecological dynamics as studied in PAF, where niches are absent and the diversity — at least initially — is stabilized by the interplay between endogenous ecological dynamical fluctuations and migration.

We are interested in understanding diversity that has existed for a very broad spectrum of evolutionary timescales, far longer than ecological or spatial mixing timescales. We thus study the extreme limit where the mutation rate is small enough that the ecological and migratory dynamics reach steady state before the introduction of each new strain. This quasistatic limit of evolution is the “hardest” for diversification. In addition, and in contrast to some previous work (Tikhonov and Monasson, 2018), we assume that global extinctions of a strain are permanent: an extinct organism cannot be resurrected even if conditions later become favorable for it. Focusing on the STC phase, we endeavor to answer: Can a highly diverse STC phase evolve? Under what conditions? Can this phase continue to diversify? Is the diversity stable to general fitness mutations that are not artificially constrained by assumptions of tradeoffs? How do statistical properties of the community change during evolution?

Summary of main results

We first summarize the main results of this work, which concern the behavior of the STC phase under serial invasion of new strains. A parent strain in the community occasionally gives rise to a mutant strain whose properties are correlated with those of its parent with correlations parametrized by ρ[0,1). Many of the behaviors are similar across this range, from independent invaders with ρ=0 to small-effect mutations with ρ close to 1. In all cases, extinctions are irreversible. Key properties of interest are the number of extant strains in the ecosystem, L, the number of successful invasions Z, and how these evolve with evolutionary time, parametrized by the number of attempted invasions, T.

We find that for a wide range of ρ (and expect for any ρ<1) the STC can enter a steadily diversifying state wherein the number of successful invasions, Z, and the number of coexisting strains, L, both increase linearly with T on average, with only small fluctuations when L is large. Whether diversification occurs, and the rate of the diversification if it does, depends on various parameters, but it is robust over a range of the parameters. If initially the ecosystem has only a modest number of strains, the evolutionary dynamics tend to cause the diversity to crash, after which is it extremely unlikely to transition into the diversifying phase. However, if the initial ecosystem is sufficiently diverse, it is highly likely to diversify further.

We then study the effects of general fitness differences that augment the average growth rate of a strain by an intrinsic amount irrespective of its interactions with the other strains. Focusing on unrelated invaders, we show that a distribution of such general fitness differences (denoted by si for strain i) can either slow down, prevent, or reverse diversification. For distributions of the si whose tails decay faster than exponentially, the diversifying phase still exists, but with the diversification rate gradually slowing down: L increases only as a power of log(T). If the distribution of the si has a broader-than-exponential tail, the diversity decreases and crashes.

The key property of a strain, in terms of which one can understand its behavior, is its bias: defined as the rate at which its population would change when at low-abundance and without migration (Figure 1B). The crucial effect of migration in the STC is to stabilize many strains with negative bias which would have gone extinct without migration. Only if its bias is strongly negative will a strain go globally extinct. The bias of strain i, denoted by ξi, has an intrinsic contribution from its general fitness si, and an extrinsic contribution from the interactions of strain i with all other strains, which includes both an i-dependent part determined by the interactions with the other strains in the community, and an i-independent part that keeps the total population constant. The extrinsic part of the bias of each strain changes as the community evolves, but its intrinsic part says the same.

Building upon the theoretical understanding of the STC phase, developed in PAF, we first analyze evolution in the simplest case of unrelated invaders (ρ=0) with no general fitnesses (si=0). The bias of each strain undergoes a random walk on evolutionary timescales, and we find that for large communities, the number of strains changes at a steady rate. For a range of parameters, this diversification rate is positive, yielding a steadily diversifying phase with the distribution of biases scaling with 1/L, as observed in numerics. We then extend our analysis of the changing bias distribution to include the effects of general fitness differences. This yields predictions of how the rate of diversity increase (or decrease) depends on the distribution of the si, corroborating the behaviors found in simulations.

The distributions of biases and abundances in evolved communities differ subtly from those of the initial communities that were assembled all-at-once from unrelated strains. At early stages of the evolution, most of the close-to-marginal, low-abundance strains are pushed out by the perturbations caused by the invading strains. This extinction process causes the shape of the abundance distributions of assembled and evolved communities to differ at low abundances. Later, in the steadily diversifying state, the numbers of extinctions caused by each invader has a roughly exponential distribution, which is consistent with our theoretical expectations. In contrast to the qualitative (albeit modest) changes in abundance distributions, we find that evolution has only a small effect on the statistics of the interactions between strains.


The structure of this paper is as follows: ‘Models’ introduces the main model and its relation to previous work. ‘Results’ describes the phenomenology of an evolving community in the STC phase, studying the effects of correlated mutants, interaction statistics, and general fitness differences on the ecological diversification. Then ‘Analysis’ develops the theory and analysis that are needed to understand these phenomena. Building upon the dynamical mean field theory developed in PAF, we present an approximate framework, and more general scaling arguments, for understanding the evolutionary dynamics, and compare the predictions with simulations. Finally ‘Discussion’ raises additional questions and discusses possible extensions. Many of the details and further analyses are relegated to appendices.


We here define the model, discussing the various roles played by local deterministic population dynamics, demographic stochasticity, spatial migration and evolutionary dynamics. Our notation is summarized in Table 1.

Table 1
Definitions of commonly used quantities.
STCSpatiotemporally chaotic state
KNumber of strains put into the initial assembled community
VMatrix of pairwise strain interactions
γSymmetry parameter of the interaction matrix; E[VijVji]=δij+γ
INumber of islands
NPopulation size on each island, fixed to be constant
νi,αFractional abundance of strain i on island α
ν¯iTime (or space) average of strain i abundance
νfloorMigration floor mν¯L/M: ~ lower range of local abundances
ΥαLagrange multiplier maintaining iνi,α=1; Υα=iνi,α(si+jVijνj,α)
siGeneral fitness of strain i
P(s)Probability distribution of the si
ΣCharacteristic scale of the s distribution
ψExponent characterizing tail of of P(s)exp[(s/Σ)ψ/ψ]
mMigration rate between islands
MRange of fluctuations in logν; M=log(1/m)
ρCorrelation between parent’s and mutant’s interactions with other strains
TEvolutionary time in epochs, equal to number of attempted invasions
ZNumber of successful invasions
LNumber of extant strains at any point in the evolution
L0Number of strains surviving in initial assembled community
UAverage diversification rate; U=dL/dT
s^Mean general fitness of extant strains;  s^=jνjsj
Σ^The scale of P(s) of extant strains; Σ^=[ddslogP(s)]1|s=s^
ξiBias of strain i, its growth rate at low abundance without migration; ξi=ζi+si-Υ¯, scales as  1/L
N(ξ)Mean abundance of a strain as a function of its bias
ξcCritical bias (negative in the STC) below which strains go extinct, scales as 1/L
ζiMean drive on strain i by other strains in its absence; ζi=jVijν¯ji
LEffective number of extant strains; L=1/jν¯j2; scales with  L
χiStatic response of strain i to perturbations; χi=dν¯idξi
XTotal static response;  X=iχi
ΞFragility of the community to perturbations;  Ξ=jχj2/(1-jχj2)

Ecological interactions

We first consider an assembled community of K unrelated strains, labelled by i=1,2,,K, with all possible pairwise interactions between them. A paradigmatic model for the ecological dynamics of the strain populations {ni} is the generalized Lotka-Volterra model (Goel et al., 1971), with each strain i having an intrinsic growth rate which is modulated by its interactions with all the other strains. These interactions are conveniently represented in a matrix W where Wij describes the effect of strain j on the growth rate of strain i. Since we are interested in closely related strains for which all interactions are similar, the total population will be roughly fixed at some N by the balance between the effects of positive intrinsic growth rate and negative competitive interactions. It is convenient to replace these large terms by a Lagrange multiplier Υ(t) that fixes the total population to ini=N, and work with fractional abundances, νini/N. This parameterization yields what are known as “replicator equations” (Chawanya and Tokita, 2002; Yoshino et al., 2008; Tokita, 2004).

Variations in intrinsic growth rates and net interactions on a strain can be combined to yield general fitness differences, {si}, between the strains. We parameterize the residual variations in interactions among the strains (after subtracting off Υ and si) by Vij. Since the Vij and si are sums and differences of similar magnitude terms, it is natural to approximate them as random variables with the hope that the model will yield behaviors that are robust to specific choices of their statistics: testing this assumption is one of the goals of this paper. For simplicity, we choose E[Vij]=0, and E[Vij2]=1 for ij — setting the overall ecological timescales — and choose the covariances to be zero except for, importantly, correlations between how i and j affect each other, defining E[VijVji]=γ. For convenience, we choose E[Vii2]=1+γ but this choice has negligible effect in large communities.

The parameter γ controls whether the interactions are mainly competitive (γ>0) or host-pathogen-like (γ<0), the latter being the focus of this work. We have shown in PAF that random interaction matrices with such anticorrelations behave very similarly to host-pathogen models with the appropriate block sub-matrix structure, as discussed further in ‘Bacteria-phage interactions and coevolution’.

Ecological dynamics

We study the simplest model with spatial structure: a large number, I, of identical islands (or demes) with interactions only within each island and migration between all pairs of islands. With Greek indices labeling islands, the dynamics of the abundances obey

(1) dνi,αdt=νi,α(si+j=1KVijνj,αΥα(t))+β=1I(mαβνi,βmβανi,α).

with mαβ the migration rate (per individual) from island β to island α and the local Lagrange multiplier, Υα=iνi,α(si+jVijνj,α), keeping the total population on each island fixed at N, (i.e. iνi,α=1 for each island). Here, we focus on the spatial mean field limit in which the migration rate is the same, given by m/I, between every pair of islands. The total migration of strain i into and out of island α is then simply m(ν¯i-νi,α) with ν¯i(t) the average of νi(t) across islands. As the number of islands becomes large, in steady state, each ν¯i(t) becomes constant in time — some being zero corresponding to global extinction. In the STC, the dynamics are asynchronous across islands and ergodicity implies the spatial average, ν¯i, is equal to the time-averaged abundance of strain i on a single island; this is a crucial self-consistency condition. The magnitude, m of the migration rate, is also of fundamental importance. If m is too small the migration is too rare to repopulate islands after local extinctions. If m is too large and the local dynamics is chaotic, the chaos will synchronize across the islands and the total population of each strain will fluctuate wildly, rapidly driving most strains extinct. We will focus on the wide intermediate m regime, which spans several orders of magnitude when K is large [PAF].

With large populations on each island, demographic fluctuations have little effect on the dynamics. Even when the local population of a strain is small, if it has positive growth rate, fluctuations will not matter much, while if it has negative growth rate it will go deterministically extinct which occurs when the fractional abundance drops below the extinction threshold of 1/N indicated in Figure 1A by the horizontal purple line. The value of the extinction threshold does not much affect the behavior as long as it is much below the lower limit of the abundance caused by migration — which we term the migration floor: νfloormν¯. For strains near local extinction (when the fractional abundance is close to 1/N) demographic fluctuations are potentially important. But with Nmν¯ very large, local extinctions for viable strains will be rare: thus we model the population dynamics as fully deterministic. If the fractional abundance on an island drops below 1/N, it is set equal to zero. Global extinction occurs when a strain’s bias becomes too negative, which results in it going below the extinction threshold everywhere. The choice of N does not matter much as long as mN1, to which we restrict consideration. Related details of numerical implementation are discussed in ‘Appendix 2’. In ‘Spatial structure and dynamics’, we comment on the effects of local extinctions in the context of real spatial dynamics.

The key properties of the STC phase [PAF] are chaotic coexistence of strains, desynchronized across islands, with the local abundances fluctuating over a range in logν of Mlog(1/m), which is quite wide for the typical m=10-5 that we use in simulations. Some strains go globally extinct but each persistent strain on each island occasionally has a bloom up to high abundance νM/L. These localized blooms are crucial for stabilizing a strain, as they dominate the migration to other islands needed to recover from local extinctions or near-extinctions. On a single island, at any given moment, O(L/M) strains are at high abundance. A snapshot of the abundances on a single island shows the strains distributed roughly uniformly in log(ν) down to the migration floor, only occasionally fluctuating substantially lower (Figure 1).

Evolutionary dynamics

The evolutionary process we model is much slower than the ecological and migratory dynamics. Simulations are divided into long epochs, with new strains added only at the end of an epoch. The epochs are chosen long enough that the ecological and migratory dynamics have reached a steady state, with some fraction of the strains having gone permanently extinct globally, leaving L persistent strains. A single new strain is then introduced and the process repeated.

The new strain, generically labeled A, is parameterized by its interactions with all other strains in the community, given by VAj and VjA, and its general fitness, sA. In the simplest case, a new strain is unrelated to extant (or extinct) strains. More generally, mutant strains, labelled M, are characterized by their degree of correlation, ρ[0,1), with a parent strain P chosen from the existing community with probability proportional to its mean abundance ν¯P. These correlations are realized such that Corr[VjM,VjP]=Corr[VMj,VPj]=ρ for jM,P. The detailed choices for j=M,P are given in ‘Appendix 2’. The general fitness, sM, can also be correlated with sP. Unrelated invaders are equivalent to ρ=0 and hence have no parent.

The actual process of invasion from low abundance on one island is complicated, and often leads to failure. To avoid a proliferation of such failed invasions, we instead assess whether the invader could successfully invade and persist if it were lucky initially. To do this, we set the mutant’s abundance to 1/L on all the islands at the same time (and proportionately decrease the abundances of the other strains to maintain iνi,α=1).


There are multiple timescales involved in the dynamics: these are discussed more fully in ‘Appendix 1’. The basic timescale for differential growth or decay of strains is set by the magnitude of the interactions and the number, L, of extant strains. The extant strains have average abundances of order 1/L, so the average total interaction on a strain, i, is the sum of L random terms, each of order V/L for typical interaction strength V. With the V having variance unity, the average net interactions of other strains on strain i is roughly its mean drive, ζi, which is of order 1/L, implying that the timescale for systematic population growth or decay is of order L. The mean drive is defined more precisely in Results. When there are general fitness differences, si, these also contribute to variations in average growth rates. The variations in the si within the community have substantial effects over a time 1/σs with σs roughly the width of the si distribution of the extant strains. Together, the mean drive and general fitness of a strain determine its crucial property: the bias ξi=ζi+si-Υ, with angular brackets denoting a time average. As introduced earlier, the bias of strain i is its average growth rate at low abundance in the absence of migration (Figure 1B). As we shall see, the size of the community is limited by the condition that the inter-strain variation in general fitness is no larger than variation in average drive from interactions. This means that the biases are of order 1/L.

The local population of each strain undergoes wild fluctuations over a logarithmic range M=log(1/m) which is quite large. During blooms, the instantaneous growth and decay rates of local populations are substantially larger than the systematic biases (‘Appendix 1’) and change rapidly from growth to decay as seen in Figure 1. The time for abundances to fluctuate from large to small — the duration of blooms — is of order ML with systematic and fluctuation contributions comparable.

An important timescale for studying slow evolution is the time to reach the STC steady state: the ecological relaxation time. This is determined by the strains that are just barely going extinct and is of order ML for an evolved community, as discussed in ‘Continual assembly and diversification’. We have chosen the evolutionary timescale to be much longer than all the other important timescales. Thus, each epoch between the addition of invaders is chosen to be several times the ecological relaxation time, typically 3ML, and we show in ‘Appendix 2’ that increasing this epoch length by a factor of 10 makes little difference in the diversification dynamics.


The STC is robust, with strains persisting for times that are exponentially long in the number of islands. However, evolutionary perturbations caused by an invading strain can drive strains deterministically extinct. This process can be understood in terms of the biases of the strains.

The bias of a strain, is determined by the community in its absence. It can be written precisely as ξi=si+ζi-Υi with the mean drive ζi=Vijνji, where the notation νji and Υi denote the time-averaged abundance of strain j and average Lagrange multiplier in the absence of strain i. The Υ is not much changed by the absence of the one strain, but the abundances of the other strains are affected in small but collectively essential ways by whether or not strain i is present, as discussed in ‘Dynamical mean field theory’. This negative feedback — proportional to γ — is what stabilizes strains whose abundance would otherwise keep growing.

With many strains participating in the chaos on each island, and desynchronization across islands, we expect the chaos to be ergodic, so that the i time averages and spatial averages (across islands) of all quantities are equal in the STC steady state. Therefore we will use spatial average notation ν¯ instead of time average notation ν, except when conceptually the latter is clearer. In practice, the I=40 islands used in numerics are enough that the persistence times of almost all surviving strains are very long and averages across islands of the more important quantities do not fluctuate much in steady state.

A crucial feature of the STC phase is that strains with somewhat negative bias can persist due to migration between desynchronized islands (Figure 1A). This stabilization is enabled by a nontrivial feature of the STC phase: during a bloom, the systematic changes in logνi caused by the bias are comparable to the cumulative stochastic growth and decay caused by the endogenous fluctuations — the zigs and zags in the dynamics of logνi (Figure 1). This is a manifestation of the system “self-tuning” to a special self-consistently chaotic state [PAF].

Despite the possibility of rescue from extinction via rare blooms, there is a critical negative bias, ξc, (sharp for large L and large I) below which strains no longer persist even as I. For strains with ξ below ξc (which depends on the parameters and the number of strains), blooms up to high abundance are not frequent enough to repopulate local extinctions and deterministic global extinction ensues. For large I and large L, strains with ξ<ξc go extinct, while strains with ξ>ξc persist indefinitely. Finite L and finite I effects, together with the finite time for each epoch, will round out the sharpness of the borderline between persistent and extinct. However the marginal strains involved have little effect on others and whether or not they persist does not much matter for the current epoch: we are interested in deterministic extinction caused by the introduction of new strains. Therefore, we need to study how the distribution of biases in the ecosystem evolves.

Continual assembly and diversification

The evolutionary process we study starts from an assembled collection of K unrelated strains. After the ecological and migratory dynamics have reached steady state, some of the strains will persist: we call the size of this initial persistent community L0. The K-L0 strains that have gone globally extinct are permanently removed.

When a new strain is introduced into the ecosystem, if it successfully invades it perturbs the biases of the extant strains, and can trigger extinctions of some of them by shifting their bias below ξc (Figure 1A) by an amount of order 1/L. We study the slowly evolving regime in which the ecosystem dynamics reach steady state between each introduction of a new strain — this takes time of order ML. The number of persistent strains and number of successful invasions as a function of the number of attempted invasions, L(T) and Z(T) respectively, are of fundamental interest.

We first describe the evolutionary dynamics when the general fitness differences between the strains can be neglected. For γ=-0.8 and unrelated invaders (ρ=0), multiple simulation runs starting with different sets of K=50 initial strains reveal that around one fifth of the replicates enter a steadily diversifying regime in which L increases roughly linearly with the number of attempted invasions, at a rate of around 0.25 per attempt. The remaining replicates crash down to only a few persistent strains. Subsequent invasions can cause L(T) to increase somewhat, but it quickly crashes back down and the community does not steadily diversify (Figure 2A). The low diversity regime that occurs after a crash (or with a very small initial community) is discussed further in ‘Appendix 3’.

Evolution of number of strains without general fitness differences.

(A) With γ=-0.8, m=10-5, and initial number of strains K=50, under serial invasion of unrelated strains most initial communities (red) crash and fail to recover, while others (about 20%, blue) continually diversify. Once the communities are large, around 80% of further invasions are successful and the mean number of extinctions per successful invasion is 0.7 (Appendix 3—figure 1) so that on average the number of strains in the community grows linearly with rate U0.25 per invasion attempt (dashed line). (B) Whether diversification occurs, and its rate if it does, depends on the symmetry parameter, γ, as seen here with K=400 and ρ=0. For γ close to -0.7, evolution reduces the diversity. For less negative γ, the STC breaks down and the diversity crashes immediately. For more negative γ, steady diversification occurs, fastest here with γ-0.8, though again slowing down as γ-1. (C) Evolving communities under successive introduction of mutants, each with correlation ρ with its parent (γ=-0.8). The diversification rate varies nonmonotonically with ρ, with fastest diversification for ρ0.8. There is a significant slowdown for ρ close to unity. Here K=50 and trajectories are shown conditional on not crashing, except for ρ=0.99, which renders the evolving community very susceptible to crashing from K=50. However the inset shows that even with ρ=0.99, it is possible to reach a diversifying regime starting from K=100.

The observations in Figure 2 illustrate one of the crucial findings of this work: spatiotemporally chaotic ecological dynamics can allow — but do not guarantee — gradual strain-level diversification up to arbitrarily high number of strains. The behavior depends on the symmetry parameter γ, which must be substantially negative for the STC to exist. Figure 2B shows that the average rate of diversification, UdL/dT, is nonmonotonic in γ, with slow diversification close to γ=-1 (its lower limit). As γ becomes less negative the rate of diversification increases at first. However for γ even less negative, the STC still supports chaotic coexistence of many strains (since L0 is still large), but the diversity decreases under evolutionary dynamics. The community diversifies most rapidly for γ-0.8. As we are interested in what can happen with various other additional features, we chose γ=-0.8 and m=10-5 for all further simulations as shorter runs are needed near these values. We expect that the qualitative conclusions will be similar for a range of γ and m around these.

Evolution with correlated mutants

In addition to studying independent invasions, we study evolution via mutations of existing strains. At the start of each epoch, a parent to mutate is chosen with probability proportional to its mean abundance. The interactions of the mutant with other strains are drawn from the same marginal distribution as the original interactions, but with correlation ρ with the interactions of the parent (Evolutionary dynamics). The direct interactions between the parent and mutant have to be chosen separately as specified in ‘Appendix 2’ but, as they only account for a small fraction of the total abundance in diverse communities, the specific choice is not important. As a function of ρ (with γ=-0.8), the rate of diversification is nonmonotonic being fastest for ρ0.8, and only weakly varying for smaller ρ (Figure 2C). As ρ nears 1, the mutant and parent are more similar, and it becomes harder for them to coexist, since any difference between them is likely to result in a systematic change in their relative abundance, eventually driving one of them to extinction (see ‘Appendix 9’). Since L can only increase when both the mutant and parent coexist, increasing ρ slows the rate of diversification but L(T) still increases linearly.

This observation implies that, for a large range of ρ, despite not enforcing any precise constraints or perfect tradeoffs, strains that would outcompete all extant members of the community are too rare to emerge and reduce diversity. Therefore, we conjecture that a continually diversifying phase exists even for ρ arbitrarily close to 1.

Evolutionary dynamics with general fitness differences

So far we have observed that when mutants or invaders differ only by their interactions with each other, there is robust and rapid diversification, provided that the initial diversity in the STC phase is high enough. Now we include general fitness differences si between strains and show how these affect the evolutionary dynamics.

Exponential distribution of si

We first analyze the simplest case: exponentially distributed selective differences with scale Σ and probability density P(s)=1Σes/ΣΘ(s) where Θ(s) is the Heaviside step function. We consider the evolutionary dynamics in the case of unrelated invaders. As a community evolves, the distribution of the si of the community will change, and we are particularly interested in the dynamics of the population-weighted mean s^=jν¯jsj.

The width of the distribution of si, here Σ, plays a controlling role. If one strain has a substantially higher growth rate than all other strains, it will outcompete them, driving many extinct. Thus a broad distribution of si is likely inconsistent with a diverse community. We therefore focus on narrow distributions: that is small Σ. The typical magnitude of the drive of a strain is of order 1/L; therefore, when Σ is much smaller than this, it will not matter much. On the other hand, if Σ were to be much larger than 1/L, the differences in the si would dominate over the drives and only the strains with the highest and quite similar si would survive. Thus L1/Σ2 seems inconsistent. Even for the initial community with L0 strains, we expect that L0 cannot be larger than order 1/Σ2 (although it can be much smaller if K1/Σ2). Indeed, in ‘Appendix 6’ we show that Σ sets the initial persistent community size, L0, in a particular limit of the model where γ=-1.

A natural conjecture is that for small Σ with an exponential distribution, steady diversification can occur until the breadth of the s distribution becomes important — when L1/Σ2 — and after that L will saturate, as seen in Figure 3. Thereafter, s^ will grow and invasions of unrelated strains are less and less likely to be successful — in ‘Diversification rate with a distribution of general fitnesses’ we show that the number of successful invasions Z, increases as log(T). However, successful invasions will on average drive exactly one other strain extinct. Figure 3 illustrates this behavior, including the large initial drop from K to L0 when K1/Σ2, the Σ-dependence of the steady-state L, and the linear increase of s^ with number of successful invasions.

Effects of exponentially distributed general fitnesses, si, on community evolution.

Here K=250 initial strains and P(s)es/Σ, with various Σ. (A) Community size, L, as a function of evolutionary time T (the number of attempted invasions) approaching an evolutionary steady state with LΣ-2 at long times. The dashed lines indicate 0.7×Σ-2, which captures the predicted scaling between the steady-state L and Σ. Data are averaged over 50 runs, conditional on not crashing, with the shaded region showing the standard error. Only single runs are shown for Σ=0.1 and 0.02, the former caused crashing and the latter saturation beyond the range of the simulations. For a narrow distribution of the general fitnesses (L0Σ-2), L increases linearly before saturating. For larger Σ with many initial strains, immediate extinctions drive L down to L01/Σ2 (‘Appendix 6’). (B) The average fitness of the community, s^=iν¯isi, grows linearly in the number of successful invasions, Z, with ds^/dZΣ3: the dashed lines have slope Σ2/6, indicating this expected scaling relationship. The inset shows the rate of successful invasions slowing down with attempted invasions, as it gets harder to draw a general fitness that is sufficiently far into the tail of P(s).

More general distributions of si

Building on an understanding of the case of exponentially distributed s, we consider a more general family of distributions, motivated by the expectation that the tail of the s distribution is particularly important for evolution:

(2) P(s)exp[1ψ(sΣ)ψ]Θ(s).

Anomalously small s strains are very unlikely to successfully invade, so the sharp cutoff at the lower end does not matter. Although we consider only positive si, all si can be shifted by a constant without affecting the dynamics because this constant gets absorbed into Υ(t).

As we will analyze in ‘Diversification rate with a distribution of general fitnesses’, the evolution of diversity is seen to depend crucially on ψ. If the tail of the s distribution falls off faster than a simple exponential, ψ>1, the community continually diversifies, albeit more and more slowly with L(T) increasing only as a power of log(T). Concomitantly, the mean s of the community, s^, gradually increases. But if the s distribution decays slower than a simple exponential, ψ<1, the diversity decreases (after an initial increase if Σ is sufficiently small) and eventually crashes. In the marginal case of a simple exponential tail, ψ=1, as seen above, the diversity saturates and fluctuates around a steady state value while the mean s^ increases linearly with the number of successful invasions. Therefore we conclude that for the evolutionary process in our models to continually generate higher diversity, the distribution of general fitnesses must decay sufficiently rapidly. Such rapid decrease of the distribution of available beneficial mutants with ongoing evolution roughly corresponds to “diminishing-returns epistasis”.

Mutants with correlated general fitnesses

What happens if — as one would expect — the invaders are mutants with general fitnesses correlated with their parents? With such mutants, it is possible for the evolution to proceed with less slowing down than for independent invaders. Indeed, with an exponential distribution of the si (ψ=1) analysis suggests that evolution proceeds at a constant rate, with both Z(T) and s^(T) growing linearly as in the absence of general fitness differences, but L still saturating. In ‘Appendix 7’ simulation results are shown for an exponential P(s) with correlations in both interactions and general fitnesses. The saturating value of L is quite similar in both the correlated and uncorrelated cases. However, for correlated mutants s^(T) pushes rapidly into the exponential tail — and surely toward the breakdown of the assumption of the existence of such large s mutations.

For P(s) that decays faster than exponentially, ψ>1, the behavior is more complicated. However, as discussed in ‘Appendix 7’, even with correlated mutants, evolution will eventually become very slow, as for uncorrelated invaders. With mutants instead of unrelated invaders, this is a direct example of the effects of diminishing-returns epistasis.


In this section we develop an approximate analytical theory of the evolutionary dynamics and provide heuristic understanding for most of the observed phenomena described above. The underlying basis is the dynamical mean field theory (DMFT) of the STC phase developed in PAF. This takes advantage of the large number of strains and the large number of islands in order to simplify the descriptions and analyses of the behaviors.

The natural quantities that characterize strains in the DMFT are their biases, {ξi}, and how these set their mean abundances, {ν¯i}. For a large randomly assembled or evolved community the mean abundances will be a function of the biases: ν¯iN(ξi), with the function N depending on the parameters, evolutionary history, and feedback from other strains. As shown in PAF, N(ξ) is linear for large argument and decays as ξ becomes negative, vanishing at ξ=ξc<0.

The relation between ν¯i, ξi and the total average force on strain i — from both direct effects and feedback — enables one to estimate the bias from the simulations (‘Appendix 2’). Armed with the DMFT description, we can understand how the biases of extant strains change over the course of invasions. We do this in detail for the simplest case — invasions of unrelated strains without general fitness differences — and show that the evolution causes L to change linearly with the number of invasions — decreasing or increasing depending on the parameters. A simple approximation to the evolution of the biases enables semi-quantitative results. Of particular importance is the result that in evolving communities, the density of biases vanishes linearly as ξξc. A corollary of this, as discussed in ‘Distribution of biases and number of extinctions’, is that the number of extinctions per successful invasion is roughly exponentially distributed. We then analyze the effects of general fitness differences, using our understanding of the exponential P(s) to generalize to other shapes of the tail of P(s), parametrized by ψ (Equation 2), and showing how the steepness of the tail affects the rate at which L increases or decreases.

Dynamical mean field theory

The DMFT approximation, which is exact in the limit of a large number of strains with random interactions between them, replaces the full statistical dynamics by the stochastic effects of the others on one chosen strain, with the statistical properties then determined self-consistently from the properties of the distributions over the strains. This approach was first used in the physics of disordered systems such as spin glasses (Sompolinsky and Zippelius, 1982), but has been applied to ecological dynamics in a number of subsequent works (Diederich and Opper, 1989; Opper and Diederich, 1992; Galla, 2006; Roy et al., 2019; Rieger, 1989; Yoshino et al., 2008).

When strain i has very low abundance, its effects on the others are very small and the forces of the others on it, jVijνj, are comprised of roughly independent random variables and thus act like gaussian noise with correlations C(t,t)=jνj(t)νj(t). However, when it rises to substantial abundance, it will weakly affect the other strains. Because of the correlation γ between Vij and Vji, these feedback terms add coherently, resulting in a contribution to growth rate of i of form γ0tR(t,t)νi(t)dt, with R a response function determined by feedback from the total impact on the community of the strain’s own past history (see ‘Appendix 5’). The DMFT allows one to recast the generalized Lotka Volterra equations as an effective single-strain problem, with self consistency conditions on the bias correlations and response function.

With the dynamical mean field understanding of an assembled STC phase in hand, we can proceed to describe the evolutionary process in terms of the distributions of properties of the extant and newly invading strains — in particular their biases and consequent mean abundances. The distribution of biases is perturbed by the introduction of new strains and this can push some of the extant biases below ξc, which itself depends on the bias distribution as modified by prior evolution, and on the number of extant strains, L. It is convenient to define an effective community size L=1/iν¯i2 which controls the variances of the mean-drive part of the bias; L scales with the actual L but discounts strains that are close to extinction.

Evolution without general fitness differences

We first analyze invasion of unrelated strains without general fitness differences. The bias of an attempted invader, labelled A, is given by ξA=jVAjν¯j-Υ¯, in terms of its interactions, {VAj} with the extant strains: for unrelated invaders, this is gaussian distributed with mean -Υ¯ and standard deviation of 1/L independent of correlations among the extant strains, though correlations in the existing community will affect the N(ξ) and hence L.

Strain A can successfully invade the community and persist if ξA>ξc. The probability of successful invasion is thus Φ[(ξc+Υ¯)L], with Φ the standard normal cumulative distribution function. In the initial assembled community, ξcL and Υ¯L are independent of L. We make the Ansatz that after a long period of evolution the distribution of extant biases, scaled by L, reaches a steady state — albeit a different state than the initial assembled community. Then the probability of successful invasion will become independent of L for large L. If the mean number of extinctions per successful invasion also reaches a steady state value which is less than unity, this explains the steady linear growth of L(T) seen in Figure 2.

With the one-by-one introduction of new strains, the bias of each extant strain undergoes some kind of random walk, and the strain goes extinct if its bias ventures below ξc. In Figure 4A, we show the evolutionary trajectories of the biases of 5 individual strains that started from similar initial values in a simulation where the community diversified from 50 to 500 strains. Extinctions are caused by ξ being pushed below ξc by an invading strain. For finite L, the sharpness of ξc will be smeared by an amount of order 1/L due to variability in the dynamic noise from strain to strain, which we have not explored.

Trajectories of biases of persistent strains (normalized by 1/L) under the influence of successive unrelated invaders with si=0.

(A) Bias trajectories of individual strains that invaded and persisted for a number of epochs. Extinctions (shown by a vertical line), occur when the bias goes below the critical bias, seen here to be around -2.5/L. The horizontal dashed line shows -Υ¯L. (B) Bias trajectories for all strains binned into groups by their starting value and averaged within bins for as long as the strains persist. Without conditioning on success, the biases of new invaders have mean -Υ¯ and standard deviation 1/L. However, conditioned on survival, the biases converge to and fluctuate around a larger value. Data in (A) are from a single simulation where the community diversifies from 50 to 500 strains, and in (B) data are pooled from 10 replicates of the same process.

To numerically investigate any systematic components of the random walk of biases, we average over a large number of strains, binning them according to their initial values normalized by 1/L. We observe a strong tendency of anomalously positive and negative biases to regress toward an intermediate value. In this plot, as evolution proceeds, the asymptotic average bias conditioned on survival is larger than -Υ¯L. (Figure 4B). This is likely due to conditioning on survival of the strains: those that persist for many epochs tend to have larger-than-average (but still negative) bias.

In ‘Appendix 8’, we carry out an analysis of the bias dynamics by approximating these by a Markov process in which the dynamics of the biases depend only on their current values. This analysis shows that when there is a successful invasion, each strain’s drive undergoes both a systematic and random change, consistent with our numerical results.

It is convenient to work with the mean drives, ν¯i (when si=0, this is just ξi+Υ). Both the systematic and the random changes in the drive are proportional to the average abundance of the invading strain, which is of order 1/L. The stochastic change is δζi±1/L, but the systematic change in the drive is smaller and depends on its current value: E[δζi|ζi]ζi/L=O(1/L3/2). In the Markovian approximation, there is a simple Langevin equation for the change of the drive of a strain due to invasions:

(3) dζidT=BζiL(T)+2D1L(T)ηi(T),

with ηi(T) approximately gaussian with mean zero and unit variance — an approximation that should be good if one coarse-grains over a substantial range of T (but with range much smaller than L). From our analysis in ‘Appendix 8’, we see that B and D are order-unity coefficients which respectively characterize the average and mean-squared response of the bias to the invasion of a new strain. Both are proportional to E[νA2] times the fragility, Ξ, of the extant community which is given by Ξ=jχj2/(1-jχj2) in terms of the individual susceptibilities of strains to changes in their biases, χi=dν¯i/dξi. This fragility characterizes the mean-square response of the system to a random perturbation applied simultaneously on all the strains — precisely the effect of a successful invasion. The Langevin equation for the drives must be supplemented by a boundary condition that if ξ=ζ-Υ goes below ξc, the strain disappears.

Analysis of the Langevin equation, which can be converted into a Fokker Planck equation for the distribution of the drives (‘Appendix 8’), shows that there is an eigenvalue-like condition which determines whether the diversification rate of the community is negative or positive, and that the coefficients B and D play a role in determining whether the community diversifies or not. This is consistent with our numerical results, which show that certain parameter regimes allow diversification and other regimes do not — even in the absence of general fitness differerences.

Distribution of biases and number of extinctions

The approximate model of the evolution of the biases makes predictions about the shape of the bias distribution as a function of attempted invasions. Before the onset of evolution, the distribution of biases is a truncated gaussian, with a lower cutoff set by ξc. However as the distribution evolves according to Equation A8.3 with the absorbing boundary condition at the critical bias, it smooths out near this cutoff, going linearly to 0 as ζξc+Υ¯ (or, equivalently, as ξξc).

Simulations confirm the expectation that the typical bias scales as 1/L over one order of magnitude in L (Appendix 3—figure 4). As predicted, one observes a smoothing out of the bias distribution toward ξc when comparing evolved and assembled communities of the same L (Figure 5B), and our analysis allows us to obtain the theory curve for the evolved ecosystem in Figure 5B as the solution of the approximate boundary value problem. However, the critical bias is sufficiently negative that the number of strains affected by the differences between the initial and evolved communities is small and the distinctions hard to see numerically.

Distributions of mean abundances and biases before and after evolution with ρ=0, si=0.

(A) Mean abundances: all these communities have L500, but evolved communities have diversified from K=50 initial strains so that they lose memory of their initial assembly conditions, while assembled communities had K=650 strains with L0500 surviving after the initial epoch. Data are pooled across 10 simulation runs of each. Though the distributions are mostly similar, there is a marked depletion in both rare and abundant strains in the evolved community. (B) The low end of the bias distribution changes from a truncated gaussian for the initial unevolved community (blue), to a linearly vanishing function (orange) because of evolution-driven extinctions of close-to-marginal strains. Bars show histograms from simulation, and solid lines show theory as detailed in ‘Appendix 8’. Inset shows the normalized bias by rank order, illustrating the smoothing of the lower end of the distribution caused by the evolution.

However, the density of biases near ξc determines the response of the community to evolutionary perturbations, since these low-bias strains are the ones most susceptible to extinction. In particular, the predicted linearly vanishing density of biases determines the distribution of the number, , of extinctions per successful invasion (Appendix 3—figure 1B). To estimate this distribution — particular the probability that is large — we use the fact that an invader will perturb the extant strains' biases by a random amount of order 1/L and proportional to the mean abundance of the invader. The positive tail of the invader’s mean abundance, ν¯A, is gaussian, since for positive ξinv, ν¯invξinv and the invader’s bias ξinv is itself gaussian distributed. The number of strains whose biases are within ν¯A of ξc is proportional to L2ν¯A2, because the distribution of strains' biases vanishes linearly at ξc. Thus for fixed ν¯A, the number of strains that are driven extinct is Poisson distributed with mean proportional to L2ν¯A2: this is of order one for large L as expected. That the tail of the distribution of ν¯A is gaussian implies ν¯A2 is approximately exponentially distributed in its tail. Integrating the Poisson distribution over this yields, for large , P()eβ/ (with β an order-unity coefficient) which is close to exponential as observed in Appendix 3—figure 1B.

A similar analysis for the initial randomly-assembled community shows that for fixed ν¯A, the mean number of extinctions triggered by the first successful invasion is of order L03/2ν¯AL0; much larger than after evolution has proceeded for a while. As the Poisson with this mean has a narrow distribution, the probability of an anomalously large number of extinctions will be dominated by the gaussian tail of the distribution of ν¯A and hence itself be roughly gaussian, though unless the initial L0 is huge, the tail is unlikely to still be in the asymptotic regime. The transient caused by a set of early invasions will likely cause a total of order L0 strains — with a relatively small coefficient — to go extinct before L starts steadily increasing, and this will occur over of order L0 invasions. For γ=-0.8 this effect appears to be very small — the critical bias is quite negative — but for smaller or larger γ the effects are noticeable (Figure 2).

The distribution of mean abundances, ν¯i, is related to that of the biases via the function N(ξ): therefore, we expect this also to evolve as the community diversifies. In particular, there should be a reduction in the number of strains at low mean abundance, since these correspond to those with close-to-marginal bias. In Figure 5A, we see that the mean abundances in an evolved community are more narrowly distributed than in an assembled community, with both fewer highly abundant and fewer rare strains. This is consistent with our picture of the bias distribution being smoothed out toward ξc due to invasion-triggered extinctions, resulting in the depletion of low-abundance strains. The depletion of abundant strains is likely due to the kill-the-winner dynamics which rewards invading strains that push the most abundant extant strains down.

Although the mean abundances are not broadly distributed on a log scale, the snapshot abundance distributions are, as seen in Figure 1. Note that most widely-used measures of diversity are not really informative for these kinds of logarithmically broad distributions. For example, the Shannon entropy would weight mostly the highly abundant strains, while the “species richness” would be highly sensitive to the lower cutoff in observable abundance.

Diversification rate with a distribution of general fitnesses

Armed with understanding of the scaling of the bias and mean drives with L, we can build upon the analysis of the simple exponential distribution of general fitness (‘Evolutionary dynamics with general fitness differences’) to analyze the evolution when P(s) decays faster or slower than exponentially. A heuristic understanding of how the dynamics of L depend on P(s) follows from the fact that without general fitness differences, the biases are distributed with characteristic scale 1/L. As L increases, the distribution of these biases gets narrower, and the system becomes progressively more “neutral" with overall differences in strain biases becoming smaller. The contribution of the general fitnesses is to add a random extra piece to each bias, broadening the distribution of extant ξi. In the limit of many invasions, the width of the drive distribution becomes comparable to the width of the distribution of the extant si and cannot decrease further. Thereafter, the shape of P(s) determines both the width of the bias distribution, and the number of coexisting strains.

If L1/Σ2 initially, the si play little role and the population-weighted mean fitness, s^, only increases gradually. But once s^ is a few times Σ, the tail of P(s) will determine the rate of increase of s^. Henceforth, s^ will grow steadily with subsequent successful invasions, since strains with sis^ are very unlikely to persist, and strains with sis^ are unlikely to have yet occurred. Thus, the range of extant si will become much narrower than Σ. This implies that the distribution over the currently relevant range can be approximated by an exponential distribution P(s)es/Σ^ with the effective width of the extant si distribution given by

(4) Σ^=1d[logP(s)]/ds|s=s^=Σψs^ψ1 ,

where the second equality is for the specific models we study (Equation 2). Provided evolution has proceeded long enough that no strains with si smaller than the original scale Σ survive, Σ^ will vary slowly for a range of evolutionary time and this sets the scale for variations of si of both extant strains and of potentially-successful invaders. This suggests that understanding the general behavior at long evolutionary times can be built on understanding the case of the simple exponential distribution (ψ=1 for which Σ^=Σ). The main difference is that now the community size will change as L1/Σ^2, with Σ^ changing as s^ increases: this will govern how L changes as invasions are attempted and occasionally occur.

In the slow evolution regime at long times, successful invaders must have sA comparable to s^. A simple argument gives an upper bound on how fast s^ can increase with T. In order to get a mean of s^ in a community of L strains, at least L attempted invaders must have occurred with sAs^. This requires a number of invasion attempts, T, such that T0s^(T)P(s)dsL(T). But if s^ were substantially smaller than this upper bound, many strains would have already occurred with si-s^Σ^ and these would persist for a long time, driving s^ up. We thus make the Ansatz, justified by the analysis and simulation data of ‘Appendix 6’ and Figure 6, that s^(T) grows with log(T) at a rate asymptotically given by this upper bound. For the distributions of interest we then have

(5) TP(s^)L(T)s^Σ[ψlog(T/L)]1/ψandΣ^Σ=(Σs^)ψ1LΣ2[log(TΣ2)]22/ψ,

where the last implication is due to LΣ^-2. These scalings become valid once the distribution of the extant si is pushed into the tail of P(s). To crudely take into account the effect of a large initial number of strains when K1/Σ2, log(T) can be replaced by log(K+T), as for the plots of L(T) and s^(T) in Figure 6. Figure 6B shows the theoretical prediction for s^(T) in the simplest case of ψ=1, and we see that s^(T)/Σ is reduced from the log[(T+K)Σ2] prediction by only an O(1) constant, as expected.

Evolutionary dynamics for unrelated invaders (ρ=0), with general fitnesses drawn from distributions parametrized by various values of ψ: faster-than-exponential decay for ψ>1 and slower-than-exponential for ψ<1.

Data are shown averaged over 50 replicates, conditional on not crashing, starting from K=300 initial strains, with the shaded region showing the standard error. For ψ=0.8, which results in decreasing diversity and crashes, a few individual trajectories are shown instead of an average. (A) Size of community as a function of the total number of strains introduced, T+K. For very long evolutionary times, we expect L[log(T+K)]2-2/ψ, but transients due to initial conditions are substantial. In order to push up into the tails of the s distributions, the parameters of P(s) are chosen differently for each ψ for ψ=0.8,1,2,3 respectively. (B) Increase of the community-average s^ with T, shown pushing into the tail of P(s). For ψ=1 the dotted line shows the theory prediction s^/Σlog[(T+K)Σ2], with deviations from this expected to be an O(1) constant for large T.

At long times, the probability of successful invasion decreases very rapidly with s^ and, as we show in ‘Appendix 6’, the cumulative number of successful invasions for ψ1 grows very slowly, with ZΣ-2(logT)3-2/ψ. Nevertheless, for ψ>1, the number of strains grows without bound albeit as a sub-linear power of the cumulative number of successful invasions. The average number of extinctions per successful invasion gradually decreases towards one as P(s) in the pertinent range, ss^, becomes closer and closer to exponential.

For longer-than-exponential tails, ψ<1, the diversity will decrease (possibly after an initial increase if K1/Σ2) and eventually — in practice rather soon — crash as seen in Figure 6A.


In this paper, we have answered an important issue of principle: Without any assumption of niche-like differences between strains, can diversity continually grow under slow evolution? We have found that this can indeed occur if the community forms a spatiotemporally chaotic phase that we have studied previously in PAF. As new strains are introduced — either separately evolved invaders or mutants of extant strains — some successfully invade, potentially driving extinctions of strains in the community. In a range of parameters, the size of the persistent community continually grows on average, while for other parameters, the diversity decreases and eventually crashes. How fast the diversification proceeds depends on the statistical properties of the strains. If each strain has a different general fitness, si, then as evolution proceeds the average si of the population gradually increases and pushes into the tail of the si distribution. If this tail falls off faster than exponentially, the community continues to diversify but more and more slowly, since fewer new strains will have sufficiently large si to invade. For broader-than-exponential distributions, the diversity eventually crashes as the general fitness differences dominate over the effects of interactions with other strains.

Building on an analytic and scaling understanding of the STC phase for an assembled community, we have developed a substantial understanding of the dynamics of the diversification or de-diversification. However even for the simple models on which we have focused, there are aspects that we do not understand.

Unresolved issues with the simple island models

Development of correlations

Even with invaders uncorrelated with the extant strains, subtle correlations build up in the interaction matrix and — although they appear rather weak (‘Appendix 4’) — the memory of earlier evolution will affect the way strain abundances change under further evolution, potentially mandating a better treatment of the evolution than the Markovian approximation we have used in ‘Evolution without general fitness differences’. With several complicating features — mutants, correlated general fitness differences, and substantial-sized initial communities — included, there are a number of crossovers that we have not attempted to analyze (‘Appendix 7’). These, and which aspects promote, slow down, or prevent, continual diversification, are likely to be quantitative and strongly model-dependent.

Nucleation of diversifying “phase’

An observation from the simulations (‘Continual assembly and diversification’) gives rise to a broader question: Why is it so hard to nucleate the diversifying STC phase? And, concomitantly, why do initially diverse communities so often crash unless the diversity is rather large? It is likely that the limited number of strains that dominate on each island over any short time interval — of order L/M — plays a role, but unclear how. Whether the difficulty of nucleating a diverse STC community is special to the structure of the models and spatial dynamics assumed, or is true more generally, certainly needs further investigation.

Spectrum of mutants and coexistence of parents and mutants

When invaders are mutants of extant strains that differ from their parent only very slightly, (with correlation coefficient ρ very close to unity), we have found that the parent and mutant coexist surprisingly frequently. Understanding this, even for the first mutant, requires analyzing the dynamics of strains with strongly correlated noise which we have not carried out, although we suspect that the very large local abundance variations that occur with low-migration rate give rise to a small decorrelation scale needed for coexistence. In each simulation, we have considered only mutants with a fixed level of correlation with their parents, leaving a number of natural questions: What are the effects of a distribution of magnitudes of mutational differences? How do these affect the invasion, coexistence, and subsequent properties of the evolving communities?

Invasion dynamics

Because of the local chaos and low migration, the invasion of a potentially-successful new strain is complex. To avoid this complication, we have introduced new strains at substantial abundance and on all islands simultaneously. In actuality, most initial invasion attempts on an island will fail: only if the strain arrives when the conditions are ripe for it to bloom, can it avoid quick extinction and send out enough offspring to other islands, which — if also sufficiently good timing — allow it to spread. How this process depends on the relatedness of mutant and parent complicates matters greatly because of the boom-bust dynamics. Strains are most likely to beget mutant offspring when their abundances are high, but at that stage of a bloom, a crash in the local population will soon follow. Therefore, although many mutants may arise when a parent strain is doing well, the correlation between their dynamics and those of their parent means that they are likely to quickly go extinct when their parent crashes down from high abundance. In contrast, mutants that emerge right before a parent blooms up to high abundance can ride the bloom and establish more readily, but would have to arise in a small parental population. Understanding the balance between these effects and their consequence for invasion probabilities is a challenge for future work — especially with real spatial structure and dynamics, discussed below.

Spatial structure and dynamics

While for some microbial populations — for example common human gut commensals — a collection of connected “islands” without much spatial structure may be a rough caricature, for most populations there is spatial structure that makes dispersal from one location to another dependent on the distance in one, two or three dimensions. Thus, instead of having all pairs of islands connected by migration, one could model a d-dimensional array of islands with nearest-neighbor migration; a spatial continuum with diffusive dispersal; or a mixture of long and short distance dispersal events as driven by wind, ocean currents, or hitchhiking on migrant animals (Hallatschek and Fisher, 2014). With real spatial structure, local sub-populations are much more prone to extinction and cannot be as readily rescued by migration from another location where the strain is blooming. Thus, in contrast to the regime we have worked in for this paper, recovery from local extinctions must play a crucial role. The dynamics of invasions, extinctions, and repopulation is very different than in the spatial mean field model: if the underlying dynamics is diffusive, invasion and repopulation will occur by propagating Fisher-Kolmogorov-Petrovsky-Piscunov (FKPP) fronts (Fisher, 1937). The properties of FKPP waves are known to be highly sensitive to dynamics at the wavefront, and the effects of demographic fluctuations have been investigated (Korolev et al., 2010). But the approximately multiplicative “noise” from the ecological interactions will surely change this, and even for a single wave understanding the impact of these larger fluctuations is still an open question (Rocco et al., 2002; Rocco et al., 2000).

With long-range dispersal over a multitude of length scales, the dynamics of invasion, extinction and repopulation will be very different, as already occurs for a single successful invader without ecological variations (Hallatschek and Fisher, 2014). Generally, understanding of the STC phase will have to build on better understanding of repopulation dynamics in the presence of large ecological fluctuations, and then understanding the evolution of communities on top of that. We leave investigations of this for future work. But we conjecture that a continually diversifying STC phase can still occur with more realistic spatial dynamics.

Bacteria-phage interactions and coevolution

An obvious weakness of the Lotka-Volterra models studied here is that the strains do not carry their own phenotypes, but are characterized by their interaction with all possible other strains. Furthermore, the antisymmetric correlations in the interaction matrix (especially without substantial general fitness differences) are rather unnatural for multiple strains of a single species. Thus, the most interesting extension of this work is to much more natural models: multiple strains of a phage species that prey on multiple strains of a bacterial species, with varying effectiveness that is a function of phenotypic properties of the particular phage and bacterial strains. Of particular importance is the interaction between a phage tail and bacterial receptor, as modelled in Weitz et al., 2005. We showed previously [PAF] that the block-antisymmetrically-correlated structure of the interaction matrix with the bacteria having no niche-structure (differing only in the way they interact with the phages) can give rise to an STC phase that is very similar to that of the antisymmetrically correlated Lotka-Volterra model studied here: a similar model was further explored in Martis, 2022. Such a bacteria-phage model can naturally accommodate general fitness advantages through phenotypic changes, eliminating the need to introduce them on separate footing. In ongoing work, we show that much of the basic phenomenology we have found here also occurs in evolving bacteria-phage phenotype models — at this stage only roughly and qualitatively.

For bacteria phage models, studying phylogenies and relatedness questions are natural. Whether more specialist phages tend to evolve, making the interaction matrix sparser and perhaps more hierarchical — and if so under what circumstances — is a particularly interesting question.

Concluding questions

We have studied evolution of communities of many closely related strains in the limit that the evolutionary dynamics is slow compared to ecological and spatial dynamics. For a class of models, and in a particular ecological “phase”, evolution drives continual diversification, provided there is sufficient diversity initially. However mutations that change general fitness of strains tend to strongly slow down or even reverse the diversification. Thus we ask: How ubiquitous is diversification in the absence of any niche-like structure? Are there models in which a diversifying phase is easier to nucleate? Will the diversification always tend to be limited or strongly-slowed by general-fitness mutational effects? Or might “entropic” effects associated with difficulty of finding such general fitness mutations — for example from discrete genomes rather than continuous phenotypic parameters, or from soft tradeoffs — counter this slowdown, or perhaps produce evolutionary dynamics that lead to sparse interaction matrices and broader distributions of biases? Conversely, if strains are initially separated in “niche space” but then start to overlap and interact as the number of strains increases, how does the behavior differ? Is continual diversification easier to nucleate? Are the statistical properties of the phylogenies resulting from this evolutionary process — here driven entirely by “selection” in the broad sense, with ecological interactions creating a balance between the many extant strains — similar to a known class of coalescent trees?

What happens when, as in large microbial populations, evolutionary processes are not slow? Faster evolution is likely to make diversification easier, but understanding this even in simple models will require much better understanding of the invasion probabilities of mutants. Other than our scenario in which spatiotemporal chaos is the key to stabilizing coexisting diversity, what other robust continually diversifying scenarios are there? And of course, most crucially, what observable features of the strain, sub-strain, and sub-sub-strain level diversity in a microbial population (or interacting populations) could provide hints to the underlying causes of extensive diversity?

Appendix 1

Spectrum of timescales

The primary timescale that governs the interplay of evolution and ecology is the ecological relaxation time which is O(ML), with Mlog(1/m), as discussed in ‘Timescales’. However the STC exhibits a number of other timescales. Although it does not play a role in the current work, there is a short timescale associated with the dynamic fluctuations: roughly the time that a strain spends near the peak of a bloom, which is similar to the inverse of its instantaneous growth or decay rate. This is dominated by the net effect of its interactions with the small subset of O(L/M) strains that happen to be abundant at that time: the variance of this is jνj2M/L, which makes the dynamic fluctuation time L/M. Between this dynamic fluctuation timescale and the timescale for blooms, correlations in the growth rates decay as a power of the time difference. During a bloom, each strain experiences multiple reversals from growth to decay: this is a special property of the self-organized chaotic state.

The time to go extinct for a strain destined to do so depends logarithmically on the extinction threshold 1/N, but as long as Nm is very large, whether extinctions occur is not strongly dependent on N, as analyzed in PAF. In our simulations we choose for convenience Nm1 so that logN is a few times M.

The timescale for migration to be effective would, if there were no differences between the strains, be of order 1/m which is very long. However, the spatial dynamics are much faster than this because of the exponential growth of local populations when they happen to be in a favorable community. This makes the timescale for exponential spread across islands of a successful invader be on the order of the bloom time on a single island, ML. This is analogous to the rapid spread of a Fisher wave driven by selection, even when the spatial dynamics is diffusive Fisher, 1937; Korolev et al., 2010.

The timescale of demographic fluctuations — even if some strains were phenotypically identical — would be very slow Nτgen with τgen1 (in our units) a generation time. In practice, these fluctuations only matter when the populations happens to be very small and is being driven extinct on the deterministic dynamic timescale L/M. But as long as Nν¯mτgen1, there are many migrants arriving and the dynamics is essentially deterministic. If the island-average ν¯ drops enough then the migrations become stochastic with time intervals between them of order 1/(Nν¯m). But when this occurs for substantial time, the global population is likely to be on the way to extinction. We focus on NmL for which the population dynamics of the persistent strains are not strongly influenced by the stochastic or migratory demographic fluctuations.

Appendix 2


Parameters and integration

For all numerics, parameter values unless otherwise mentioned are I=40, m=10-5, γ=-0.8, N=1020. In order to integrate the dynamics, we use an adaptive forward Euler step, with the time step chosen so that the maximum fractional change in the abundance of any strain is no greater than 3/4. This means that the time step scales with the dynamic fluctuation time (‘Appendix 1’) which is O(L/M). Therefore, in a single time step an abundance can change to anywhere from 1/4 to 7/4 of its current value. Because of the wide range over which abundances vary, such large changes do not cause numerical problems. Extinctions and invasions are treated deterministically, with local extinctions occurring for strain i on island α if νi,α<1/N. In this case νi,α is set to 0. An irreversible global extinction of strain i occurs if νi,α=0 for all α simultaneously. Recolonization of a locally extinct strain happens when Nmν¯dt>1, where dt is the time step of the integration. Since this time step is chosen to reflect the basic timescale of the abundance dynamics, this choice of the recolonization threshold is consistent.

At the start of a new epoch, one strain is introduced into the community: this is done deterministically, with fractional abundance of the new strain set at 1/L on all islands and the other strains” abundances adjusted proportionately to maintain the overall population constraint. To save computational time, some incoming strains are rejected because they are very unlikely to successfully invade. This can be estimated by calculating their biases from the properties of the community before their introduction. By roughly estimating the critical bias from earlier successful invasions, we can conservatively reject some incoming strains without having to run the dynamics. With general fitness differences this can yield substantial speed-up by, for a community with population mean s^, rejecting strains with s-s^ sufficiently negative.

Since the epoch length is O(ML) (set by the ecological relaxation time) and the time step scales as L/M, with each timestep requiring O(L2) computations to compute the instantaneous growth rates of the strains, the runtime of a single epoch scales as O(M3/2L5/2).

Correlations between mutant and parent interactions

Mutation of a parent strain to create an invader is comprised of two parts: changes in the parent’s interactions with the other strains, and a change in the parent’s general fitness. To generate the interaction part of a mutation, we append a new row and column to V which parameterizes the interaction between the new strain M and the parent strain P, each with another strain, k, according to E[VMkVPk]=E[VkPVkM]=ρ for M,Pk. In order to preserve the correlation γ between across-diagonal entries of V, we take VMk=ρVPk+1-ρ2Z1 and VkM=ρVkP+1-ρ2(γZ1+1-γ2Z2) where the Zi are i.i.d. standard normal random variables. This preserves the desired correlations, with E[VMkVkM]=γ.

However, we have to treat the direct interactions between the parent and mutant, k=M,P, more carefully, since it is not always possible to preserve γ while also having the desired correlations between VPP, VPM and VMP. Defining D, H, J and K as


the symmetry conditions we enforce are

(A2.2) E[D2]=E[K2]=1+γ,E[H2]=E[J2],E[DJ]=E[DH]=E[KJ]=E[KH]=ρ(1+γ).

In addition, we require that D=H=J=K in the limit ρ1, and E[HJ]=γ as ρ0. The parameterization that we choose which respects these conditions is

(A2.3) H=ρD+1ρ2Z1
(A2.4) J=ρD+1ρ2(γZ1+1γ2Z2)
(A2.5) K=ρ2D+ρ(H+J)+(1ρ2)1+γZ3.

where the Zi are i.i.d. standard normal random variables. Now we have E[HJ]=ρ2+γ(1-ρ2), which does not precisely preserve the desired correlation E[HJ]=γ, but it has the correct limit as ρ0. Our parameterization gives E[DK]=ρ2(1+γ), but it is possible to choose an alternative parameterization which has a different value of E[DK] but still respects the desired symmetries.

The dominant interaction between the parent and the mutant is mediated by the other strains in the ecosystem, with the direct influence of the parent on the growth rate of the mutant and vice versa smaller than the contributions mediated by all the other strains by a factor of L. Therefore the choice of the direct parent-mutant interaction ought not to be of much importance — unless sufficiently strong correlations develop under continuing evolution, for which we do not see evidence.

Epoch length in the diversifying phase

At the beginning of each simulation run, an ensemble of K randomly drawn strains is assembled and run to reach a chaotic steady state under the dynamics of Equation 1 for time O(MK). Between the introductions of new strains the ecological dynamics are run for epochs of length cepochML with cepoch=3. In addition, for most simulations we intersperse regular-length epochs with longer epochs to get rid of any marginal strains which may be barely surviving. These longer epochs (longer than the others by a factor of 3) occur every 100 epochs.

The duration of the initial epoch with an assembled community of K strains was usually chosen to be shorter than later epochs. But for the evolutions this does not much matter as the distribution of close-to-marginal strains is different than in later epochs and these are in any case quite likely to go extinct in subsequent epochs. However, for the simulations for which we compared assembled and evolved communities (e.g. Figure 5), the dynamics of the initial assembled communities were run for a longer time 12ML.

In order to check that the evolutionary dynamics are in the quasistatic limit, where the strains are introduced slowly enough that their introduction rate does not matter, we look at how the diversification rate depends on the length of the epochs between strain introductions. In Appendix 2—figure 1 we vary cepoch over one order of magnitude and find similar evolutionary dynamics for all epoch lengths. In the case of ρ=0, there is a longer transient for the long epochs, though once in the diversifying regime, L increases at a similar rate for all epoch lengths. Based on this, all data presented in the paper were generated with cepoch=3, which makes numerics faster, and should not change qualitative conclusions about the slowly diversifying regime.

Appendix 2—figure 1

(A) ρ=0. (B) ρ=0.95. Rates of diversification (measured per attempted invasion) are insensitive to a factor-of-10 change in the intervals between the introduction of the new types: the epoch durations shown are cepochML with several values of cepoch. Although it can take more invasions to get into the steadily diversifying regime for long epochs, the rates of diversity increase depend little on the epoch lengths. Both dashed lines both have slope 0.2 (diversification rates happen to be similar for ρ=0 and ρ=0.95).

Estimation of drive and bias

The mean drive and bias are emergent quantities from the mean field analysis (‘Appendix 5’), and are therefore not manifest in our direct numerical integration of the dynamics: however we can extract them from numerics as detailed below. Here, we use the notation ν¯ij for the average abundance of strain i in the absence of strain j. The averaged growth rate of a strain, excluding migration, can be measured directly. In the mean field approximation, this is decomposed as the sum of its bias and feedback from it perturbing the other strains::

(A2.6) r¯i=si+jVijν¯jΥ¯
(A2.7) =si+ζiΥ¯+γXν¯i
(A2.8) =ξi+γXν¯i,

where the drive, ζi=jVijν¯ji, can be written in terms of mean field quantities as ζi=jVijν¯j-γXν¯i. The term γXν¯i captures the feedback of strain i back onto itself via the other strains. Therefore, if we know the susceptibility, X, we can calculate the drive and bias of a strain: since r¯i and ν¯i are measurable in our simulation, we need only to subtract off the γXν¯i term from r¯i to get ξi. We can find X using the self consistency condition that X=idν¯idξi=idν¯idζi.

First, we make a guess for the susceptibility X which allows a provisional calculation of the ζi. Then, using our numerical simulation data, we fit parameters a, b and c to a functional form giving the mean abundance in terms of the drive and s:

(A2.9) N(ζi)=blog[1+exp(ζi+sica)].

The justification for this form is that it gives ν¯iζi+si-c for ζi+sic>0. This is the expectation that ν¯i becomes proportional to ξi for positive ξi — though here we include an offset via the parameter c which roughly represents the expected correction to the linear function N(ξ) at large bias. The quantity ζi+si-c is similar to ξi but not quite the same due to a systematic difference between c and Υ¯. For large negative bias, the fitted form captures the roughly exponential decrease of ν¯ due to the rareness of blooms.

The data are fit by this functional form quite well (Appendix 2—figure 2). We observe that c tends to be smaller than the measured Υ¯ within a given epoch, which could be due to the effect of migration which elevates ν¯i for strains with negative ξi. The parameter combination b/a is roughly equal to -1/γX, as expected from the mean field analysis. From the hypothesized form of N(ζi) we can recalculate the susceptibility via X=iΔν¯iΔζi (where the derivative is calculated numerically). We then update our value of X by averaging the old guess with the new estimate.

Appendix 2—figure 2
Inference of the drives, ζi, for (A) communities evolving with ρ=0 and (B) communities evolving with ρ=0.95.

In both cases, all si=0. The black lines shows the fit of the functional form for N(ζ) after X has converged to a self-consistent value. Data are pooled (both for fitting and for plotting) over 5 consecutive epochs, each with 300 extant strains. For ρ=0.95 it seems that correlations in the interactions that have accumulated due to the relatedness are affecting the relationship between ν¯ and ζ. Certain strains are visible outliers from the average dependence of ν¯ on ζ, and appear multiple times on plots since data are pooled across 5 epochs.

By iterating over the susceptibility until it converges to a value where the assumed drives ζi=jVijν¯j-γXν¯i reproduce the susceptibility via the self-consistency condition, we can get an estimate for the susceptibility and therefore the drive as well. Then the bias is ξiζi+si-Υ¯, where we have subtracted the Lagrange multiplier (neglecting effects of order Υ¯-Υ¯i which are smaller by of order 1/L) and added back in the general fitness. We can validate that our bias estimator is reasonably accurate by comparing the inferred biases to the biases of strains in their first epoch of invasion — the latter can be measured directly from simulations since they are simply the invasion eigenvalues for incoming strains. As this works well, we use it to infer the biases for all the other strains in the community. (Note that in contrast to the assembled communities studied in PAF, we cannot use the condition that the drive averaged over all the initial K strains, is zero. PAF used this — along with the expected truncated gaussian shape with variance iν¯i2 and lower limit ξc+Υ¯ — as a check in calculating the drives. This no longer works here because conditioning on evolution means that the mean drive is no longer and the drive distribution is no longer a truncated gaussian.)

To increase the quantity of data on which to perform our fit for N(ζ) in a given epoch, we pool data from up to 20 epochs around the focal epoch, provided their L is within 5% of the L in the focal epoch, with the hypothesis that these communities are statistically similar and therefore have a similar N(ζ).

The function N(ξ) for the evolved communities with ρ=0 looks quite similar to the function for assembled communities However for communities evolved with ρ=0.95, N(ξ) looks substantially different. This indicates that the build-up of correlations in the community requires modifications of the independence assumptions of the mean field theory.

After inferring the function N(ζ) and the drive for each strain, we can calculate the fragility, defined as Ξ=jχj2/(1-jχj2). The closer jχj2 is to 1, the more unstable the community is to perturbations. Indeed we find that for simulations diversifying from 50 to 500 strains, jχj20.6 for the evolved community of 500 strains, while early in the simulations, when L<100, we find jχj20.75. When comparing jχj2 for assembled and evolved communities, both of 500 strains, there is not a clear difference. Note that the close-to-marginal strains, which are more abundant in the assembled community, only contribute small amounts to Ξ, so this is not surprising.

Appendix 3

Dynamics of diversification

Extinctions, diversity crashes, and nucleation from low diversity

Appendix 3—figure 1
Crashes and distribution of number of extinctions per invasion.

(A) The fate of an evolving community depends on the number of initial strains K. The characteristic size below which crashes dominate is 50 for the parameters γ=-0.8 and ρ=0 as shown here. (B) The distribution of the number of extinctions per successful invasion depends on the size of the community. Data are aggregated across a range of K, each run with 100 replicates. For small community size, L, invasions occur in which a substantial fraction of the strains go extinct. But for large L, multiple extinction events are very rare and the distribution is close to geometric: the dotted line is P()=(1α)α for the probability of extinctions, with α0.41, corresponding to an average of 0.7 extinctions per successful invasion as shown in Figure 2A.

A crucial question is whether it is possible to build up a highly diverse community from a small initial number of strains and, if so, on what this depends. Starting simulations at values of K between 10 and 90, we observe a crossover size, K*, of the initial number of strains above which diversification is robust, and below which the diversity typically crashes under the evolutionary dynamics and does not again increase substantially for the duration of the simulations, (Appendix 3—figure 1A). The presence of a crossover K* suggests that once the system is sufficiently diverse, it tends toward further diversity. Thus the main obstacle to diversification in these models is going from a single strain to 50 strains. We define K* heuristically as the lowest K that diversifies with probability more than 1/2 (see Appendix 3—figure 2). For each K, there is a corresponding L0(K): the number of strains that persists after the initial drop in diversity. This somewhat-variable L0 appears to be the primary determinant of a community’s stability to evolutionary perturbations. However this stability also depends on the community’s history. For example, if the community has evolved gradually from a smaller number of strains to some current size LE, its response to continued evolution may be different than if it has been assembled from an initial number K that dropped through its initial evolution to an ecologically stable community of size LE. This motivates definition of a crossover size L*, different from L0(K*), given by the minimum size of an evolved community for which the probability of entering the diversifying phase is greater than 1/2. We discuss the differences between assembled and gradually evolved communities further in ‘Appendix 4’, but leave an analysis of how these properties affect the nucleation probability for future work. The size L* is a natural quantity to work with, as it is the relevant one when considering how a community might nucleate from a small number of strains into a diversifying regime. As we discuss later, we expect L<L0(K) since just-assembled communities are more susceptible to perturbations than gradually assembled ones.

Appendix 3—figure 2
The probability of establishing the steadily diversifying regime increases with the initial number of strains and decreases with the migration rate.

Although changing the migration rate from 10–5 to 10–10 changes the establishment probability significantly, further reducing m has little noticeable effect, indicating that the dependence of establishment probability on m is rather weak over several orders of magnitude. Data are averaged over 100 simulations for each set of parameter values. If m=0, then migration cannot stabilize the spatiotemporal chaos, and so the establishment probability should vanish — but this is a very singular limit.

We can quantify how likely the system is to diversify or crash by inspecting the number of extinctions, , that occur for each invader that successfully enters the ecosystem. The distribution of is dependent on L and shows that more diverse ecosystems suffer on average fewer extinctions triggered by invasion of each new strain (Appendix 3—figure 1B). For small communities L of order L* or less, extinctions of a substantial fraction of the strains occur, and cascades of extinctions in response to successive invasions cause L to crash. The fragility of the communities with LL* to evolutionary perturbations and their strong tendency not to recover after a crash (Figure 2A), implies that the process by which a low diversity community could diversify is very different than the steady diversification of already large communities. The transition from the low diversity to the diversifying regimes must be mediated by a very rare nucleation event in which L becomes roughly larger than L*. We discuss the low diversity regime and speculate about such nucleations in ‘Appendix 3’.

In the steadily diversifying regime the distribution of extinctions per successful invasion is close to exponential — a result that we derive in ‘Distribution of biases and number of extinctions’. For large L, the chances that a substantial fraction of the strains go extinct is extremely small and decreases exponentially as L increases further: thus for large L the continual diversification is essentially deterministic.

Probability of diversification as a function of initial community size

The initial community size beyond which diversification becomes likely depends m, γ and ρ, and is defined as K*. For some range of γ (Figure 2B), K* is infinite: the system always eventually crashes. The role of m in determining K* can be partially understood from the mechanism of diversity crashes due to synchronization of dynamics between islands. This occurs more frequently when m is larger, but its dynamics are subtle because of local blooms up to high abundance which can can have outsize effects. Though we have not studied the synchronization and diversity crash process here, it suggests an interpretation of our simulation results (Appendix 3—figure 2) which show that decreasing m lowers K*, presumably because it makes the STC more difficult to spatially synchronize and hence less likely to crash.

Nucleation from low diversity

Here we further investigate the evolutionary dynamics when the diversity is low. We are particularly interested in the nucleation process by which an ecosystem could transition from the LL* regime to the steadily diversifying regime, where L* is the crossover beyond which continual diversification becomes likely. For the parameters we have investigated, this turns out to be extremely rare. Even in runs of 105 attempted invasions starting from a single strain, the community always remains small. We observe that the number of extinctions per successful invasion is broadly distributed (as when starting with a larger community LL* in Appendix 3—figure 1) and the evolutionary dynamics of L typically proceeds by incremental increases punctuated by large decreases. In this regime, the steady state distribution of L is observed to be very close to Poisson (Appendix 3—figure 3) — a result that we do not have a theory for. From this Poisson behavior, one can attempt to extrapolate the probability of reaching L* and transitioning to steady diversification by a lucky fluctuation in L. Taking L*L0(K*70)55 based on Figures Appendix 3—figure 2 and Appendix 3—figure 1, this estimate suggests 1060 attempted invasions would be needed. However, the probability distribution of L will surely deviate from Poisson long before L*, so our prediction for the chance of nucleation is likely too small by many orders of magnitude, depending on the value of L for which the probability distribution starts to be significantly higher than the Poisson extrapolation. At this point we do not have even conjectures for the shape of this distribution, nor for how the value of L* changes with γ, ρ or m. At the marginal point γ=-1, we have the special antisymmetric model (‘Appendix 6’), and diversification becomes impossible, so near there we expect that L*, along with the fragility, diverges.

Appendix 3—figure 3
Distribution of community size L starting from a single strain, over 105 attempted invasions with m=10-5 and γ=-0.8.

The solid line shows a Poisson fit with mean 1.6, excluding L=0. This fit is remarkably good. The inset shows the trajectory of L as it reached its maximum value which occurred only once, followed by a crash back down.

Results from ‘Appendix 3’ suggest that the nucleation process may be related to desynchronization which is necessary to enter the steadily diversifying regime. In fact we observed that decreasing the migration rate m increases the mean of the Poisson-like distribution of L in the low-diversity regime, and also reduces L* (Appendix 3—figure 2), thereby perhaps substantially increasing the probability of nucleation — although it is still too improbable to observe.

Scaling relationship between bias distribution and L

Appendix 3—figure 4
Scaling relationship between the effective community size, L1/iν¯i2, and the inverse variance of the bias distribution, 1/σξ2, for unrelated invaders with s=0.

The dashed line has slope 0.47. Data are pooled over 3 runs from L=50 to 500 strains (inset) so that when L becomes large the initial conditions have been forgotten. Inset shows that L and L are proportional as expected. (The bend in L versus 1/σξ2 toward the end of the curve is likely due to stopping the simulation the first time L reaches 500, which truncates the curve asymmetrically.).

The central scaling relationship from which our theoretical results follow is that the width of the bias distribution scales as 1/L, which scales the same way as 1/Liν¯i2. In Appendix 3—figure 4, we confirm these relationships by calculating the standard deviation of the bias distribution in a diversifying community where L ranges over one order of magnitude. We find the the relationship holds for both L and L, which are proportional to each other with L0.74×L.

Diversifying phase with ρ close to 1

One of the important results of this work is that there is a continually diversifying phase of eco-evolutionary dynamics both for unrelated invaders and for small-effect mutations, across the whole range of ρ[0,1). In Appendix 3—figure 5, we illustrate the robustness of this phenomenon for highly correlated parents and mutants with ρ=0.97, running a large number of replicate simulations and seeing that about 10% of these enter the diversifying phase, when K=50. Although the community is less likely to nucleate into the diversifying phase than when ρ=0, L still increases linearly with attempted invasions once the community is nucleated.

Appendix 3—figure 5
The diversifying phase persists for ρ close to 1.

Here we show 99 simulation replicates, starting from K=50. Out of these, 11 nucleated into the diversifying phase, albeit with diversity increasing at a slower rate than with ρ smaller (Figure 2). Simulations which crashed are shown in red while those that entered the diversifying phase are shown in blue.

Appendix 4

Statistics of V conditioned on evolution

How do the statistics of the interaction matrix Vij change when conditioned on evolutionary history? How does this compare to conditioning on non-extinction in an assembled community? We ran a set of 10 replicates without general fitness differences over a period where the community grows from 50 initial strains to 500 strains, for both ρ=0 and ρ=0.95. From these, we obtained ensembles of interaction matrices conditioned on gradual evolution with ρ both 0 and close to 1.

We also examined the interaction matrices of communities that evolved with a steady state L, by using an exponential distribution of s with Σ=0.07, starting from K=50 and reaching steady state with L150 while evolving over the course of 4000 successful invasions. In this case one might be more likely to see correlations building up, since new strains might not dilute the correlations building up in V. For comparison, we studied assembled (but not evolved) communities starting with sufficiently large K to have L0 similar to the evolved condition.

Appendix 4—table 1
Definition of some statistics of the interaction matrix of a community. The standard deviation is calculated as Stdevi[{xi}]=1Li=1L(xi-1Lj=1Lxj)2 (the population standard deviation) and Corr[X,Y] is the Pearson correlation coefficient calculated for a sample (the covariance normalized by the product of standard deviations).

In Appendix 4—table 1, we list the summary statistics of the interaction matrices that we analyzed: γ^ is the empirical cross-diagonal correlation in the V matrix of extant strains, σV^ is the empirical width of the interaction strength distribution after evolution, μ^V is the average of the interactions and Δ^ is the standard deviation, over i, of 1LjVij, which acts like an effective si, crudely approximating the width of the bias distribution. Aside from γ^, all the quantities are defined including the diagonal terms in V.

Appendix 4—table 2
Means (up to corrections due to diagonal entries) and standard deviations of the chosen statistics from 500 realizations of the matrix V from the original ensemble with γ=-0.8, and two values of L pertaining to the different evolutionary conditions. Here we do not condition on either non-extinction or evolution.
mean (L = 500)−0.8100.0446
standard deviation (L = 500)0.0010.0020.0010.0015
mean (L = 150)−0.8100.0816
standard deviation (L = 150)0.00360.0060.00310.0046
Appendix 4—table 3
Statistics of the interaction matrix for both evolved and assembled communities. Evolved communities grew to 500 strains from an initial K=50 and assembled communities had L0500 after a single epoch of duration 12ML, both with si=0. The last row is from a simulation with exponential P(s) and Σ=0.07, which results in L150 at steady state. The mean of each statistic is reported across 10 simulation replicates. In Appendix 4—table 4 we show these statistics appropriately normalized with respect to the original V ensemble. Bold entries correspond to those more than 3 standard deviations away from the original ensemble mean.
assembled (Σ=0)−0.80141.00120.00260.0337
diversifying (Σ=0, ρ=0)−0.80151.00280.00610.0289
diversifying (Σ=0, ρ=0.95)−0.80821.04270.14220.0328
steady state (Σ=0.07, ρ=0)−0.80741.02050.01650.0977

To check for statistical significance, we normalize each statistic of the evolved interaction matrices by the standard deviation of the estimator of the corresponding statistic from matrices drawn from the original ensemble without any conditioning. This indicates how atypical the measured statistics are for the matrices of size L×L. We thus define quantities γ~, σ~V, μ~V and Δ~ to be the empirical matrix statistics measured in units of standard deviations away from their mean in the original V ensemble.

Appendix 4—table 4
Properties of the interaction matrix for evolved and assembled communities, displayed in terms of number of standard deviations from the mean in the original V ensemble. These are the same data as in Appendix 4—table 3, but normalized according to appropriate scale of deviations calculated from original V ensemble. Numbers in boldface have magnitude greater than 3, and are thus clearly statistically significant.
assembled (Σ=0)−−7.2
diversifying (Σ=0, ρ=0)−−10.4
diversifying (Σ=0, ρ=0.95)−8.221.4142.2−7.9
steady state (Σ=0.07, ρ=0)−

Although a considerable number of statistics are significantly different from their original ensemble values, few incur substantial changes, with the largest change by far occurring in the μV for the simulations with ρ=0.95; however even this is a relatively small effect.

Our results show that both evolution and assembly tend to make γ^ slightly more negative than the γ of the original ensemble, and that these processes also favor — as might have been expected — an increasing mean interaction strength, μ^V. Of particular interest, given previous discussion, is Δ^, which is similar to the width of the bias distribution. For L500, the statistics of the matrices from the original ensemble have Δ^1/5000.0447. However, once we condition on assembly or evolution, the value of Δ^ can be expected to decrease, because we eliminate those strains with negative bias, meaning that strains i with an anomalously small value of jVij will go extinct, making the range of jVij narrower. This is indeed observed in Appendix 4—table 3, at least for simulations without general fitness differences. Interestingly we see that Δ^ decreases with evolution more for ρ=0 than for ρ=0.95, which could be because, in the latter case, the elements in a single row of V are correlated with those of a high-abundance parent, which could increase the width of the jVij distribution. Nonetheless, we still expect that Δ^ scales as 1/L for evolved communities.

For the evolution with L roughly constant from exponential P(s), we see that Δ^ after evolution is actually larger than in an unconditioned matrix. However, here the bias has an extra contribution from the general fitness, so Δ^ is no longer as good an estimate of the bias distribution width (though it should still be similar).

We can compare these results to previous results of Bunin, 2016 in which he calculated the statistics of the interaction matrix conditioned on assembly in the stable fixed-point phase where niche interactions — large negative diagonal terms in V — stabilize the diversity instead of spatial structure. Bunin observed that, conditional on assembly, correlations between the effect of strains k and l both on strain i become negative, that is Corr[Vik,Vil]<0 with small magnitude of order 1/K. This is consistent with our observation that the distribution of the sum of these interactions for each strain, jVij, will have a smaller width than in the original ensemble. In addition, Bunin finds an O(1/K) shift in the mean of the interaction matrix toward less competition, which is consistent with μ^V>0. We also see slight but systematic changes in γ^ and σV^, particularly in the case of evolution of correlated mutants with ρ=0.95.

Although many of the changes of the interaction statistics of assembled or evolved communities are statistically significant, they are mostly very small, with the possible exception of μ^V for ρ=0.95. We conjecture that for ρ=0 all the effects are small by some power of L. But for ρ near one, this is less clear: whether for L large compared to some inverse power of 1-ρ, the correlations induced by evolution will still be small, or whether the correlations will persist for arbitrarily large L, we leave as an open question.

Appendix 5

Dynamical mean field theory

In this section we provide some more background on the dynamical mean field theory used as a basis for the heuristics and scaling arguments throughout the paper. The main idea of dynamical mean field theory is to reduce an interacting many-body problem into a single-body stochastic problem, with the statistics of the stochasticity to be determined self consistently. As discussed in PAF, a dynamical mean field theory analysis of Equation 1 results in autonomous stochastic integro-differential equations for each strain independently. Since the strains are statistically equivalent on each island, we drop the island subscript:

(A5.1) dνidt=νi(si+ζi(t)+γtR(t,t)νi(t)dtΥ(t))+m(ν¯i(t)νi).

Here, the correlation function C(t,t) of ζi(t), the response function R(t,t), and the island-averaged abundance ν¯i, must be determined self-consistently. In order to use this description for the evolved (in addition to assembled) communities we have assumed that the statistics of the V do not change substantially in evolved communities: we show in ‘Appendix 4’ that our simulation results are consistent with this Ansatz. Here, we have made the time-dependence in ζi(t) explicit — but elsewhere in the paper when we refer to the time-averaged or mean drive ζ¯i, we have dropped the overbar for readability.

In ecological steady state, some of the strains will have gone extinct, and the correlations and responses of the persistent strains will only depend on time differences t-t, and quantities will have well-behaved time averages equal to island averages, denoted by overbars. The average growth rate of strain i on a single island excluding migration is

(A5.2) r¯i=si+jVijν¯jforcefromallstrainsΥ¯=si+ζi+γXν¯idriveandfeedbackΥ¯=ξi+γXν¯i,

where the static susceptibility, X, is the time-integrated total response function X=tR(t,t)dt=iχi. Here, we define the individual strain static susceptibility as χidν¯idξi. In steady state with migration, dlogνidt=0; therefore r¯i<0 must be exactly compensated by the average effect of the migration term m(ν¯i/νi-1), which is always positive (by the Cauchy-Schwarz inequality) and can be much larger than m when νi is small.

The stochastic DMFT equations cannot be reduced to equations for the correlation and response function and must be treated directly. Their self-consistent solutions in the STC phase were analyzed in our previous work [PAF], by asymptotics in the large parameter M. This analysis yields super-diffusive random walks of the log-abundances, persistence around the migration floor, and the statistics of the occasional blooms that occur for all persistent strains with negative bias. These aspects together determine how the mean abundance of a strain depends on its bias, via ν¯N(ξ) with corrections smaller by factors of O(1/L).

When the bias of a strain is positive the average input from migration will be relatively small, and hence r¯i0 which implies that ν¯iξi/(-γX). For strains with negative bias, the migration is essential and the statistics of the blooms, which are rarer the more negative the bias, make ν¯i exponentially small in -ξi [PAF]. For strains with more negative bias, ξi<ξc<0, the migration cannot sustain their local populations and they go globally extinct with their spatially averaged frequency decaying exponentially in time.

The distribution of the biases of a community plays a crucial role in determining its response to invasions. In a randomly assembled community of unrelated strains, the biases are essentially independent up to corrections smaller by O(1/L), with the mean drives, ζi, gaussian distributed with average zero and standard deviation of order 1/L. More precisely, since ξi is determined by the community in the absence of i, the Vij that determine it are independent of the abundances in the community. The variance of ζ is 1/Ljν¯j2 (with only small corrections from the neglect of the strain itself): this simple result is explicit in the self-consistency condition of the DMFT correlation function. We therefore use L as a measure of the effective size of a community which weights the contributions of strains by their mean abundances and, unlike L, is insensitive to whether the close to marginal strains have or have not gone extinct. The distribution of the biases of the persistent strains in the randomly assembled community is a truncated gaussian with lower limit ξc, which is of order 1/L — the magnitude of ξc and L must be determined self-consistently. Crucially, these will depend on the dynamics — especially the blooms — not just the mean quantities.

Appendix 6

Analysis of general fitnesses

Number of persistent strains with exponential P(s) in perfectly antisymmetric model

The antisymmetric model with γ=-1 and no migration (considering only a single island) provides a situation in which we can calculate various quantities analytically and anchor the DMFT analysis. As analyzed in PAF, there is a unique uninvadable fixed point characterized by abundances {νi*} with half the strains having gone extinct, and this fixed point is marginally stable. Indeed, there is a family of chaotic steady states, parametrized a temperature-like quantity Θ which is roughly the range of the fluctuations in logνi, with the average abundances ν¯i=νi*.

Here we analyze the effect of introducing independent exponentially distributed si for each strain, and find the number of strains that persist at the fixed point. Since there is no migration, the relationship between the biases and abundances is especially simple, as detailed below. The critical bias is ξc=0 and the mean s is simply s^=Υ¯. In the limit KΣ21, the distribution of the drives with scale 1/K will be much narrower than that of the si. Let us define xiζi+si for convenience. The scale of p(x), the xi distribution, will be Σ and it will decay exponentially for large x. We can get the distribution of x from the convolution of an exponential and gaussian distribution, with the drives having variance σζ2:

(A6.1) p(x)eσζ2/2Σ2Σex/Σ, for x>σζ(σζΣ+O(1)).

With the Lagrange multiplier Υ¯, the average abundances ν¯i at the mean field fixed point will be ν¯i=max(0,xi-Υ¯X)=max(0,ξi/X) with X the total static susceptibility given by iχi. Let us define ϕ as the fraction of strains that persists with positive bias. When Υ¯>σζ[σζ/Σ+O(1)], we can make the approximation

(A6.2) ϕ=Υ¯p(x)dxeσζ2/2Σ2eΥ¯/Σ.

Self consistency requires that the average abundance over the initial K strains is

(A6.3) Υ¯(xΥ¯X)p(x)dx=1Keσζ2/2Σ2ΣXeΥ¯/Σ1KΣXϕ1K.

Now we combine this relation with the self consistency condition

(A6.4) X=idν¯idξi=KϕXX=Kϕ,

to obtain ϕ1KΣ2. Therefore, the number of persistent strains, when KΣ21 is Σ-2, with coefficient one.

To obtain the behavior for arbitrary K, one must solve the mean field self consistency equations, with the exact distribution

(A6.5) p(x)=1Σexp(σζ22xΣ2Σ2)Φ(xΣσζ2Σσζ),

where Φ is the standard normal cumulative distribution function. There are three mean field equations for the three unknowns Υ¯, σζ and X.

(A6.6) 1K=Υ¯(xΥ¯X)p(x)dx,
(A6.7) σζ2K=Υ¯(xΥ¯X)2p(x)dx,
(A6.8) X=KXΥ¯p(x)dx.

The number of persistent strains is KΥ¯p(x)dx, where Υ¯ solves the above equations. Since for perfectly antisymmetric interactions, s^=Υ¯, the solution of these allows us to check that s^/Σ increases with Σ, as seen in Figure 3. Though one could solve the mean field equations numerically, we can obtain asymptotic results using the approximation in Equation A6.1.

For KΣ21, we have Υ¯Σ(1+log(KΣ2)), X1/Σ and σζ2Σ, which is consistent with the asymptotic Ansatz for p(x) for the present strains: since Υ¯-σζ2/Σσζ all persistent strains have large enough x for the approximation to be valid. Then the distribution of ζ for the persistent strains is

(A6.9) p(ζ|persistent)12πΣe(ζ/2Σ1)2,

valid except for ζ>Υ¯=Σ(1+log(KΣ2)) which is a negligible fraction of the distribution. The width of both the ζ and the s distributions for the persistent strains is of order Σ. As expected more generally, these are comparable once s^ is in the tail of P(s). This calculation could be done for any choice of P(s), and should give the behavior of Υ¯ with K, which is the same as the behavior of s^ with T in the evolving phase of the STC: namely s^Σ(logT)1/ψ.

Joint distribution of ζ and s of successful invaders

Here, the statistical properties of uncorrelated invaders, needed to understand the evolution of the STC state, are analyzed. Every successful invader A must have ξA>ξc. Since the invader’s bias is a sum of sA, ζA and -Υ¯, we can calculate the distributions of sA and ζA conditional on invasion.

We analyze the case with an exponential distribution of the si with scale Σ. Since the distribution of the invader ζA and sA are independent, we can write their joint distribution as

(A6.10) p(ζA,sA)=1ΣL2πexp[sAΣL2ζA2]Θ(sA).

Conditioning on successful invasion multiplies by a factor of Θ(ζA+sA-ξc-Υ¯), which enforces that ξA>ξc. We are interested in the large time limit when Υ¯=s^+O(Σ)Σ. In this limit, the restriction that sA>0 has negligible effect and the analysis is simple.

The invasion probability is

(A6.11) pinvade=0p(ζA,sA)Θ(ζA+sAξcΥ¯)dζAdsAexp[12Σ2LΥ¯+ξcΣ],

which decays exponentially as e-s^/Σ, as expected. The marginal distribution of sA conditioned on invasion is

(A6.12) p(sA|invasion)1ΣpinvadeesA/ΣΦ[L(sAΥ¯ξc)],

where Φ is the standard normal cdf. Similarly, the distribution of ζA conditioned on invasion is

(A6.13) p(ζA|invasion)1pinvadeL2πexp[ζAΣL2ζA2ξc+Υ¯Σ],

and the bias distribution conditioned on invasion is

(A6.14) p(ξA|invasion)1Σexp[(ξAξc)Σ]Θ(ξAξc).

Of particular interest is the mean sA:

(A6.15) E(sA|invasion)Υ¯+ξc+Σ1ΣLs^+O(Σ),

since E(ζA|invasion)1/(ΣL) and E(ξA|invasion)ξc+Σ.

Effects of shape of tail of general fitness distribution

Here, we analyze the evolutionary dynamics for the class of general-fitness distributions with high-s tails parametrized by ψ, and show support from numerical simulations for the scaling Ansatz and results given in the main text.

For attempted invasions by unrelated strains, as the evolution progresses, s^ increases past Σ and pushes up into the tail of P(s). The rate of successful invasions then decreases rapidly since invaders need general fitness (at least) comparable to s^ in order to invade: the probability of successful invasion is of order pinvadeΣ^P(s^) with the prefactor from the width of the distribution of sA near s^. In this regime, the analysis of ‘Appendix 6’ can be carried over with Σ replaced by Σ^. The drive on the invader from the extant community, ζA1/L, will be of order Σ^, yielding the key scaling prediction LΣ^-2 which is tested in Appendix 6—figure 1A.

Appendix 6—figure 1
Consistency of predicted scaling Ansatz with simulations.

(A) The combination LΣ^2 is predicted to approach a constant independent of ψ for long evolutionary times. Solid lines indicate the mean value and shaded regions indicate standard error over 50 replicates, conditional on the diversity not crashing. The curves for ψ=0.8 are shown individually as this value of ψ results in crashing and the fluctuations are large. The dashed horizontal line is at 0.7, the value found for ψ=1 from the data of Figure 3A. Note that for ψ=3, transients are still substantial. (B) Theory predicts that ds^/dZ scales as Σ^/L. The dashed line has slope 1/9 with evolutionary time increasing along the direction of the arrow. Data for ψ=0.8 are quite noisy and not shown. Here the curves are smoothed by a moving average over 1000 successful invasions for both ds^/dZ and Σ^/L, and further averaged over 50 replicates conditional on not crashing. Transients from the initial conditions are observable at the upper right.

Successful invaders will have sA=s^+O(Σ^), since Σ^ is the typical amount by which the invader’s general fitness is likely to be larger (or smaller) than the s^ of the extant community. Therefore, we expect that the change in the community-mean s^ per successful invasion will be ds^/dZΣ^/L (obtainable from ‘Appendix 6’): this scaling is tested in Appendix 6—figure 1B. Substituting for LΣ^-2, we have ds^/dZΣ^3. Integrating this equation and using Σ^=Σψ/s^ψ-1 one obtains scaling laws quoted in the main text, valid for Z1/Σ2 and LZ2ψ-23ψ-2 with the latter exponent less than one for the relevant range of ψ.

For ψ>1, both L and s^ increase sub-linearly with Z, but dL/dZ decreases steadily with Z since the evolution gets closer to an average of one extinction per successful invasion, which obtains exactly in the steady state for the exponential distribution. By contrast, for ψ<1, L decreases with Z and s^ increases super-linearly with Z, as the distribution does not decay fast enough to prevent increasingly fit mutants from emerging and outcompeting more than one strain per successful invasion.

The behaviors of the scaling combination LΣ^2 is shown versus Z in Appendix 6—figure 1A. this combination is predicted to approach a ψ-independent constant for large Z. The value of this constant is consistent with the 0.7 found for ψ=1 in Figure 3A. The other key scaling law, ds^/dZΣ^/L, is tested in Appendix 6—figure 1B, with the predicted coefficient of 1/9 from Figure 3. Although these show convincing evidence for the predicted asymptotic scaling forms, much of the observed evolution is not yet in the late-time asymptotic regime.

To get the growth of s^ — and hence the other quantities — with the number of attempted invasions, we write; ds^/dTpinvadeΣ^/L and integrate this using the scaling relations to get

(A6.16) s^Σ[ψlog(TΣ2)+(33ψ)loglog(TΣ2)+O(1)]1/ψ,

with the unknown O(1) correction reflecting an additive uncertainty of order Σ^ in s^. Up to this correction, the expression has an identical form to the upper limit from pinvadeTL as explained in the main text, and quoted there without the loglog correction (which vanishes for ψ=1). This prediction, with T replaced by T+K (without the loglog part or the O(1) correction), is plotted in Figure 6B. Slow crossovers from the initial non-universal behavior towards the predicted behaviors are observable with the asymptotic predictions thus not well-testable in the time-dependences.

Appendix 7

Correlated general fitness mutations

We have shown that including general fitnesses with ψ1 causes evolution with unrelated invaders to slow down, with Z increasing only as power of logT. Here we analyze what happens if the general fitnesses of mutants are correlated with those of their parent. In order to generate such mutant general fitnesses that are still drawn from the overall P(s), we can take the mutant fitness to be sampled from a Markov chain starting at the parent, with stationary measure given by P(s), and a “time” between samples that depends on the desired correlation between mutant and parent fitnesses.

Exponential distribution of s

The case of exponential P(s) is simple: one can choose the mutant fitness to be gaussian distributed with mean given by the parent’s general fitness minus a constant offset which can be tuned so that the stationary distribution is exponential with the desired Σ. The variance of the gaussian determines the degree to which the parent and mutant are correlated, parametrized by ρs, which need not be the same ρ defined by the correlations between parent and mutant interactions. The crucial difference from uncorrelated s is that the probability of ss^ is now independent of s^ rather than decreasing exponentially with s^/Σ. This enables s^ and Z to increase linearly in T. However, L still saturates to a value of order 1/Σ2, with a coefficient that depends on the correlations. Appendix 7—figure 1 shows simulations that confirm these predictions.

Appendix 7—figure 1
Simulation results for correlated mutants with an exponential P(s).

The primary difference from the independent-invaders case (Figure 3) is that the number of successful invasions is proportional to the number of attempted invasions, since a mutant s is correlated with that of its parent, so invasions do not slow down with increasing s^. Here ρ=0.8 for the interactions and for the general fitnesses the correlation is ρs0.9. Note that for the same Σ with unrelated invaders (ρ=0 and ρs=0) the steady state values of L are quite similar to those shown here.

Non-exponential s distributions with correlated mutants

Incorporating correlations in the general fitnesses of parent and mutant for other distributions of s gives rise to other complications. In the ψ>1 case of particular interest, as s^ grows the mean s-s^ of mutants becomes more and more negative and overwhelms the random part of sM-sP. There is a then a scale of s that diverges as ρs1, such that for s^ above this scale, successful invasions become very unlikely and the diversification slows down even with strongly correlated mutants. This gives rise to complicated crossovers that we do not fully understand and are hard to disentangle in the simulations.

Appendix 8

Dynamics of drive distribution

Markov approximation

In ‘Evolution without general fitness differences’, we describe how the serial invasion of strains causes the biases of extant strains to undergo a random walk with an absorbing boundary condition. In general the statistics of this random walk are complicated and depend on the conditioning on evolutionary history. In order to make analytical progress, we make a Markov approximation of the effects of the evolutionary history and assume that the distribution of the biases at epoch T+1 only depends on itself at epoch T. The average drive of probe strain 0 before the successful addition of strain A is ζ0=jV0jν¯j, while after the invasion it changes to

(A8.1) ζ0=ζ0+δζ0=jV0j(ν¯jA+δν¯j)+V0Aν¯A,

where ν¯A is the mean abundance at which strain A establishes and ν¯jA is the average abundance of strain j before A invades. The sum on j does not include 0 or A and the δ’s denote changes that result from the invasion. Since ν¯A is of order 1/L, its direct effect on the drive is smaller than ζ0 by a factor of 1/L. However A will also change all of the other biases by similar-magnitude random amounts with uncorrelated signs, causing δν¯jχjδζj. Thus E[δζ02]=ν¯A2+jV0j2χj2δζj2 with the expectation over the interactions with the probe strain. Since the properties of the extant strains depend neither on the V0j nor the perturbations via VjA, the averages over each of the squared factors can be performed separately, leaving jχj2 times the average over j of δζj2. But now we can consider each strain separately to be the probe strain so that the average over the extant strains must be the same on both sides. This yields Estrains[δζ02]=ν¯A2/(1-jχj2). The more directly measurable quantity is the total squared abundance change of the extant strains caused by a random change in growth rate of each. With the perturbation caused by the invader, we have j(δν¯j)2=Ξν¯A2 in terms of the nonlinear response, which we call the fragility: Ξ=jχj2/(1-jχj2). This quantity is analogous to the spin-glass susceptibility of random magnets Fisher and Huse, 1988. The fragility is of order unity and the form of the denominator shows that it can diverge: such divergence indicates an instability of the community and breakdown of the DMFT Ansatz. Note that the fragility is a general measure of the sensitivity of a community to perturbations. In the stable niche-phase of the Lotka Volterra model with large negative diagonal interactions of magnitude Q=-E[Vii], the fragility diverges as K increases to the stability boundary at KcritQ2, indicating an instability in the mean field solution Bunin, 2017. The fragility is infinite in the perfectly antisymmetric model (where Q=0), corresponding to being exactly marginal and highly sensitive to added strains.

In addition to the random change of a ζj induced by a successful invader, there is also a systematic change, as seen in Figure 4. This is because δζ0 and ζ0 involve the same set of random V0j and are thus correlated. We can use the fact that E[X|Y=y]=Cov(X,Y)Var(Y)y for zero mean gaussian random variables X and Y to see that the conditional expectation is E[δζ0|ζ0]=ζ0jν¯jδν¯j/jν¯j2. Since each ν¯j and δν¯j are correlated as a consequence of ζj and δζj being correlated, we again have to self-consistently determine this correlation. Unfortunately this is more complicated as the δν¯j also have a contribution from small changes in the function N(ζ) caused by the invaded strain. However one can understand the form of the conditional expectations by noting that after the new strain has invaded,

(A8.2) L(T+1)1=L(T)1+j(δν¯j)2+2jν¯jδν¯j.

This implies that the last term must be negative if L and L are increasing or staying constant. Therefore the correlation between ζj and δζj caused by the invader strain is negative (this is true also if only considering the direct effect of A on j). The additional effects on j from the changes in other strains is enhanced by the fragility, Ξ. If the fragility is high, one can show that 2jν¯jδν¯j is proportional to -ΞE[νA2], since the last two terms in Equation A8.2 should be of similar magnitude to preserve L: in this limit, as argued below, the diversity will decrease.

The above Markovian analysis is an approximation because δζj for each extant strain is correlated with its whole past trajectory — as it involves many of the same V0j — and its conditional expectation will be a weighted sum over the full past ζj(T). Furthermore, it is possible that conditioning on all this past could suppress E[δζj2] by a substantial factor (in the Markovian approximation the analogous correction is smaller by a 1/L factor). However the basic structure of the stochastic changes is correctly captured by the Markovian approximation. This approximation is summarized in Equation 3, which writes the change in bias after one invasion attempt as the sum of a systematic and stochastic part.

If neither extinctions nor invasions occurred, there would be a trivial steady state with a gaussian distribution of the biases with width 1/L as expected for an assembled community. With a source of new strains, which come in with drives that are gaussian distributed with width 1/L, the effective community size is LCL with C<1. We thence obtain a Fokker Planck equation for the number density N of the drives:

(A8.3) TN(ζ,T)=Bζ[ζLN(ζ,T)]+DL2ζ2N(ζ,T)+CL2π eζ2CL/2Θ(ζξc(T)Υ¯(T)),

where we impose the absorbing boundary condition at the critical bias corresponding to N(ξc+Υ¯,T)=0 with ξc and Υ¯ both scaling as 1/L. The last term in Equation A8.3 is the truncated-gaussian drive distribution of successfully invading strains with scale 1/L. The fraction of invasions that are successful, dZ/dT, is just the integral of the truncated gaussian, 1-Φ[L(ξc+Υ¯)]. The number of strains, L(T), is changing and at this point unknown: it is determined self-consistently from L(T)=dζN(ζ,T). It is straightforward to show, as we do in ‘Appendix 8’, that this Fokker-Planck equation admits a scaling solution with the Ansatz that all the quantities scale as 1/L. The rate of diversification, U=dL/dT, is of order one and determined by an eigenvalue-like condition: it can be either positive or negative depending on ξc+Υ¯ and the coefficients, B,D,C. If the community is very fragile (Ξ large), then the B and D terms will dominate over the input and the loss of strains per successful invasion will be large: in this regime the diversification rate is negative. On the other hand if -(ξc+Υ¯)L is large, few strains will go extinct and the community will diversify. Of course, these quantities are determined self-consistently, depending on both the distribution of biases and the function N(ξ), which itself will evolve, reaching a scaling form in the steadily diversifying state.

How good is the Markov approximation? The average change in the drive should really be a weighted integral over the history of the drive, ζ(T), over all T<T. If one makes the Ansatz of a scaling solution with L steadily increasing, then for large T, LT and in the variables scaled by T the weighting function should be a function of of q=log(T/T), so that the integral over the history is a convolution in q. In general, the steady state could, of course, not be found even if this weighting function were known, and the integro-differential generalization of the Fokker-Planck equation would have to be analyzed numerically. But with the widths of the distributions of the extant and invader biases being similar, as observed, the details of the weighting function might not much affect the bias distribution.

Fokker-Planck equation without general fitness differences

Without general fitness differences or correlations, the evolution of the average drive of a strain in the Markovian approximation obeys a Langevin equation (Equation 3). With a source term corresponding to incoming strains, this yields a Fokker Planck equation (Equation A8.3) for the number density N(ζ,T) whose integral at any time gives us the number of strains. The crucial Ansatz is that all the quantities (ζ, ξc and Υ¯) scale as 1/L, and we now take L to be a function of T. Then we can hypothesize a scaling solution N(ζ,T)L3/2g(ζL) for some function g so that N(ζ,T)dζL(T), as the integral over the distribution of drives gives us the total number of strains. For ease of notation, we define the similarity variable uζL, and derivatives with respect to u by primes: then Equation A8.3 becomes

(A8.4) dLdT(32g(u)+12ug(u))=Dg(u)+C2πeCu2/2Θ(uu0)+Bg(u)+Bug(u),

where we define u0L(ξc+Υ¯). We can make an Ansatz of the form

(A8.5) L(T)=UT+L0,

where the coefficient U is the average rate of increase or decrease of the diversity per invasion attempt.

The equation for g(u) then becomes

(A8.6) (B3U2)g+(BU2)ug+Dg+C2πeCu2/2Θ(uu0)=0.

We can solve this equation numerically with boundary conditions g(u0)=0 and g()=0 (exact solutions are available for certain values of B,D and U). The first of these boundary conditions comes from the extinction criterion that enforces ζi>ξc+Υ¯ for every extant strain, and second comes from the need to have L<. The solution depends on the parameters u0, B, D and U. By enforcing u0g(u)du=1, we can find the value of U in terms of u0,B and D. One can see, numerically, that for B and D both large, the consistent value of U is negative, implying loss of diversity.

In Figure 5B, we use the measured values of U and C, and adjust the values of B and D to give a normalized function g(u). The blue curve is a truncated unit-variance gaussian for the initial assembled community, with mean and lower limit fit by hand, and the orange curve is a numerical solution to Equation A8.6 with U=dL/dT0.24 and B=2.11, D=1.79 adjusted to enforce normalization. Note that C=L/L0.74. The mode of the truncated gaussian is Υ¯L as expected from Figure 4. The bias distributions for assembled and evolved communities are very similar except near ξc, which is what one expects from the Markov approximation. For -ζc=-ξc-Υ1/L, when the ratio of extant to invading widths, DC/B0.79, is close to unity, the bulk of the drive distribution is close to gaussian with mean zero, which is reflected in Figure 5B.

Fokker-Planck equation with independent general fitness differences

If incoming strains have independently drawn general fitnesses, si, then we can write an evolution equation for the joint distribution of the drives and the si:

(A8.7) ZN(ζ,s,Z)=Bζ(ζLN)+DL2ζ2N+P>(s)CL2πeζ2CL/2Θ(ζ+sξcΥ¯),

where we abbreviate P(s|s>ξc+Υ¯ζ) by P>(s). However, unlike the case with general fitness differences, now Υ¯ does not scale as 1/L. Indeed, the scaling of Υ¯ with Z determines how L depends on Z. If we specify the form of dΥ¯dZ, then we can look for solutions which travel up in the s direction at the same speed as Υ¯.

We can define similarity variables uζL and vL(s-Υ¯). Then we look for solutions of the form N(ζ,s,Z)L2g(u,v) so that the integral of N over ζ and s is proportional to L. The boundary condition is that g vanishes along the line u+v=ucrit. Plugging this into our Fokker Planck equation, and defining ucritLξc, we obtain

(A8.8) dLdZ(2g+u2gu+v2gv)dΥ¯dZL3/2gv=Dguu+C2πP>(v)eCu2/2Θ(u+vucrit)+Bg+Bugu,

where subscripts denote derivatives of g. We can see that the P>(s) distribution generically breaks the scaling — because it depends on s=v/L, not only on v. But we can look for solutions which obey the scaling with 1/L. Take v/L as exponentially distributed with scale Σ. Then in order for our scaling Ansatz to be valid, the quantity vΣL must be independent of L, so we must have L1/ΣLΣ2. This requires dΥ¯dZΣ3 in order to make the left hand side of Equation A8.8 independent of L. Furthermore, we have dL/dZ0 since LΣ-2. Therefore, we have

(A8.9) dLdZ=gv+Dguu+C2πP>(v)eCu2/2Θ(u+vucrit)+Bg+Bugu2g+u2gu+v2gv.

Setting this to 0 would allow us to solve for g(u,v) and therefore find the scaling solution for the joint distribution of mean drives and the si.

Appendix 9

Coexistence of correlated strains

The simulations show that even strains which interact very similarly with the rest of the community can coexist. Here, we present some of the numerical data and show how simple approximations fail to explain this behavior.

Lotka-Volterra approximation ignoring migration

A natural way to understand very closely related strains (e.g. a mutant and parent) is via an effective Lotka Volterra model for the dynamics of just these two strains in the presence of all the others. For a correlation ρ between the parent and mutant interactions, one can obtain a pair of DMFT equations describing the correlated dynamics of parent and mutant in the form of Equation A5.1. Here, we consider the crude approximation of neglecting the effects of migration and averaging the dynamics of logν to obtain equations for parent and mutant frequencies νP and νM:

(A9.1) νP˙=νP[ξP-(-γX)(νP+ρνM)]
(A9.2) νM˙=νM[ξM-(-γX)(νM+ρνP)].

This analysis is appropriate only for strains with positive bias — and even then approximate – but the following results provide a null model against which one can compare the dynamics of parent-mutant replacement in the evolving STC.

In the absence of the mutant, ν¯PξP/(-γX) (which is the correct behavior for large positive ξP). The condition for invasion of a mutant into an environment containing only the parent is that its invasion eigenvalue, suppressed by the effects of the parent, is greater than the critical bias which amounts to ξMρξP>ξc and the condition for invasion of the parent into the mutant is likewise ξPρξM>ξc. In terms of the drives, these conditions are ζMρζP>(1ρ)Υ¯+ξc for the mutant to invade the parent, and the same condition with M and P flipped for the parent to invade the mutant. With the mutant drive parameterized as ζM=ρζP+1-ρ2Z1 with si a gaussian random variable with mean 0 and scale ζP, the probabilities of coexistence between parent and mutant, and of replacement of the parent, are

(A9.3) pcoexist=Φ[(1ξ~c)/ρρ+Υ~(11/ρ)1ρ2]Φ[Υ~(1ρ)+ξ~c1ρ2]
(A9.4) preplace=1Φ[(1ξ~c)/ρρ+Υ~(11/ρ)1ρ2]

where Φ is the standard normal cdf, and Υ~ and ξ~c are the Lagrange multiplier and critical bias, both normalized by the scale of the ζ distribution. Note that the expressions are only valid for Υ~<1ξ~c1ρ2: otherwise, in this approximation, the coexistence probability vanishes and all invasions are replacements. The probability of mutant invasion is pinvade=1-Φ[Υ~(1-ρ)+ξ~c1-ρ2]. For 1-ρ1 and ξ~c=0, the invasion probability is pinvade12(1-Υ~(1-ρ)/π) and the coexistence probability predicted by this simple approximation is pcoexist(1-Υ~)(1-ρ)/π. With ρ=0.99 and Υ~=0 (in the STC Υ~ should be nonnegative, and setting it to 0 maximizes the coexistence probability), this predicts that the probability of coexistence conditional on the mutant invading is pcoexist/pinvade0.12 which is much smaller than the observed pcoexist|invade>0.5 (Appendix 9—figure 1).

Appendix 9—figure 1
The coexistence probability between parent and mutant conditioned on successful invasion of the mutant (with si=0).

The probability is averaged over 4 simulations of a diversifying community for each value of ρ, with error bars showing the standard error. Here ρ is shown on a logit scale to emphasize the difference between the points as ρ1. Even for ρ=0.999, the coexistence probability is order 1/2 over he epoch of duration 30ML, likely long enough to enable small differences between closely related strains to have an effect.

The main problem with this approximation is that it needs to be used when the biases are close to the critical ξc which is substantially negative, while the approximation of ignoring migration certainly breaks down in this regime. Properly understanding the coexistence of replacement of related mutant and parent requires a fuller understanding of their coupled dynamics, including migration.

Coexistence probability

Of particular interest is the probability with which mutant and parent coexist when they are correlated by amount ρ. The effects of the persistent chaos on this coexistence probability are subtle, and here we show that the simulation results cannot be captured by the null model of ‘Appendix 9’, which neglects migration and blooms.

Comparing numerical results from simulation of the dynamics to those of ‘Appendix 9’, we find that the probability of coexistence between parent and mutant, conditioned on mutant invasion, is substantially larger than predicted by Equation A9.1 and Equation A9.2. This is likely due to the fact that in the STC, a majority of persistent strains have negative bias and even for those with positive bias, the effect of migration is not small. The boom-bust dynamics with migration makes coexistence much easier, but understanding this even semi-quantitatively is challenging because of the correlations in the dynamical drive which makes parental blooms suppress the mutants and vice versa. Nevertheless, we expect that if ρ is extremely close to unity, the coexistence probability of parent and mutant will go to zero as 1-ρ, as in the simple approximation from ‘Appendix 9’, but with a small coefficient that might be of order some inverse power of M. In addition one should note that the natural timescale over which differences between the parent and mutant will accumulate is of order (1-ρ2)-1/2ML with the coefficient only 20 for ρ=0.999. Therefore even for ρ this close to 1, the typical differences between parent and mutant will still be felt on timescales of order ML over which the close-to-marginal strains rise or die out.

A dynamical analysis of a parent and mutant including migration becomes far more complicated once the community has evolved long enough that a substantial fraction of the strains have turned over (or coexist with relatives). The underlying simplification that makes DMFT valid for assembled communities, the approximate independence of the biases, will break down because of correlations among many of the strains due to their common ancestry. While some of the heuristic behavior of such communities and how they evolve may not change much even for ρ quite close to 1, there will certainly be quantitative changes.

The approximate Markovian analysis of the evolution of the bias distribution (‘Evolution without general fitness differences’), can be modified to roughly account for correlations between mutants and parents, ρ>0. Such correlations will modify the effects of invading strains on the bias of the extant strains. This is most readily understood in the limit that ρ is very close to unity. In this case, one expects that the total abundance of the parent and mutant after the invasion will be similar to that of the parent before the invasion and, because of the high correlations, their effects on the other strains will be very similar. This suggests that both the stochastic and systematic changes of the biases should be multiplied by a factor of order 1-ρ2 on the right hand side of Equation 3. For the first few mutants, the analysis could be carried through similarly. However once the number of successfully invaded mutants is comparable to the original L0, the correlations between the biases of the extant strains can not be ignored and the approximation of independent ζi breaks down. We leave analysis of this for future work.

Data availability

The current manuscript is a computational study, so no data have been generated for this manuscript. Simulations use only standard algorithms: details in paper.


  1. Book
    1. MacArthur RH
    2. Wilson EO
    The Theory of Island Biogeography
    Princeton University Press.
    1. Rocco A
    2. Ebert U
    3. Saarloos W
    (2000) Subdiffusive fluctuations of "pulled" fronts with multiplicative noise
    Physical Review. E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 62:R13–R16.

Article and author information

Author details

  1. Aditya Mahadevan

    Department of Physics, Stanford University, Stanford, United States
    Software, Formal analysis, Investigation, Visualization, Writing - original draft
    Contributed equally with
    Michael T Pearce
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0009-0000-5571-9993
  2. Michael T Pearce

    Department of Physics, Stanford University, Stanford, United States
    Present address
    Meta Data Science, Menlo Park, United States
    Software, Formal analysis, Investigation, Visualization, Methodology
    Contributed equally with
    Aditya Mahadevan
    Competing interests
    No competing interests declared
  3. Daniel S Fisher

    Department of Applied Physics, Stanford University, Stanford, United States
    Conceptualization, Formal analysis, Supervision, Writing – review and editing
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-5559-2491


National Science Foundation (PHY-160760)

  • Aditya Mahadevan

National Institutes of Health (R01AI13699201)

  • Aditya Mahadevan

Simons Foundation (Sabbatical Fellowship)

  • Daniel S Fisher

National Science Foundation (PHY-2210386)

  • Aditya Mahadevan

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.


We thank Pankaj Mehta for useful discussions and comments on the manuscript. This work was supported in part by the National Science Foundation via PHY-1607606,and PHY-2210386, the National Institutes of Health via R01AI13699201, the Simons Foundation via a sabbatical Fellowship to DSF, and the Stanford University Sherlock Cluster for the computations.

Version history

  1. Preprint posted: May 28, 2022 (view preprint)
  2. Received: August 16, 2022
  3. Accepted: April 27, 2023
  4. Accepted Manuscript published: April 28, 2023 (version 1)
  5. Version of Record published: June 28, 2023 (version 2)


© 2023, Mahadevan, Pearce et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.


  • 1,002
  • 221
  • 7

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Aditya Mahadevan
  2. Michael T Pearce
  3. Daniel S Fisher
Spatiotemporal ecological chaos enables gradual evolutionary diversification without niches or tradeoffs
eLife 12:e82734.

Share this article

Further reading

    1. Ecology
    Xueyou Li, William V Bleisch ... Xue-Long Jiang
    Research Article

    Spatial and temporal associations between sympatric species underpin biotic interactions, structure ecological assemblages, and sustain ecosystem functioning and stability. However, the resilience of interspecific spatiotemporal associations to human activity remains poorly understood, particularly in mountain forests where anthropogenic impacts are often pervasive. Here, we applied context-dependent Joint Species Distribution Models to a systematic camera-trap survey dataset from a global biodiversity hotspot in eastern Himalayas to understand how prominent human activities in mountain forests influence species associations within terrestrial mammal communities. We obtained 10,388 independent detections of 17 focal species (12 carnivores and five ungulates) from 322 stations over 43,163 camera days of effort. We identified a higher incidence of positive associations in habitats with higher levels of human modification (87%) and human presence (83%) compared to those located in habitats with lower human modification (64%) and human presence (65%) levels. We also detected a significant reduction of pairwise encounter time at increasing levels of human disturbance, corresponding to more frequent encounters between pairs of species. Our findings indicate that human activities can push mammals together into more frequent encounters and associations, which likely influences the coexistence and persistence of wildlife, with potential far-ranging ecological consequences.

    1. Ecology
    Lan Pang, Gangqi Fang ... Jianhua Huang
    Research Article

    The success of an organism depends on the molecular and ecological adaptations that promote its beneficial fitness. Parasitoids are valuable biocontrol agents for successfully managing agricultural pests, and they have evolved diversified strategies to adapt to both the physiological condition of hosts and the competition of other parasitoids. Here, we deconstructed the parasitic strategies in a highly successful parasitoid, Trichopria drosophilae, which parasitizes a broad range of Drosophila hosts, including the globally invasive species D. suzukii. We found that T. drosophilae had developed specialized venom proteins that arrest host development to obtain more nutrients via secreting tissue inhibitors of metalloproteinases (TIMPs), as well as a unique type of cell—teratocytes—that digest host tissues for feeding by releasing trypsin proteins. In addition to the molecular adaptations that optimize nutritional uptake, this pupal parasitoid has evolved ecologically adaptive strategies including the conditional tolerance of intraspecific competition to enhance parasitic success in older hosts and the obligate avoidance of interspecific competition with larval parasitoids. Our study not only demystifies how parasitoids weaponize themselves to colonize formidable hosts but also provided empirical evidence of the intricate coordination between the molecular and ecological adaptations that drive evolutionary success.