Spatiotemporal ecological chaos enables gradual evolutionary diversification without niches or tradeoffs
Abstract
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. Hostpathogenlike interactions drive the community into a spatiotemporally chaotic state characterized by continual, spatiallylocal, 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 dynamicalmean fieldtheory 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 finescale 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.
https://doi.org/10.7554/eLife.82734.sa0Introduction
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 (BonillaRosso 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 “microniches” 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 piconiches, 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 multiniche 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 (IgnacioEspinoza et al., 2020; MartinPlatero 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 closetoperfect 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; RodriguezValera 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 withinspecies diversity. Understanding this process theoretically is a longterm goal, towards which the present work is a step. To make progress, we need to distill the general phenomenon of finescale 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 nichelike 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 wellmixed 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 oneonone: 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 bacteriaphage diversity, to build up an understanding of the complex ecoevolutionary dynamics, we focus in this paper on simpler models that — as we have shown in PAF — capture many of the key features.
Hostpathogen, 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 commonlyinvoked 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.
Complementary work (Roy et al., 2020; Roy et al., 2019) suggests the generality of the STC beyond anticorrelated pairwise interactions, although in the LotkaVolterra models studied in these works of Roy et al., the diversity is limited by the strength of selfinteractions, 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 interspecies interactions (Bunin, 2017). Here we approach evolutionary dynamics in similar generalized LotkaVolterra models, but from the opposite starting point: all interactions are of comparable magnitude which makes the effects of selfinteractions 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 closetorigid 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 consumerresource 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 ecoevolutionary 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, boombust 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 contextdependent — 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 $\rho \in [0,1)$. Many of the behaviors are similar across this range, from independent invaders with $\rho =0$ to smalleffect mutations with $\rho $ 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 $\rho $ (and expect for any $\rho <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 $s}_{i$ for strain $i$) can either slow down, prevent, or reverse diversification. For distributions of the $s}_{i$ 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 $\mathrm{log}(T)$. If the distribution of the $s}_{i$ has a broaderthanexponential 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 lowabundance 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 ${\xi}_{i}$, has an intrinsic contribution from its general fitness $s}_{i$, 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 ($\rho =0$) with no general fitnesses (${s}_{i}=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/\sqrt{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 $s}_{i$, 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 allatonce from unrelated strains. At early stages of the evolution, most of the closetomarginal, lowabundance 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.
Outline
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.
Models
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.
Ecological interactions
We first consider an assembled community of $K$ unrelated strains, labelled by $i=1,2,\mathrm{\dots},K$, with all possible pairwise interactions between them. A paradigmatic model for the ecological dynamics of the strain populations $\{{n}_{i}\}$ is the generalized LotkaVolterra 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 ${W}_{ij}$ 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 $\mathrm{{\rm Y}}(t)$ that fixes the total population to ${\sum}_{i}{n}_{i}=N$, and work with fractional abundances, ${\nu}_{i}\equiv {n}_{i}/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, $\{{s}_{i}\}$, between the strains. We parameterize the residual variations in interactions among the strains (after subtracting off $\mathrm{{\rm Y}}$ and $s}_{i$) by ${V}_{ij}$. Since the ${V}_{ij}$ and $s}_{i$ 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 $\mathrm{E}[{V}_{ij}]=0$, and $\mathrm{E}[{V}_{ij}^{2}]=1$ for $i\ne j$ — 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 $\mathrm{E}[{V}_{ij}{V}_{ji}]=\gamma $. For convenience, we choose $\mathrm{E}[{V}_{ii}^{2}]=1+\gamma $ but this choice has negligible effect in large communities.
The parameter $\gamma $ controls whether the interactions are mainly competitive ($\gamma >0$) or hostpathogenlike ($\gamma <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 hostpathogen models with the appropriate block submatrix structure, as discussed further in ‘Bacteriaphage 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
with ${m}_{\alpha \beta}$ the migration rate (per individual) from island $\beta $ to island $\alpha $ and the local Lagrange multiplier, ${\mathrm{{\rm Y}}}_{\alpha}={\sum}_{i}{\nu}_{i,\alpha}({s}_{i}+{\sum}_{j}{V}_{ij}{\nu}_{j,\alpha})$, keeping the total population on each island fixed at $N$, (i.e. ${\sum}_{i}{\nu}_{i,\alpha}=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 $\alpha $ is then simply $m({\overline{\nu}}_{i}{\nu}_{i,\alpha})$ with ${\overline{\nu}}_{i}(t)$ the average of ${\nu}_{i}(t)$ across islands. As the number of islands becomes large, in steady state, each ${\overline{\nu}}_{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, ${\overline{\nu}}_{i}$, is equal to the timeaveraged abundance of strain $i$ on a single island; this is a crucial selfconsistency 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: $\nu}_{\mathit{\text{floor}}}\sim m\overline{\nu$. For strains near local extinction (when the fractional abundance is close to $1/N$) demographic fluctuations are potentially important. But with $Nm\overline{\nu}$ 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 $mN\gg 1$, 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 $\mathrm{log}\nu$ of $\mathcal{M}\equiv \mathrm{log}(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 $\nu \sim \mathcal{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 nearextinctions. On a single island, at any given moment, $\mathcal{O}(L/\mathcal{M})$ strains are at high abundance. A snapshot of the abundances on a single island shows the strains distributed roughly uniformly in $\mathrm{log}(\nu )$ 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 ${V}_{Aj}$ and ${V}_{jA}$, and its general fitness, ${s}_{A}$. 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, $\rho \in [0,1)$, with a parent strain $P$ chosen from the existing community with probability proportional to its mean abundance ${\overline{\nu}}_{P}$. These correlations are realized such that $\mathrm{Corr}[{V}_{jM},{V}_{jP}]=\mathrm{Corr}[{V}_{Mj},{V}_{Pj}]=\rho $ for $j\ne M,P$. The detailed choices for $j=M,P$ are given in ‘Appendix 2’. The general fitness, $s}_{\mathit{\text{M}}$, can also be correlated with $s}_{\mathit{\text{P}}$. Unrelated invaders are equivalent to $\rho =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 ${\sum}_{i}{\nu}_{i,\alpha}=1$).
Timescales
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, ${\zeta}_{i}$, which is of order $1/\sqrt{L}$, implying that the timescale for systematic population growth or decay is of order $\sqrt{L}$. The mean drive is defined more precisely in Results. When there are general fitness differences, $s}_{i$, these also contribute to variations in average growth rates. The variations in the $s}_{i$ within the community have substantial effects over a time $\sim 1/{\sigma}_{s}$ with ${\sigma}_{s}$ roughly the width of the $s}_{i$ distribution of the extant strains. Together, the mean drive and general fitness of a strain determine its crucial property: the bias ${\xi}_{i}={\zeta}_{i}+{s}_{i}\u27e8\mathrm{{\rm Y}}\u27e9$, 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 interstrain variation in general fitness is no larger than variation in average drive from interactions. This means that the biases are of order $1/\sqrt{L}$.
The local population of each strain undergoes wild fluctuations over a logarithmic range $\mathcal{M}=\mathrm{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 $\mathcal{M}\sqrt{L}$ 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 $\mathcal{M}L$ 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 $3\mathcal{M}L$, and we show in ‘Appendix 2’ that increasing this epoch length by a factor of 10 makes little difference in the diversification dynamics.
Results
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 ${\xi}_{i}={s}_{i}+{\zeta}_{i}\u27e8{\mathrm{{\rm Y}}}_{\setminus i}\u27e9$ with the mean drive ${\zeta}_{i}={V}_{ij}\u27e8{\nu}_{j\setminus i}\u27e9$, where the notation $\u27e8{\nu}_{j\setminus i}\u27e9$ and $\u27e8{\mathrm{{\rm Y}}}_{\setminus i}\u27e9$ denote the timeaveraged abundance of strain $j$ and average Lagrange multiplier in the absence of strain $i$. The $\mathrm{{\rm Y}}$ 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 $\gamma $ — 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 $\overline{\nu}$ instead of time average notation $\u27e8\nu \u27e9$, 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 $\mathrm{log}{\nu}_{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 $\mathrm{log}{\nu}_{i}$ (Figure 1). This is a manifestation of the system “selftuning” to a special selfconsistently chaotic state [PAF].
Despite the possibility of rescue from extinction via rare blooms, there is a critical negative bias, ${\xi}_{\mathrm{c}}$, (sharp for large $L$ and large $I$) below which strains no longer persist even as $I\to \mathrm{\infty}$. For strains with $\xi $ below ${\xi}_{\mathrm{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 $\xi <{\xi}_{\mathrm{c}}$ go extinct, while strains with $\xi >{\xi}_{\mathrm{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 $L}_{0$. The $K{L}_{0}$ 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 ${\xi}_{\mathrm{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 $\mathcal{M}L$. 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 $\gamma =0.8$ and unrelated invaders ($\rho =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’.
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 strainlevel diversification up to arbitrarily high number of strains. The behavior depends on the symmetry parameter $\gamma $, which must be substantially negative for the STC to exist. Figure 2B shows that the average rate of diversification, $U\equiv \u27e8dL/dT\u27e9$, is nonmonotonic in $\gamma $, with slow diversification close to $\gamma =1$ (its lower limit). As $\gamma $ becomes less negative the rate of diversification increases at first. However for $\gamma $ even less negative, the STC still supports chaotic coexistence of many strains (since $L}_{0$ is still large), but the diversity decreases under evolutionary dynamics. The community diversifies most rapidly for $\gamma \approx 0.8$. As we are interested in what can happen with various other additional features, we chose $\gamma =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 $\gamma $ 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 $\rho $ 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 $\rho $ (with $\gamma =0.8$), the rate of diversification is nonmonotonic being fastest for $\rho \approx 0.8$, and only weakly varying for smaller $\rho $ (Figure 2C). As $\rho $ 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 $\rho $ slows the rate of diversification but $L(T)$ still increases linearly.
This observation implies that, for a large range of $\rho $, 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 $\rho $ 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 $s}_{i$ between strains and show how these affect the evolutionary dynamics.
Exponential distribution of $s}_{i$
We first analyze the simplest case: exponentially distributed selective differences with scale $\mathrm{\Sigma}$ and probability density $\mathcal{P}(s)=\frac{1}{\mathrm{\Sigma}}{e}^{s/\mathrm{\Sigma}}\mathrm{\Theta}(s)$ where $\mathrm{\Theta}(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 $s}_{i$ of the community will change, and we are particularly interested in the dynamics of the populationweighted mean $\widehat{s}={\sum}_{j}{\overline{\nu}}_{j}{s}_{j}$.
The width of the distribution of $s}_{i$, here $\mathrm{\Sigma}$, 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 $s}_{i$ is likely inconsistent with a diverse community. We therefore focus on narrow distributions: that is small $\mathrm{\Sigma}$. The typical magnitude of the drive of a strain is of order $1/\sqrt{L}$; therefore, when $\mathrm{\Sigma}$ is much smaller than this, it will not matter much. On the other hand, if $\mathrm{\Sigma}$ were to be much larger than $1/\sqrt{L}$, the differences in the $s}_{i$ would dominate over the drives and only the strains with the highest and quite similar $s}_{i$ would survive. Thus $L\gg 1/{\mathrm{\Sigma}}^{2}$ seems inconsistent. Even for the initial community with $L}_{0$ strains, we expect that $L}_{0$ cannot be larger than order $1/{\mathrm{\Sigma}}^{2}$ (although it can be much smaller if $K\ll 1/{\mathrm{\Sigma}}^{2}$). Indeed, in ‘Appendix 6’ we show that $\mathrm{\Sigma}$ sets the initial persistent community size, $L}_{0$, in a particular limit of the model where $\gamma =1$.
A natural conjecture is that for small $\mathrm{\Sigma}$ with an exponential distribution, steady diversification can occur until the breadth of the $s$ distribution becomes important — when $L\sim 1/{\mathrm{\Sigma}}^{2}$ — and after that $L$ will saturate, as seen in Figure 3. Thereafter, $\widehat{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 $\mathrm{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 $L}_{0$ when $K\gg 1/{\mathrm{\Sigma}}^{2}$, the $\mathrm{\Sigma}$dependence of the steadystate $L$, and the linear increase of $\widehat{s}$ with number of successful invasions.
More general distributions of $s}_{i$
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:
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 $s}_{i$, all $s}_{i$ can be shifted by a constant without affecting the dynamics because this constant gets absorbed into $\mathrm{{\rm Y}}(t)$.
As we will analyze in ‘Diversification rate with a distribution of general fitnesses’, the evolution of diversity is seen to depend crucially on $\psi $. If the tail of the $s$ distribution falls off faster than a simple exponential, $\psi >1$, the community continually diversifies, albeit more and more slowly with $L(T)$ increasing only as a power of $\mathrm{log}(T)$. Concomitantly, the mean $s$ of the community, $\widehat{s}$, gradually increases. But if the $s$ distribution decays slower than a simple exponential, $\psi <1$, the diversity decreases (after an initial increase if $\mathrm{\Sigma}$ is sufficiently small) and eventually crashes. In the marginal case of a simple exponential tail, $\psi =1$, as seen above, the diversity saturates and fluctuates around a steady state value while the mean $\widehat{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 “diminishingreturns 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 $s}_{i$ ($\psi =1$) analysis suggests that evolution proceeds at a constant rate, with both $Z(T)$ and $\widehat{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 $\mathcal{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 $\widehat{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 $\mathcal{P}(s)$ that decays faster than exponentially, $\psi >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 diminishingreturns epistasis.
Analysis
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, $\{{\xi}_{i}\}$, and how these set their mean abundances, $\{{\overline{\nu}}_{i}\}$. For a large randomly assembled or evolved community the mean abundances will be a function of the biases: ${\overline{\nu}}_{i}\approx \mathcal{N}({\xi}_{i})$, with the function $\mathcal{N}$ depending on the parameters, evolutionary history, and feedback from other strains. As shown in PAF, $\mathcal{N}(\xi )$ is linear for large argument and decays as $\xi $ becomes negative, vanishing at $\xi ={\xi}_{\mathrm{c}}<0$.
The relation between ${\overline{\nu}}_{i}$, ${\xi}_{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 semiquantitative results. Of particular importance is the result that in evolving communities, the density of biases vanishes linearly as $\xi \to {\xi}_{\mathrm{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 $\mathcal{P}(s)$ to generalize to other shapes of the tail of $\mathcal{P}(s)$, parametrized by $\psi $ (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 selfconsistently 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, ${\sum}_{j}{V}_{ij}{\nu}_{j}$, are comprised of roughly independent random variables and thus act like gaussian noise with correlations $C(t,{t}^{\prime})={\sum}_{j}{\nu}_{j}(t){\nu}_{j}({t}^{\prime})$. However, when it rises to substantial abundance, it will weakly affect the other strains. Because of the correlation $\gamma $ between ${V}_{ij}$ and ${V}_{ji}$, these feedback terms add coherently, resulting in a contribution to growth rate of $i$ of form $\gamma {\int}_{0}^{t}R(t,{t}^{\mathrm{\prime}}){\nu}_{i}({t}^{\mathrm{\prime}})d{t}^{\mathrm{\prime}}$, 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 singlestrain 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 ${\xi}_{\mathrm{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 $\mathcal{L}=1/{\sum}_{i}{\overline{\nu}}_{i}^{2}$ which controls the variances of the meandrive part of the bias; $\mathcal{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 ${\xi}_{A}={\sum}_{j}{V}_{Aj}{\overline{\nu}}_{j}\overline{\mathrm{{\rm Y}}}$, in terms of its interactions, $\{{V}_{Aj}\}$ with the extant strains: for unrelated invaders, this is gaussian distributed with mean $\overline{\mathrm{{\rm Y}}}$ and standard deviation of $1/\sqrt{\mathcal{L}}$ independent of correlations among the extant strains, though correlations in the existing community will affect the $\mathcal{N}(\xi )$ and hence $\mathcal{L}$.
Strain $A$ can successfully invade the community and persist if $\xi}_{A}>{\xi}_{\mathrm{c}$. The probability of successful invasion is thus $\mathrm{\Phi}[({\xi}_{\mathrm{c}}+\overline{\mathrm{{\rm Y}}})\sqrt{\mathcal{L}}]$, with $\mathrm{\Phi}$ the standard normal cumulative distribution function. In the initial assembled community, ${\xi}_{\mathrm{c}}\sqrt{\mathcal{L}}$ and $\overline{\mathrm{{\rm Y}}}\sqrt{\mathcal{L}}$ are independent of $\mathcal{L}$. We make the Ansatz that after a long period of evolution the distribution of extant biases, scaled by $\sqrt{\mathcal{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 onebyone 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 ${\xi}_{\mathrm{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 $\xi $ being pushed below ${\xi}_{\mathrm{c}}$ by an invading strain. For finite $L$, the sharpness of ${\xi}_{\mathrm{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.
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/\sqrt{\mathcal{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 $\overline{\mathrm{{\rm Y}}}\sqrt{\mathcal{L}}$. (Figure 4B). This is likely due to conditioning on survival of the strains: those that persist for many epochs tend to have largerthanaverage (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, ${\overline{\nu}}_{i}$ (when ${s}_{i}=0$, this is just ${\xi}_{i}+\mathrm{{\rm Y}}$). 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 $\delta {\zeta}_{i}\sim \pm 1/L$, but the systematic change in the drive is smaller and depends on its current value: $\mathrm{E}[\delta {\zeta}_{i}{\zeta}_{i}]\sim {\zeta}_{i}/L=\mathcal{O}(1/{L}^{3/2})$. In the Markovian approximation, there is a simple Langevin equation for the change of the drive of a strain due to invasions:
with ${\eta}_{i}(T)$ approximately gaussian with mean zero and unit variance — an approximation that should be good if one coarsegrains 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 orderunity coefficients which respectively characterize the average and meansquared response of the bias to the invasion of a new strain. Both are proportional to $\mathrm{E}[{\nu}_{A}^{2}]$ times the fragility, $\mathrm{\Xi}$, of the extant community which is given by $\mathrm{\Xi}={\sum}_{j}{\chi}_{j}^{2}/(1{\sum}_{j}{\chi}_{j}^{2})$ in terms of the individual susceptibilities of strains to changes in their biases, ${\chi}_{i}=d{\overline{\nu}}_{i}/d{\xi}_{i}$. This fragility characterizes the meansquare 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 $\xi =\zeta \mathrm{{\rm Y}}$ goes below ${\xi}_{\mathrm{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 eigenvaluelike 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 ${\xi}_{\mathrm{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 $\zeta \to {\xi}_{\mathrm{c}}+\overline{\mathrm{{\rm Y}}}$ (or, equivalently, as $\xi \to {\xi}_{\mathrm{c}}$).
Simulations confirm the expectation that the typical bias scales as $1/\sqrt{L}$ over one order of magnitude in $L$ (Appendix 3—figure 4). As predicted, one observes a smoothing out of the bias distribution toward ${\xi}_{\mathrm{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.
However, the density of biases near ${\xi}_{\mathrm{c}}$ determines the response of the community to evolutionary perturbations, since these lowbias strains are the ones most susceptible to extinction. In particular, the predicted linearly vanishing density of biases determines the distribution of the number, $\mathrm{\ell}$, of extinctions per successful invasion (Appendix 3—figure 1B). To estimate this distribution — particular the probability that $\mathrm{\ell}$ 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, ${\overline{\nu}}_{A}$, is gaussian, since for positive ${\xi}_{\text{inv}}$, ${\overline{\nu}}_{\text{inv}}\sim {\xi}_{\text{inv}}$ and the invader’s bias ${\xi}_{\text{inv}}$ is itself gaussian distributed. The number of strains whose biases are within ${\overline{\nu}}_{A}$ of ${\xi}_{\mathrm{c}}$ is proportional to ${L}^{2}{\overline{\nu}}_{A}^{2}$, because the distribution of strains' biases vanishes linearly at ${\xi}_{\mathrm{c}}$. Thus for fixed ${\overline{\nu}}_{A}$, the number of strains that are driven extinct is Poisson distributed with mean proportional to ${L}^{2}{\overline{\nu}}_{A}^{2}$: this is of order one for large $L$ as expected. That the tail of the distribution of ${\overline{\nu}}_{A}$ is gaussian implies ${\overline{\nu}}_{A}^{2}$ is approximately exponentially distributed in its tail. Integrating the Poisson distribution over this yields, for large $\mathrm{\ell}$, $\mathcal{P}(\ell )\sim {e}^{\beta \ell}/\sqrt{\ell}$ (with $\beta $ an orderunity coefficient) which is close to exponential as observed in Appendix 3—figure 1B.
A similar analysis for the initial randomlyassembled community shows that for fixed ${\overline{\nu}}_{A}$, the mean number of extinctions triggered by the first successful invasion is of order ${L}_{0}^{3/2}{\overline{\nu}}_{A}\sim \sqrt{{L}_{0}}$; 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 ${\overline{\nu}}_{A}$ and hence itself be roughly gaussian, though unless the initial $L}_{0$ 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 $L}_{0$ strains — with a relatively small coefficient — to go extinct before $L$ starts steadily increasing, and this will occur over of order $L}_{0$ invasions. For $\gamma =0.8$ this effect appears to be very small — the critical bias is quite negative — but for smaller or larger $\gamma $ the effects are noticeable (Figure 2).
The distribution of mean abundances, ${\overline{\nu}}_{i}$, is related to that of the biases via the function $\mathcal{N}(\xi )$: 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 closetomarginal 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 ${\xi}_{\mathrm{c}}$ due to invasiontriggered extinctions, resulting in the depletion of lowabundance strains. The depletion of abundant strains is likely due to the killthewinner 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 widelyused 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 $\mathcal{P}(s)$ decays faster or slower than exponentially. A heuristic understanding of how the dynamics of $L$ depend on $\mathcal{P}(s)$ follows from the fact that without general fitness differences, the biases are distributed with characteristic scale $1/\sqrt{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 ${\xi}_{i}$. In the limit of many invasions, the width of the drive distribution becomes comparable to the width of the distribution of the extant $s}_{i$ and cannot decrease further. Thereafter, the shape of $\mathcal{P}(s)$ determines both the width of the bias distribution, and the number of coexisting strains.
If $L\ll 1/{\mathrm{\Sigma}}^{2}$ initially, the $s}_{i$ play little role and the populationweighted mean fitness, $\widehat{s}$, only increases gradually. But once $\widehat{s}$ is a few times $\mathrm{\Sigma}$, the tail of $\mathcal{P}(s)$ will determine the rate of increase of $\widehat{s}$. Henceforth, $\widehat{s}$ will grow steadily with subsequent successful invasions, since strains with ${s}_{i}\ll \widehat{s}$ are very unlikely to persist, and strains with ${s}_{i}\gg \widehat{s}$ are unlikely to have yet occurred. Thus, the range of extant $s}_{i$ will become much narrower than $\mathrm{\Sigma}$. This implies that the distribution over the currently relevant range can be approximated by an exponential distribution $\mathcal{P}(s)\sim {e}^{s/\hat{\mathrm{\Sigma}}}$ with the effective width of the extant $s}_{i$ distribution given by
where the second equality is for the specific models we study (Equation 2). Provided evolution has proceeded long enough that no strains with $s}_{i$ smaller than the original scale $\mathrm{\Sigma}$ survive, $\widehat{\mathrm{\Sigma}}$ will vary slowly for a range of evolutionary time and this sets the scale for variations of $s}_{i$ of both extant strains and of potentiallysuccessful invaders. This suggests that understanding the general behavior at long evolutionary times can be built on understanding the case of the simple exponential distribution ($\psi =1$ for which $\widehat{\mathrm{\Sigma}}=\mathrm{\Sigma}$). The main difference is that now the community size will change as $L\sim 1/{\widehat{\mathrm{\Sigma}}}^{2}$, with $\widehat{\mathrm{\Sigma}}$ changing as $\widehat{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 ${s}_{A}$ comparable to $\widehat{s}$. A simple argument gives an upper bound on how fast $\widehat{s}$ can increase with $T$. In order to get a mean of $\widehat{s}$ in a community of $L$ strains, at least $L$ attempted invaders must have occurred with ${s}_{A}\gtrsim \widehat{s}$. This requires a number of invasion attempts, $T$, such that $T{\int}_{0}^{\hat{s}(T)}\mathcal{P}(s)ds\ge L(T)$. But if $\widehat{s}$ were substantially smaller than this upper bound, many strains would have already occurred with ${s}_{i}\widehat{s}\gg \widehat{\mathrm{\Sigma}}$ and these would persist for a long time, driving $\widehat{s}$ up. We thus make the Ansatz, justified by the analysis and simulation data of ‘Appendix 6’ and Figure 6, that $\widehat{s}(T)$ grows with $\mathrm{log}(T)$ at a rate asymptotically given by this upper bound. For the distributions of interest we then have
where the last implication is due to $L\sim {\widehat{\mathrm{\Sigma}}}^{2}$. These scalings become valid once the distribution of the extant $s}_{i$ is pushed into the tail of $\mathcal{P}(s)$. To crudely take into account the effect of a large initial number of strains when $K\gtrsim 1/{\mathrm{\Sigma}}^{2}$, $\mathrm{log}(T)$ can be replaced by $\mathrm{log}(K+T)$, as for the plots of $L(T)$ and $\widehat{s}(T)$ in Figure 6. Figure 6B shows the theoretical prediction for $\widehat{s}(T)$ in the simplest case of $\psi =1$, and we see that $\widehat{s}(T)/\mathrm{\Sigma}$ is reduced from the $\mathrm{log}[(T+K){\mathrm{\Sigma}}^{2}]$ prediction by only an $\mathcal{O}(1)$ constant, as expected.
At long times, the probability of successful invasion decreases very rapidly with $\widehat{s}$ and, as we show in ‘Appendix 6’, the cumulative number of successful invasions for $\psi \ge 1$ grows very slowly, with $Z\sim {\mathrm{\Sigma}}^{2}{(\mathrm{log}T)}^{32/\psi}$. Nevertheless, for $\psi >1$, the number of strains grows without bound albeit as a sublinear power of the cumulative number of successful invasions. The average number of extinctions per successful invasion gradually decreases towards one as $\mathcal{P}(s)$ in the pertinent range, $s\approx \widehat{s}$, becomes closer and closer to exponential.
For longerthanexponential tails, $\psi <1$, the diversity will decrease (possibly after an initial increase if $K\gg 1/{\mathrm{\Sigma}}^{2}$) and eventually — in practice rather soon — crash as seen in Figure 6A.
Discussion
In this paper, we have answered an important issue of principle: Without any assumption of nichelike 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, $s}_{i$, then as evolution proceeds the average $s}_{i$ of the population gradually increases and pushes into the tail of the $s}_{i$ 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 $s}_{i$ to invade. For broaderthanexponential 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 dediversification. 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 substantialsized 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 modeldependent.
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/\mathcal{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 $\rho $ 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 lowmigration 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 potentiallysuccessful 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 boombust 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 nearestneighbor 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 subpopulations 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 FisherKolmogorovPetrovskyPiscunov (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 longrange 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.
Bacteriaphage interactions and coevolution
An obvious weakness of the LotkaVolterra 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 blockantisymmetricallycorrelated structure of the interaction matrix with the bacteria having no nichestructure (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 LotkaVolterra model studied here: a similar model was further explored in Martis, 2022. Such a bacteriaphage 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 bacteriaphage 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 nichelike structure? Are there models in which a diversifying phase is easier to nucleate? Will the diversification always tend to be limited or stronglyslowed by generalfitness 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, substrain, and subsubstrain 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 $\mathcal{O}(\mathcal{M}L)$, with $\mathcal{M}\equiv \mathrm{log}(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 $\mathcal{O}(L/\mathcal{M})$ strains that happen to be abundant at that time: the variance of this is ${\sum}_{j}{\nu}_{j}^{2}\sim \mathcal{M}/L$, which makes the dynamic fluctuation time $\sim \sqrt{L/\mathcal{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 selforganized 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 $Nm\gg 1$ so that $\mathrm{log}N$ is a few times $\mathcal{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, $\mathcal{M}\sqrt{L}$. 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 $\sim N{\tau}_{\mathit{\text{gen}}}$ with ${\tau}_{\mathit{\text{gen}}}\ll 1$ (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 $\sim \sqrt{L/\mathcal{M}}$. But as long as $N\overline{\nu}m{\tau}_{\mathit{\text{gen}}}\gg 1$, there are many migrants arriving and the dynamics is essentially deterministic. If the islandaverage $\overline{\nu}$ drops enough then the migrations become stochastic with time intervals between them of order $1/(N\overline{\nu}m)$. But when this occurs for substantial time, the global population is likely to be on the way to extinction. We focus on $Nm\gg L$ for which the population dynamics of the persistent strains are not strongly influenced by the stochastic or migratory demographic fluctuations.
Appendix 2
Numerics
Parameters and integration
For all numerics, parameter values unless otherwise mentioned are $I=40$, $m={10}^{5}$, $\gamma =0.8$, $N={10}^{20}$. 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 $\mathcal{O}(\sqrt{L/\mathcal{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 $\alpha $ if ${\nu}_{i,\alpha}<1/N$. In this case ${\nu}_{i,\alpha}$ is set to 0. An irreversible global extinction of strain $i$ occurs if ${\nu}_{i,\alpha}=0$ for all $\alpha $ simultaneously. Recolonization of a locally extinct strain happens when $Nm\overline{\nu}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 speedup by, for a community with population mean $\widehat{s}$, rejecting strains with $s\widehat{s}$ sufficiently negative.
Since the epoch length is $\mathcal{O}(\mathcal{M}L)$ (set by the ecological relaxation time) and the time step scales as $\sqrt{L/\mathcal{M}}$, with each timestep requiring $\mathcal{O}({L}^{2})$ computations to compute the instantaneous growth rates of the strains, the runtime of a single epoch scales as $\mathcal{O}({\mathcal{M}}^{3/2}{L}^{5/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 $\mathrm{E}[{V}_{Mk}{V}_{Pk}]=\mathrm{E}[{V}_{kP}{V}_{kM}]=\rho $ for $M,P\ne k$. In order to preserve the correlation $\gamma $ between acrossdiagonal entries of $V$, we take ${V}_{Mk}=\rho {V}_{Pk}+\sqrt{1{\rho}^{2}}{Z}_{1}$ and ${V}_{kM}=\rho {V}_{kP}+\sqrt{1{\rho}^{2}}(\gamma {Z}_{1}+\sqrt{1{\gamma}^{2}}{Z}_{2})$ where the ${Z}_{i}$ are i.i.d. standard normal random variables. This preserves the desired correlations, with $\mathrm{E}[{V}_{Mk}{V}_{kM}]=\gamma $.
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 $\gamma $ while also having the desired correlations between ${V}_{PP}$, ${V}_{PM}$ and ${V}_{MP}$. Defining $D$, $H$, $J$ and $K$ as
the symmetry conditions we enforce are
In addition, we require that $D=H=J=K$ in the limit $\rho \to 1$, and $\mathrm{E}[HJ]=\gamma $ as $\rho \to 0$. The parameterization that we choose which respects these conditions is
where the ${Z}_{i}$ are i.i.d. standard normal random variables. Now we have $\mathrm{E}[HJ]={\rho}^{2}+\gamma (1{\rho}^{2})$, which does not precisely preserve the desired correlation $\mathrm{E}[HJ]=\gamma $, but it has the correct limit as $\rho \to 0$. Our parameterization gives $\mathrm{E}[DK]={\rho}^{2}(1+\gamma )$, but it is possible to choose an alternative parameterization which has a different value of $\mathrm{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 $\sqrt{L}$. Therefore the choice of the direct parentmutant 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 $\mathcal{O}(\mathcal{M}\sqrt{K})$. Between the introductions of new strains the ecological dynamics are run for epochs of length ${c}_{\mathrm{epoch}}\mathcal{M}L$ with ${c}_{\mathrm{epoch}}=3$. In addition, for most simulations we intersperse regularlength 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 $\sim 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 closetomarginal 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 $12\mathcal{M}L$.
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 $c}_{\text{epoch}$ over one order of magnitude and find similar evolutionary dynamics for all epoch lengths. In the case of $\rho =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 ${c}_{\mathrm{epoch}}=3$, which makes numerics faster, and should not change qualitative conclusions about the slowly diversifying regime.
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 ${\overline{\nu}}_{i\setminus j}$ 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::
where the drive, ${\zeta}_{i}={\sum}_{j}{V}_{ij}{\overline{\nu}}_{j\setminus i}$, can be written in terms of mean field quantities as ${\zeta}_{i}={\sum}_{j}{V}_{ij}{\overline{\nu}}_{j}\gamma X{\overline{\nu}}_{i}$. The term $\gamma X{\overline{\nu}}_{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 $\overline{r}}_{i$ and ${\overline{\nu}}_{i}$ are measurable in our simulation, we need only to subtract off the $\gamma X{\overline{\nu}}_{i}$ term from ${\overline{r}}_{i}$ to get ${\xi}_{i}$. We can find $X$ using the self consistency condition that $X={\sum}_{i}\frac{d{\overline{\nu}}_{i}}{d{\xi}_{i}}={\sum}_{i}\frac{d{\overline{\nu}}_{i}}{d{\zeta}_{i}}.$
First, we make a guess for the susceptibility $X$ which allows a provisional calculation of the ${\zeta}_{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$:
The justification for this form is that it gives ${\overline{\nu}}_{i}\sim {\zeta}_{i}+{s}_{i}c$ for ${\zeta}_{i}+{s}_{i}c>0$. This is the expectation that ${\overline{\nu}}_{i}$ becomes proportional to ${\xi}_{i}$ for positive ${\xi}_{i}$ — though here we include an offset via the parameter $c$ which roughly represents the expected correction to the linear function $\mathcal{N}(\xi )$ at large bias. The quantity ${\zeta}_{i}+{s}_{i}c$ is similar to ${\xi}_{i}$ but not quite the same due to a systematic difference between $c$ and $\overline{\mathrm{{\rm Y}}}$. For large negative bias, the fitted form captures the roughly exponential decrease of $\overline{\nu}$ 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 $\overline{\mathrm{{\rm Y}}}$ within a given epoch, which could be due to the effect of migration which elevates ${\overline{\nu}}_{i}$ for strains with negative ${\xi}_{i}$. The parameter combination $b/a$ is roughly equal to $1/\gamma X$, as expected from the mean field analysis. From the hypothesized form of $\mathcal{N}({\zeta}_{i})$ we can recalculate the susceptibility via $X={\sum}_{i}\frac{\mathrm{\Delta}{\overline{\nu}}_{i}}{\mathrm{\Delta}{\zeta}_{i}}$ (where the derivative is calculated numerically). We then update our value of $X$ by averaging the old guess with the new estimate.
By iterating over the susceptibility until it converges to a value where the assumed drives ${\zeta}_{i}={\sum}_{j}{V}_{ij}{\overline{\nu}}_{j}\gamma X{\overline{\nu}}_{i}$ reproduce the susceptibility via the selfconsistency condition, we can get an estimate for the susceptibility and therefore the drive as well. Then the bias is ${\xi}_{i}\approx {\zeta}_{i}+{s}_{i}\overline{\mathrm{{\rm Y}}},$ where we have subtracted the Lagrange multiplier (neglecting effects of order $\overline{\mathrm{{\rm Y}}}{\overline{\mathrm{{\rm Y}}}}_{\setminus i}$ which are smaller by of order $1/\sqrt{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 ${\sum}_{i}{\overline{\nu}}_{i}^{2}$ and lower limit ${\xi}_{\mathrm{c}}+\overline{\mathrm{{\rm Y}}}$ — 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 $\mathcal{N}(\zeta )$ 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 $\mathcal{N}(\zeta )$.
The function $\mathcal{N}(\xi )$ for the evolved communities with $\rho =0$ looks quite similar to the function for assembled communities However for communities evolved with $\rho =0.95$, $\mathcal{N}(\xi )$ looks substantially different. This indicates that the buildup of correlations in the community requires modifications of the independence assumptions of the mean field theory.
After inferring the function $\mathcal{N}(\zeta )$ and the drive for each strain, we can calculate the fragility, defined as $\mathrm{\Xi}={\sum}_{j}{\chi}_{j}^{2}/(1{\sum}_{j}{\chi}_{j}^{2})$. The closer ${\sum}_{j}{\chi}_{j}^{2}$ is to 1, the more unstable the community is to perturbations. Indeed we find that for simulations diversifying from 50 to 500 strains, ${\sum}_{j}{\chi}_{j}^{2}\cong 0.6$ for the evolved community of $\cong 500$ strains, while early in the simulations, when $L<100$, we find ${\sum}_{j}{\chi}_{j}^{2}\cong 0.75$. When comparing ${\sum}_{j}{\chi}_{j}^{2}$ for assembled and evolved communities, both of 500 strains, there is not a clear difference. Note that the closetomarginal strains, which are more abundant in the assembled community, only contribute small amounts to $\mathrm{\Xi}$, so this is not surprising.
Appendix 3
Dynamics of diversification
Extinctions, diversity crashes, and nucleation from low diversity
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 $\cong 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 ${L}_{0}(K)$: the number of strains that persists after the initial drop in diversity. This somewhatvariable $L}_{0$ 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 $L}_{\mathit{\text{E}}$, 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 $L}_{\mathit{\text{E}}$. This motivates definition of a crossover size ${L}^{*}$, different from ${L}_{0}({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}^{\ast}<{L}_{0}({K}^{\ast})$ since justassembled communities are more susceptible to perturbations than gradually assembled ones.
We can quantify how likely the system is to diversify or crash by inspecting the number of extinctions, $\mathrm{\ell}$, that occur for each invader that successfully enters the ecosystem. The distribution of $\mathrm{\ell}$ 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 $L\cong {L}^{*}$ 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$, $\gamma $ and $\rho $, and is defined as ${K}^{*}$. For some range of $\gamma $ (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 $L\ll {L}^{*}$ 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 10^{5} 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 $L\sim {L}^{*}$ 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}^{*}\approx {L}_{0}({K}^{*}\approx 70)\approx 55$ based on Figures Appendix 3—figure 2 and Appendix 3—figure 1, this estimate suggests $\sim {10}^{60}$ 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 $\gamma $, $\rho $ or $m$. At the marginal point $\gamma =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.
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 Poissonlike distribution of $L$ in the lowdiversity 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$
The central scaling relationship from which our theoretical results follow is that the width of the bias distribution scales as $1/\sqrt{L}$, which scales the same way as $1/\sqrt{\mathcal{L}}\equiv \sqrt{{\sum}_{i}{\overline{\nu}}_{i}^{2}}$. 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 $\mathcal{L}$, which are proportional to each other with $\mathcal{L}\approx 0.74\times L$.
Diversifying phase with $\rho $ close to 1
One of the important results of this work is that there is a continually diversifying phase of ecoevolutionary dynamics both for unrelated invaders and for smalleffect mutations, across the whole range of $\rho \in [0,1)$. In Appendix 3—figure 5, we illustrate the robustness of this phenomenon for highly correlated parents and mutants with $\rho =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 $\rho =0$, $L$ still increases linearly with attempted invasions once the community is nucleated.
Appendix 4
Statistics of $V$ conditioned on evolution
How do the statistics of the interaction matrix ${V}_{ij}$ change when conditioned on evolutionary history? How does this compare to conditioning on nonextinction 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 $\cong 500$ strains, for both $\rho =0$ and $\rho =0.95$. From these, we obtained ensembles of interaction matrices conditioned on gradual evolution with $\rho $ 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 $\mathrm{\Sigma}=0.07$, starting from $K=50$ and reaching steady state with $L\cong 150$ while evolving over the course of $\cong 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 $L}_{0$ similar to the evolved condition.
In Appendix 4—table 1, we list the summary statistics of the interaction matrices that we analyzed: $\widehat{\gamma}$ is the empirical crossdiagonal correlation in the $V$ matrix of extant strains, $\widehat{{\sigma}_{V}}$ is the empirical width of the interaction strength distribution after evolution, ${\widehat{\mu}}_{V}$ is the average of the interactions and $\widehat{\mathrm{\Delta}}$ is the standard deviation, over $i$, of $\frac{1}{L}{\sum}_{j}{V}_{ij}$, which acts like an effective s_{i}, crudely approximating the width of the bias distribution. Aside from $\widehat{\gamma}$, all the quantities are defined including the diagonal terms in $V$.
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\times L$. We thus define quantities $\stackrel{~}{\gamma}$, ${\stackrel{~}{\sigma}}_{V}$, ${\stackrel{~}{\mu}}_{V}$ and $\stackrel{~}{\mathrm{\Delta}}$ to be the empirical matrix statistics measured in units of standard deviations away from their mean in the original $V$ ensemble.
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 ${\mu}_{V}$ for the simulations with $\rho =0.95$; however even this is a relatively small effect.
Our results show that both evolution and assembly tend to make $\widehat{\gamma}$ slightly more negative than the $\gamma $ of the original ensemble, and that these processes also favor — as might have been expected — an increasing mean interaction strength, ${\widehat{\mu}}_{V}$. Of particular interest, given previous discussion, is $\widehat{\mathrm{\Delta}}$, which is similar to the width of the bias distribution. For $L\cong 500$, the statistics of the matrices from the original ensemble have $\widehat{\mathrm{\Delta}}\cong 1/\sqrt{500}\cong 0.0447$. However, once we condition on assembly or evolution, the value of $\widehat{\mathrm{\Delta}}$ can be expected to decrease, because we eliminate those strains with negative bias, meaning that strains $i$ with an anomalously small value of ${\sum}_{j}{V}_{ij}$ will go extinct, making the range of ${\sum}_{j}{V}_{ij}$ narrower. This is indeed observed in Appendix 4—table 3, at least for simulations without general fitness differences. Interestingly we see that $\widehat{\mathrm{\Delta}}$ decreases with evolution more for $\rho =0$ than for $\rho =0.95$, which could be because, in the latter case, the elements in a single row of $V$ are correlated with those of a highabundance parent, which could increase the width of the ${\sum}_{j}{V}_{ij}$ distribution. Nonetheless, we still expect that $\widehat{\mathrm{\Delta}}$ scales as $1/\sqrt{L}$ for evolved communities.
For the evolution with $L$ roughly constant from exponential $\mathcal{P}(s)$, we see that $\widehat{\mathrm{\Delta}}$ after evolution is actually larger than in an unconditioned matrix. However, here the bias has an extra contribution from the general fitness, so $\widehat{\mathrm{\Delta}}$ 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 fixedpoint 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 $\mathrm{C}\mathrm{o}\mathrm{r}\mathrm{r}[{V}_{ik},{V}_{il}]<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, ${\sum}_{j}{V}_{ij}$, will have a smaller width than in the original ensemble. In addition, Bunin finds an $\mathcal{O}(1/K)$ shift in the mean of the interaction matrix toward less competition, which is consistent with ${\hat{\mu}}_{V}>0$. We also see slight but systematic changes in $\widehat{\gamma}$ and $\widehat{{\sigma}_{V}}$, particularly in the case of evolution of correlated mutants with $\rho =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 ${\widehat{\mu}}_{V}$ for $\rho =0.95$. We conjecture that for $\rho =0$ all the effects are small by some power of $L$. But for $\rho $ near one, this is less clear: whether for $L$ large compared to some inverse power of $1\rho $, 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 manybody problem into a singlebody 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 integrodifferential equations for each strain independently. Since the strains are statistically equivalent on each island, we drop the island subscript:
Here, the correlation function $C(t,{t}^{\prime})$ of ${\zeta}_{i}(t)$, the response function $R(t,{t}^{\prime})$, and the islandaveraged abundance ${\overline{\nu}}_{i}$, must be determined selfconsistently. 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 timedependence in ${\zeta}_{i}(t)$ explicit — but elsewhere in the paper when we refer to the timeaveraged or mean drive ${\overline{\zeta}}_{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}^{\prime}$, and quantities will have wellbehaved time averages equal to island averages, denoted by overbars. The average growth rate of strain $i$ on a single island excluding migration is
where the static susceptibility, $X$, is the timeintegrated total response function $X={\int}_{\mathrm{\infty}}^{t}R(t,{t}^{\mathrm{\prime}})d{t}^{\mathrm{\prime}}=\sum _{i}{\chi}_{i}$. Here, we define the individual strain static susceptibility as ${\chi}_{i}\equiv \frac{d{\overline{\nu}}_{i}}{d{\xi}_{i}}$. In steady state with migration, $\u27e8\frac{d\mathrm{log}{\nu}_{i}}{dt}\u27e9=0$; therefore ${\overline{r}}_{i}<0$ must be exactly compensated by the average effect of the migration term $m(\u27e8{\overline{\nu}}_{i}/{\nu}_{i}\u27e91)$, which is always positive (by the CauchySchwarz inequality) and can be much larger than $m$ when ${\nu}_{i}$ is small.
The stochastic DMFT equations cannot be reduced to equations for the correlation and response function and must be treated directly. Their selfconsistent solutions in the STC phase were analyzed in our previous work [PAF], by asymptotics in the large parameter $\mathcal{M}$. This analysis yields superdiffusive random walks of the logabundances, 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 $\overline{\nu}\approx \mathcal{N}(\xi )$ with corrections smaller by factors of $\mathcal{O}(1/\sqrt{L})$.
When the bias of a strain is positive the average input from migration will be relatively small, and hence ${\overline{r}}_{i}\approx 0$ which implies that ${\overline{\nu}}_{i}\approx {\xi}_{i}/(\gamma 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 ${\overline{\nu}}_{i}$ exponentially small in ${\xi}_{i}$ [PAF]. For strains with more negative bias, ${\xi}_{i}<{\xi}_{\mathrm{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 $\mathcal{O}(1/\sqrt{L})$, with the mean drives, ${\zeta}_{i}$, gaussian distributed with average zero and standard deviation of order $1/\sqrt{L}$. More precisely, since ${\xi}_{i}$ is determined by the community in the absence of $i$, the ${V}_{ij}$ that determine it are independent of the abundances in the community. The variance of $\zeta $ is $1/\mathcal{L}\equiv {\sum}_{j}{\overline{\nu}}_{j}^{2}$ (with only small corrections from the neglect of the strain itself): this simple result is explicit in the selfconsistency condition of the DMFT correlation function. We therefore use $\mathcal{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 ${\xi}_{\mathrm{c}}$, which is of order $1/\sqrt{\mathcal{L}}$ — the magnitude of ${\xi}_{\mathrm{c}}$ and $\mathcal{L}$ must be determined selfconsistently. 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 $\mathcal{P}(s)$ in perfectly antisymmetric model
The antisymmetric model with $\gamma =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 $\{{\nu}_{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 temperaturelike quantity $\mathrm{\Theta}$ which is roughly the range of the fluctuations in $\mathrm{log}{\nu}_{i}$, with the average abundances ${\overline{\nu}}_{i}={\nu}_{i}^{*}$.
Here we analyze the effect of introducing independent exponentially distributed $s}_{i$ 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 ${\xi}_{\mathrm{c}}=0$ and the mean $s$ is simply $\widehat{s}=\overline{\mathrm{{\rm Y}}}$. In the limit $K{\mathrm{\Sigma}}^{2}\gg 1$, the distribution of the drives with scale $1/\sqrt{K}$ will be much narrower than that of the $s}_{i$. Let us define ${x}_{i}\equiv {\zeta}_{i}+{s}_{i}$ for convenience. The scale of $p(x)$, the x_{i} distribution, will be $\mathrm{\Sigma}$ 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 ${\sigma}_{\zeta}^{2}$:
With the Lagrange multiplier $\overline{\mathrm{{\rm Y}}}$, the average abundances ${\overline{\nu}}_{i}$ at the mean field fixed point will be ${\overline{\nu}}_{i}=\mathrm{max}(0,\frac{{x}_{i}\overline{\mathrm{{\rm Y}}}}{X})=\mathrm{max}(0,{\xi}_{i}/X)$ with $X$ the total static susceptibility given by ${\sum}_{i}{\chi}_{i}$. Let us define $\varphi $ as the fraction of strains that persists with positive bias. When $\overline{\mathrm{{\rm Y}}}>{\sigma}_{\zeta}[{\sigma}_{\zeta}/\mathrm{\Sigma}+\mathcal{O}(1)]$, we can make the approximation
Self consistency requires that the average abundance over the initial $K$ strains is
Now we combine this relation with the self consistency condition
to obtain $\varphi \approx \frac{1}{K{\mathrm{\Sigma}}^{2}}$. Therefore, the number of persistent strains, when $K{\mathrm{\Sigma}}^{2}\gg 1$ is ${\mathrm{\Sigma}}^{2}$, with coefficient one.
To obtain the behavior for arbitrary $K$, one must solve the mean field self consistency equations, with the exact distribution
where $\mathrm{\Phi}$ is the standard normal cumulative distribution function. There are three mean field equations for the three unknowns $\overline{\mathrm{{\rm Y}}}$, ${\sigma}_{\zeta}$ and $X$.
The number of persistent strains is $K{\int}_{\overline{\mathrm{{\rm Y}}}}^{\mathrm{\infty}}p(x)dx$, where $\overline{\mathrm{{\rm Y}}}$ solves the above equations. Since for perfectly antisymmetric interactions, $\widehat{s}=\overline{\mathrm{{\rm Y}}}$, the solution of these allows us to check that $\widehat{s}/\mathrm{\Sigma}$ increases with $\mathrm{\Sigma}$, 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{\mathrm{\Sigma}}^{2}\gg 1$, we have $\overline{\mathrm{{\rm Y}}}\approx \mathrm{\Sigma}(1+\mathrm{log}(K{\mathrm{\Sigma}}^{2}))$, $X\approx 1/\mathrm{\Sigma}$ and ${\sigma}_{\zeta}\approx \sqrt{2}\mathrm{\Sigma}$, which is consistent with the asymptotic Ansatz for $p(x)$ for the present strains: since $\overline{\mathrm{{\rm Y}}}{\sigma}_{\zeta}^{2}/\mathrm{\Sigma}\gg {\sigma}_{\zeta}$ all persistent strains have large enough $x$ for the approximation to be valid. Then the distribution of $\zeta $ for the persistent strains is
valid except for $\zeta >\overline{\mathrm{{\rm Y}}}=\mathrm{\Sigma}(1+\mathrm{log}(K{\mathrm{\Sigma}}^{2}))$ which is a negligible fraction of the distribution. The width of both the $\zeta $ and the $s$ distributions for the persistent strains is of order $\mathrm{\Sigma}$. As expected more generally, these are comparable once $\widehat{s}$ is in the tail of $\mathcal{P}(s)$. This calculation could be done for any choice of $\mathcal{P}(s)$, and should give the behavior of $\overline{\mathrm{{\rm Y}}}$ with $K$, which is the same as the behavior of $\widehat{s}$ with $T$ in the evolving phase of the STC: namely $\widehat{s}\sim \mathrm{\Sigma}{(\mathrm{log}T)}^{1/\psi}$.
Joint distribution of $\zeta $ 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 $\xi}_{A}>{\xi}_{\mathrm{c}$. Since the invader’s bias is a sum of ${s}_{A}$, ${\zeta}_{A}$ and $\overline{\mathrm{{\rm Y}}}$, we can calculate the distributions of ${s}_{A}$ and ${\zeta}_{A}$ conditional on invasion.
We analyze the case with an exponential distribution of the $s}_{i$ with scale $\mathrm{\Sigma}$. Since the distribution of the invader ${\zeta}_{A}$ and ${s}_{A}$ are independent, we can write their joint distribution as
Conditioning on successful invasion multiplies by a factor of $\mathrm{\Theta}({\zeta}_{A}+{s}_{A}{\xi}_{\mathrm{c}}\overline{\mathrm{{\rm Y}}})$, which enforces that $\xi}_{A}>{\xi}_{\mathrm{c}$. We are interested in the large time limit when $\overline{\mathrm{{\rm Y}}}=\hat{s}+\mathcal{O}(\mathrm{\Sigma})\gg \mathrm{\Sigma}$. In this limit, the restriction that ${s}_{A}>0$ has negligible effect and the analysis is simple.
The invasion probability is
which decays exponentially as ${e}^{\widehat{s}/\mathrm{\Sigma}}$, as expected. The marginal distribution of ${s}_{A}$ conditioned on invasion is
where $\mathrm{\Phi}$ is the standard normal cdf. Similarly, the distribution of ${\zeta}_{A}$ conditioned on invasion is
and the bias distribution conditioned on invasion is
Of particular interest is the mean ${s}_{A}$:
since $\mathrm{E}({\zeta}_{A}\mathrm{invasion})\approx 1/(\mathrm{\Sigma}\mathcal{L})$ and $\mathrm{E}({\xi}_{A}\mathrm{invasion})\approx {\xi}_{\mathrm{c}}+\mathrm{\Sigma}$.
Effects of shape of tail of general fitness distribution
Here, we analyze the evolutionary dynamics for the class of generalfitness distributions with high$s$ tails parametrized by $\psi $, 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, $\widehat{s}$ increases past $\mathrm{\Sigma}$ and pushes up into the tail of $\mathcal{P}(s)$. The rate of successful invasions then decreases rapidly since invaders need general fitness (at least) comparable to $\widehat{s}$ in order to invade: the probability of successful invasion is of order ${p}_{\text{invade}}\sim \hat{\mathrm{\Sigma}}\mathcal{P}(\hat{s})$ with the prefactor from the width of the distribution of ${s}_{A}$ near $\widehat{s}$. In this regime, the analysis of ‘Appendix 6’ can be carried over with $\mathrm{\Sigma}$ replaced by $\widehat{\mathrm{\Sigma}}$. The drive on the invader from the extant community, ${\zeta}_{A}\sim 1/\sqrt{L}$, will be of order $\widehat{\mathrm{\Sigma}}$, yielding the key scaling prediction $L\sim {\widehat{\mathrm{\Sigma}}}^{2}$ which is tested in Appendix 6—figure 1A.
Successful invaders will have ${s}_{A}=\hat{s}+\mathcal{O}(\hat{\mathrm{\Sigma}})$, since $\widehat{\mathrm{\Sigma}}$ is the typical amount by which the invader’s general fitness is likely to be larger (or smaller) than the $\widehat{s}$ of the extant community. Therefore, we expect that the change in the communitymean $\widehat{s}$ per successful invasion will be $d\widehat{s}/dZ\sim \widehat{\mathrm{\Sigma}}/L$ (obtainable from ‘Appendix 6’): this scaling is tested in Appendix 6—figure 1B. Substituting for $L\sim {\widehat{\mathrm{\Sigma}}}^{2}$, we have $d\widehat{s}/dZ\sim {\widehat{\mathrm{\Sigma}}}^{3}$. Integrating this equation and using $\widehat{\mathrm{\Sigma}}={\mathrm{\Sigma}}^{\psi}/{\widehat{s}}^{\psi 1}$ one obtains scaling laws quoted in the main text, valid for $Z\gg 1/{\mathrm{\Sigma}}^{2}$ and $L\sim {Z}^{\frac{2\psi 2}{3\psi 2}}$ with the latter exponent less than one for the relevant range of $\psi $.
For $\psi >1$, both $L$ and $\widehat{s}$ increase sublinearly 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 $\psi <1$, $L$ decreases with $Z$ and $\widehat{s}$ increases superlinearly 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{\widehat{\mathrm{\Sigma}}}^{2}$ is shown versus $Z$ in Appendix 6—figure 1A. this combination is predicted to approach a $\psi $independent constant for large $Z$. The value of this constant is consistent with the 0.7 found for $\psi =1$ in Figure 3A. The other key scaling law, $d\widehat{s}/dZ\sim \widehat{\mathrm{\Sigma}}/L$, is tested in Appendix 6—figure 1B, with the predicted coefficient of $\cong 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 latetime asymptotic regime.
To get the growth of $\widehat{s}$ — and hence the other quantities — with the number of attempted invasions, we write; $d\widehat{s}/dT\sim {p}_{\text{invade}}\widehat{\mathrm{\Sigma}}/L$ and integrate this using the scaling relations to get
with the unknown $\mathcal{O}(1)$ correction reflecting an additive uncertainty of order $\widehat{\mathrm{\Sigma}}$ in $\widehat{s}$. Up to this correction, the expression has an identical form to the upper limit from ${p}_{\text{invade}}T\ge L$ as explained in the main text, and quoted there without the $\mathrm{log}\mathrm{log}$ correction (which vanishes for $\psi =1$). This prediction, with $T$ replaced by $T+K$ (without the $\mathrm{log}\mathrm{log}$ part or the $\mathcal{O}(1)$ correction), is plotted in Figure 6B. Slow crossovers from the initial nonuniversal behavior towards the predicted behaviors are observable with the asymptotic predictions thus not welltestable in the timedependences.
Appendix 7
Correlated general fitness mutations
We have shown that including general fitnesses with $\psi \ge 1$ causes evolution with unrelated invaders to slow down, with $Z$ increasing only as power of $\mathrm{log}T$. 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 $\mathcal{P}(s)$, we can take the mutant fitness to be sampled from a Markov chain starting at the parent, with stationary measure given by $\mathcal{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 $\mathcal{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 $\mathrm{\Sigma}$. The variance of the gaussian determines the degree to which the parent and mutant are correlated, parametrized by ${\rho}_{s}$, which need not be the same $\rho $ defined by the correlations between parent and mutant interactions. The crucial difference from uncorrelated $s$ is that the probability of $s\gtrsim \widehat{s}$ is now independent of $\widehat{s}$ rather than decreasing exponentially with $\widehat{s}/\mathrm{\Sigma}$. This enables $\widehat{s}$ and $Z$ to increase linearly in $T$. However, $L$ still saturates to a value of order $1/{\mathrm{\Sigma}}^{2}$, with a coefficient that depends on the correlations. Appendix 7—figure 1 shows simulations that confirm these predictions.
Nonexponential $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 $\psi >1$ case of particular interest, as $\widehat{s}$ grows the mean $s\widehat{s}$ of mutants becomes more and more negative and overwhelms the random part of ${s}_{M}{s}_{P}$. There is a then a scale of $s$ that diverges as ${\rho}_{s}\to 1$, such that for $\widehat{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 ${\zeta}_{0}={\sum}_{j}{V}_{0j}{\overline{\nu}}_{j}$, while after the invasion it changes to
where ${\overline{\nu}}_{A}$ is the mean abundance at which strain $A$ establishes and ${\overline{\nu}}_{j\setminus A}$ is the average abundance of strain $j$ before $A$ invades. The sum on $j$ does not include 0 or $A$ and the $\delta $’s denote changes that result from the invasion. Since ${\overline{\nu}}_{A}$ is of order $1/L$, its direct effect on the drive is smaller than ${\zeta}_{0}$ by a factor of $1/\sqrt{L}$. However $A$ will also change all of the other biases by similarmagnitude random amounts with uncorrelated signs, causing $\delta {\overline{\nu}}_{j}\approx {\chi}_{j}\delta {\zeta}_{j}$. Thus $\mathrm{E}[\delta {\zeta}_{0}^{2}]={\overline{\nu}}_{A}^{2}+{\sum}_{j}{V}_{0j}^{2}{\chi}_{j}^{2}\delta {\zeta}_{j}^{2}$ with the expectation over the interactions with the probe strain. Since the properties of the extant strains depend neither on the ${V}_{0j}$ nor the perturbations via ${V}_{jA}$, the averages over each of the squared factors can be performed separately, leaving ${\sum}_{j}{\chi}_{j}^{2}$ times the average over $j$ of $\delta {\zeta}_{j}^{2}$. 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 ${\mathrm{E}}_{\text{strains}}[\delta {\zeta}_{0}^{2}]={\overline{\nu}}_{A}^{2}/(1{\sum}_{j}{\chi}_{j}^{2})$. 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 ${\sum}_{j}{(\delta {\overline{\nu}}_{j})}^{2}=\mathrm{\Xi}{\overline{\nu}}_{A}^{2}$ in terms of the nonlinear response, which we call the fragility: $\mathrm{\Xi}={\sum}_{j}{\chi}_{j}^{2}/(1{\sum}_{j}{\chi}_{j}^{2})$. This quantity is analogous to the spinglass 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 nichephase of the Lotka Volterra model with large negative diagonal interactions of magnitude $Q=\mathrm{E}[{V}_{ii}]$, the fragility diverges as $K$ increases to the stability boundary at ${K}_{\text{crit}}\sim {Q}^{2}$, 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 ${\zeta}_{j}$ induced by a successful invader, there is also a systematic change, as seen in Figure 4. This is because $\delta {\zeta}_{0}$ and ${\zeta}_{0}$ involve the same set of random ${V}_{0j}$ and are thus correlated. We can use the fact that $\mathrm{E}[XY=y]=\frac{\mathrm{Cov}(X,Y)}{\mathrm{Var}(Y)}y$ for zero mean gaussian random variables $X$ and $Y$ to see that the conditional expectation is $\mathrm{E}[\delta {\zeta}_{0}{\zeta}_{0}]={\zeta}_{0}{\sum}_{j}{\overline{\nu}}_{j}\delta {\overline{\nu}}_{j}/{\sum}_{j}{\overline{\nu}}_{j}^{2}$. Since each ${\overline{\nu}}_{j}$ and $\delta {\overline{\nu}}_{j}$ are correlated as a consequence of ${\zeta}_{j}$ and $\delta {\zeta}_{j}$ being correlated, we again have to selfconsistently determine this correlation. Unfortunately this is more complicated as the $\delta {\overline{\nu}}_{j}$ also have a contribution from small changes in the function $\mathcal{N}(\zeta )$ caused by the invaded strain. However one can understand the form of the conditional expectations by noting that after the new strain has invaded,
This implies that the last term must be negative if $L$ and $\mathcal{L}$ are increasing or staying constant. Therefore the correlation between ${\zeta}_{j}$ and $\delta {\zeta}_{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, $\mathrm{\Xi}$. If the fragility is high, one can show that $2{\sum}_{j}{\overline{\nu}}_{j}\delta {\overline{\nu}}_{j}$ is proportional to $\mathrm{\Xi}\mathrm{E}[{\nu}_{A}^{2}]$, since the last two terms in Equation A8.2 should be of similar magnitude to preserve $\mathcal{L}$: in this limit, as argued below, the diversity will decrease.
The above Markovian analysis is an approximation because $\delta {\zeta}_{j}$ for each extant strain is correlated with its whole past trajectory — as it involves many of the same ${V}_{0j}$ — and its conditional expectation will be a weighted sum over the full past ${\zeta}_{j}({T}^{\prime})$. Furthermore, it is possible that conditioning on all this past could suppress $\mathrm{E}[\delta {\zeta}_{j}^{2}]$ 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 $\sim 1/\sqrt{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/\sqrt{\mathcal{L}}$, the effective community size is $\mathcal{L}\approx CL$ with $C<1$. We thence obtain a Fokker Planck equation for the number density $N$ of the drives:
where we impose the absorbing boundary condition at the critical bias corresponding to $N({\xi}_{\mathrm{c}}+\overline{\mathrm{{\rm Y}}},T)=0$ with ${\xi}_{\mathrm{c}}$ and $\overline{\mathrm{{\rm Y}}}$ both scaling as $1/\sqrt{L}$. The last term in Equation A8.3 is the truncatedgaussian drive distribution of successfully invading strains with scale $1/\sqrt{\mathcal{L}}$. The fraction of invasions that are successful, $dZ/dT$, is just the integral of the truncated gaussian, $1\mathrm{\Phi}[\sqrt{\mathcal{L}}({\xi}_{\mathrm{c}}+\overline{\mathrm{{\rm Y}}})]$. The number of strains, $L(T)$, is changing and at this point unknown: it is determined selfconsistently from $L(T)=\int d\zeta N(\zeta ,T)$. It is straightforward to show, as we do in ‘Appendix 8’, that this FokkerPlanck equation admits a scaling solution with the Ansatz that all the quantities scale as $1/\sqrt{L}$. The rate of diversification, $U=dL/dT$, is of order one and determined by an eigenvaluelike condition: it can be either positive or negative depending on ${\xi}_{\mathrm{c}}+\overline{\mathrm{{\rm Y}}}$ and the coefficients, $B,D,C$. If the community is very fragile ($\mathrm{\Xi}$ 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 $({\xi}_{\mathrm{c}}+\overline{\mathrm{{\rm Y}}})\sqrt{\mathcal{L}}$ is large, few strains will go extinct and the community will diversify. Of course, these quantities are determined selfconsistently, depending on both the distribution of biases and the function $\mathcal{N}(\xi )$, 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, $\zeta ({T}^{\prime})$, over all ${T}^{\mathrm{\prime}}<T$. If one makes the Ansatz of a scaling solution with $L$ steadily increasing, then for large $T$, $L\sim T$ and in the variables scaled by $T$ the weighting function should be a function of of $q=\mathrm{log}(T/{T}^{\prime})$, 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 integrodifferential generalization of the FokkerPlanck 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.
FokkerPlanck 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(\zeta ,T)$ whose integral at any time gives us the number of strains. The crucial Ansatz is that all the quantities ($\zeta $, ${\xi}_{\mathrm{c}}$ and $\overline{\mathrm{{\rm Y}}}$) scale as $1/\sqrt{L}$, and we now take $L$ to be a function of $T$. Then we can hypothesize a scaling solution $N(\zeta ,T)\sim {L}^{3/2}g(\zeta \sqrt{L})$ for some function $g$ so that $\int N(\zeta ,T)d\zeta \sim 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\equiv \zeta \sqrt{L}$, and derivatives with respect to $u$ by primes: then Equation A8.3 becomes
where we define ${u}_{0}\equiv \sqrt{L}({\xi}_{\mathrm{c}}+\overline{\mathrm{{\rm Y}}})$. We can make an Ansatz of the form
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
We can solve this equation numerically with boundary conditions $g({u}_{0})=0$ and $g(\mathrm{\infty})=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 $\zeta}_{i}>{\xi}_{\mathrm{c}}+\overline{\mathrm{{\rm Y}}$ for every extant strain, and second comes from the need to have $L<\mathrm{\infty}$. The solution depends on the parameters u_{0}, $B$, $D$ and $U$. By enforcing ${\int}_{{u}_{0}}^{\mathrm{\infty}}g(u)du=1$, we can find the value of $U$ in terms of ${u}_{0},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 unitvariance 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/dT\cong 0.24$ and $B=2.11$, $D=1.79$ adjusted to enforce normalization. Note that $C=\mathcal{L}/L\cong 0.74$. The mode of the truncated gaussian is $\cong \overline{\mathrm{{\rm Y}}}\sqrt{\mathcal{L}}$ as expected from Figure 4. The bias distributions for assembled and evolved communities are very similar except near ${\xi}_{\mathrm{c}}$, which is what one expects from the Markov approximation. For ${\zeta}_{c}={\xi}_{\mathrm{c}}\mathrm{{\rm Y}}\gg 1/\sqrt{\mathcal{L}}$, when the ratio of extant to invading widths, $\sqrt{DC/B}\cong 0.79$, is close to unity, the bulk of the drive distribution is close to gaussian with mean zero, which is reflected in Figure 5B.
FokkerPlanck equation with independent general fitness differences
If incoming strains have independently drawn general fitnesses, $s}_{i$, then we can write an evolution equation for the joint distribution of the drives and the $s}_{i$:
where we abbreviate $\mathcal{P}(ss>{\xi}_{\mathrm{c}}+\overline{\mathrm{{\rm Y}}}\zeta )$ by ${\mathcal{P}}_{>}(s)$. However, unlike the case with general fitness differences, now $\overline{\mathrm{{\rm Y}}}$ does not scale as $1/\sqrt{L}$. Indeed, the scaling of $\overline{\mathrm{{\rm Y}}}$ with $Z$ determines how $L$ depends on $Z$. If we specify the form of $\frac{d\overline{\mathrm{{\rm Y}}}}{dZ}$, then we can look for solutions which travel up in the $s$ direction at the same speed as $\overline{\mathrm{{\rm Y}}}$.
We can define similarity variables $u\equiv \zeta \sqrt{L}$ and $v\equiv \sqrt{L}(s\overline{\mathrm{{\rm Y}}})$. Then we look for solutions of the form $N(\zeta ,s,Z)\sim {L}^{2}g(u,v)$ so that the integral of $N$ over $\zeta $ and $s$ is proportional to $L$. The boundary condition is that $g$ vanishes along the line $u+v={u}_{\text{crit}}$. Plugging this into our Fokker Planck equation, and defining ${u}_{\text{crit}}\equiv \sqrt{L}{\xi}_{\mathrm{c}}$, we obtain
where subscripts denote derivatives of $g$. We can see that the ${\mathcal{P}}_{>}(s)$ distribution generically breaks the scaling — because it depends on $s=v/\sqrt{L}$, not only on $v$. But we can look for solutions which obey the scaling with $1/\sqrt{L}$. Take $v/\sqrt{L}$ as exponentially distributed with scale $\mathrm{\Sigma}$. Then in order for our scaling Ansatz to be valid, the quantity $\frac{v}{\mathrm{\Sigma}\sqrt{L}}$ must be independent of $L$, so we must have $\sqrt{L}\sim 1/\mathrm{\Sigma}\phantom{\rule{thickmathspace}{0ex}}\u27f9\phantom{\rule{thickmathspace}{0ex}}L\sim {\mathrm{\Sigma}}^{2}$. This requires $\frac{d\overline{\mathrm{{\rm Y}}}}{dZ}\sim {\mathrm{\Sigma}}^{3}$ in order to make the left hand side of Equation A8.8 independent of $L$. Furthermore, we have $dL/dZ\to 0$ since $L\sim {\mathrm{\Sigma}}^{2}$. Therefore, we have
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 s_{i}.
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.
LotkaVolterra 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 $\rho $ 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 $\mathrm{log}\nu $ to obtain equations for parent and mutant frequencies ${\nu}_{P}$ and ${\nu}_{M}$:
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 parentmutant replacement in the evolving STC.
In the absence of the mutant, ${\overline{\nu}}_{P}\approx {\xi}_{P}/(\gamma X)$ (which is the correct behavior for large positive ${\xi}_{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 $\xi}_{M}\rho {\xi}_{P}>{\xi}_{\mathrm{c}$ and the condition for invasion of the parent into the mutant is likewise $\xi}_{P}\rho {\xi}_{M}>{\xi}_{\mathrm{c}$. In terms of the drives, these conditions are $\zeta}_{M}\rho {\zeta}_{P}>(1\rho )\overline{\mathrm{{\rm Y}}}+{\xi}_{\mathrm{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 ${\zeta}_{M}=\rho {\zeta}_{P}+\sqrt{1{\rho}^{2}}{Z}_{1}$ with $s}_{i$ a gaussian random variable with mean 0 and scale ${\zeta}_{P}$, the probabilities of coexistence between parent and mutant, and of replacement of the parent, are
where $\mathrm{\Phi}$ is the standard normal cdf, and $\stackrel{~}{\mathrm{{\rm Y}}}$ and ${\stackrel{~}{\xi}}_{c}$ are the Lagrange multiplier and critical bias, both normalized by the scale of the $\zeta $ distribution. Note that the expressions are only valid for $\stackrel{~}{\mathrm{{\rm Y}}}<1\frac{{\stackrel{~}{\xi}}_{c}}{1{\rho}^{2}}$: otherwise, in this approximation, the coexistence probability vanishes and all invasions are replacements. The probability of mutant invasion is ${p}_{\text{invade}}=1\mathrm{\Phi}\left[\frac{\stackrel{~}{\mathrm{{\rm Y}}}(1\rho )+{\stackrel{~}{\xi}}_{c}}{\sqrt{1{\rho}^{2}}}\right]$. For $1\rho \ll 1$ and ${\stackrel{~}{\xi}}_{c}=0$, the invasion probability is ${p}_{\text{invade}}\approx \frac{1}{2}(1\stackrel{~}{\mathrm{{\rm Y}}}\sqrt{(1\rho )/\pi})$ and the coexistence probability predicted by this simple approximation is ${p}_{\text{coexist}}\approx (1\stackrel{~}{\mathrm{{\rm Y}}})\sqrt{(1\rho )/\pi}$. With $\rho =0.99$ and $\stackrel{~}{\mathrm{{\rm Y}}}=0$ (in the STC $\stackrel{~}{\mathrm{{\rm Y}}}$ 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 ${p}_{\text{coexist}}/{p}_{\text{invade}}\cong 0.12$ which is much smaller than the observed ${p}_{\text{coexist}\text{invade}}>0.5$ (Appendix 9—figure 1).
The main problem with this approximation is that it needs to be used when the biases are close to the critical ${\xi}_{\mathrm{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 $\rho $. 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 boombust dynamics with migration makes coexistence much easier, but understanding this even semiquantitatively 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 $\rho $ is extremely close to unity, the coexistence probability of parent and mutant will go to zero as $\sqrt{1\rho}$, as in the simple approximation from ‘Appendix 9’, but with a small coefficient that might be of order some inverse power of $\mathcal{M}$. In addition one should note that the natural timescale over which differences between the parent and mutant will accumulate is of order ${(1{\rho}^{2})}^{1/2}\mathcal{M}\sqrt{L}$ with the coefficient only $\cong 20$ for $\rho =0.999$. Therefore even for $\rho $ this close to 1, the typical differences between parent and mutant will still be felt on timescales of order $\mathcal{M}L$ over which the closetomarginal 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 $\rho $ 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, $\rho >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 $\rho $ 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{\rho}^{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 $L}_{0$, the correlations between the biases of the extant strains can not be ignored and the approximation of independent ${\zeta}_{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.
References

Molecular signatures of resource competition: Clonal interference favors ecological diversification and can lead to incipient SpeciationEvolution; International Journal of Organic Evolution 75:2641–2657.https://doi.org/10.1111/evo.14315

Maintenance of plant species diversity by pathogensAnnual Review of Ecology, Evolution, and Systematics 46:305–325.https://doi.org/10.1146/annurevecolsys112414054306

Prochlorococcus: the structure and function of collective diversityNature Reviews. Microbiology 13:13–27.https://doi.org/10.1038/nrmicro3378

Ecological communities with LotkaVolterra dynamicsPhysical Review. E 95:042414.https://doi.org/10.1103/PhysRevE.95.042414

Largedimensional replicator equations with antisymmetric random interactionsJournal of the Physical Society of Japan 71:429–431.https://doi.org/10.1143/JPSJ.71.429

Macarthur’s consumerresource modelTheoretical Population Biology 37:26–38.https://doi.org/10.1016/00405809(90)90025Q

General theory of competitive coexistence in spatiallyvarying environmentsTheoretical Population Biology 58:211–237.https://doi.org/10.1006/tpbi.2000.1486

Replicators with random interactions: A solvable modelPhysical Review. A, General Physics 39:4333–4336.https://doi.org/10.1103/physreva.39.4333

The wave of advance of advantageous genesAnnals of Eugenics 7:355–369.https://doi.org/10.1111/j.14691809.1937.tb02153.x

Nonequilibrium dynamics of spin glassesPhysical Review B 38:373–385.https://doi.org/10.1103/PhysRevB.38.373

Community structure follows simple assembly rules in microbial microcosmsNature Ecology & Evolution 1:1–7.https://doi.org/10.1038/s415590170109

Random replicators with asymmetric couplingsJournal of Physics A 39:3853–3869.https://doi.org/10.1088/03054470/39/15/001

On the volterra and other nonlinear models of interacting populationsReviews of Modern Physics 43:231–276.https://doi.org/10.1103/RevModPhys.43.231

Longterm stability and red queenlike strain dynamics in marine virusesNature Microbiology 5:265–271.https://doi.org/10.1038/s415640190628x

Genetic demixing and evolution in linear stepping stone modelsReviews of Modern Physics 82:1691–1718.https://doi.org/10.1103/RevModPhys.82.1691

Phase transition and 1/f noise in a game dynamical modelPhysical Review Letters 69:1616–1619.https://doi.org/10.1103/PhysRevLett.69.1616

Metabolic tradeoffs promote diversity in a model ecosystemPhysical Review Letters 118:028103.https://doi.org/10.1103/PhysRevLett.118.028103

Solvable model of a complex ecosystem with randomly interacting speciesJournal of Physics A 22:3447–3460.https://doi.org/10.1088/03054470/22/17/011

Subdiffusive fluctuations of "pulled" fronts with multiplicative noisePhysical Review. E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 62:R13–R16.https://doi.org/10.1103/physreve.62.r13

Diffusion coefficient of propagating fronts with multiplicative noisePhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 65:012102.https://doi.org/10.1103/PhysRevE.65.012102

Explaining microbial population genomics through phage predationNature Reviews. Microbiology 7:828–836.https://doi.org/10.1038/nrmicro2235

Complex interactions can create persistent fluctuations in highdiversity ecosystemsPLOS Computational Biology 16:e1007827.https://doi.org/10.1371/journal.pcbi.1007827

Tractable models of ecological assemblyEcology Letters 24:1029–1037.https://doi.org/10.1111/ele.13702

Innovation rather than improvement: A solvable highdimensional model highlights the limitations of scalar fitnessJournal of Statistical Physics 172:74–104.https://doi.org/10.1007/s1095501819566

Species abundance patterns in complex evolutionary dynamicsPhysical Review Letters 93:178102.https://doi.org/10.1103/PhysRevLett.93.178102

Coevolution maintains diversity in the stochastic "kill the winner" modelPhysical Review Letters 119:268101.https://doi.org/10.1103/PhysRevLett.119.268101

Rank abundance relations in evolutionary dynamics of random replicatorsPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 78:031924.https://doi.org/10.1103/PhysRevE.78.031924
Article and author information
Author details
Funding
National Science Foundation (PHY160760)
 Aditya Mahadevan
National Institutes of Health (R01AI13699201)
 Aditya Mahadevan
Simons Foundation (Sabbatical Fellowship)
 Daniel S Fisher
National Science Foundation (PHY2210386)
 Aditya Mahadevan
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Pankaj Mehta for useful discussions and comments on the manuscript. This work was supported in part by the National Science Foundation via PHY1607606,and PHY2210386, 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
 Preprint posted: May 28, 2022 (view preprint)
 Received: August 16, 2022
 Accepted: April 27, 2023
 Accepted Manuscript published: April 28, 2023 (version 1)
 Version of Record published: June 28, 2023 (version 2)
Copyright
© 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.
Metrics

 1,002
 views

 221
 downloads

 7
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
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)
Further reading

 Ecology
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 contextdependent Joint Species Distribution Models to a systematic cameratrap 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 farranging ecological consequences.

 Ecology
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.