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 small-effect 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 ${\displaystyle \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 ${\displaystyle s_i}$ for strain $i$) can either slow down, prevent, or reverse diversification. For distributions of the ${\displaystyle 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 $\log (T)$. If the distribution of the ${\displaystyle s_i}$ has a broader-than-exponential tail, the diversity decreases and crashes.

The key property of a strain, in terms of which one can understand its behavior, is its bias: defined as the rate at which its population would change when at low-abundance and without migration (Figure 1B). The crucial effect of migration in the STC is to stabilize many strains with negative bias which would have gone extinct without migration. Only if its bias is strongly negative will a strain go globally extinct. The bias of strain $i$, denoted by ${\xi }_i$, has an intrinsic contribution from its general fitness ${\displaystyle 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 ${\displaystyle 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 all-at-once from unrelated strains. At early stages of the evolution, most of the close-to-marginal, low-abundance strains are pushed out by the perturbations caused by the invading strains. This extinction process causes the shape of the abundance distributions of assembled and evolved communities to differ at low abundances. Later, in the steadily diversifying state, the numbers of extinctions caused by each invader has a roughly exponential distribution, which is consistent with our theoretical expectations. In contrast to the qualitative (albeit modest) changes in abundance distributions, we find that evolution has only a small effect on the statistics of the interactions between strains.

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

We 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 Lotka-Volterra model (Goel et al., 1971), with each strain $i$ having an intrinsic growth rate which is modulated by its interactions with all the other strains. These interactions are conveniently represented in a matrix $W$ where $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 ${\displaystyle s_i}$) by $V_{ij}$. Since the $V_{ij}$ and ${\displaystyle 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 (${\displaystyle \gamma >0}$) or host-pathogen-like (${\displaystyle \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 host-pathogen models with the appropriate block sub-matrix structure, as discussed further in ‘Bacteria-phage interactions and coevolution’.

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 time-averaged abundance of strain $i$ on a single island; this is a crucial self-consistency condition. The magnitude, $m$ of the migration rate, is also of fundamental importance. If $m$ is too small the migration is too rare to repopulate islands after local extinctions. If $m$ is too large and the local dynamics is chaotic, the chaos will synchronize across the islands and the total population of each strain will fluctuate wildly, rapidly driving most strains extinct. We will focus on the wide intermediate $m$ regime, which spans several orders of magnitude when $K$ is large [PAF].

With large populations on each island, demographic fluctuations have little effect on the dynamics. Even when the local population of a strain is small, if it has positive growth rate, fluctuations will not matter much, while if it has negative growth rate it will go deterministically extinct which occurs when the fractional abundance drops below the extinction threshold of $1/N$ indicated in Figure 1A by the horizontal purple line. The value of the extinction threshold does not much affect the behavior as long as it is much below the lower limit of the abundance caused by migration — which we term the migration floor: ${\displaystyle {\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 ${\displaystyle \log \nu }$ of ${\displaystyle \mathcal{M}\equiv \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 ${\displaystyle \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 near-extinctions. On a single island, at any given moment, ${\displaystyle \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 $\log (\nu )$ down to the migration floor, only occasionally fluctuating substantially lower (Figure 1).

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, ${\displaystyle s_{\mathit{\text{M}}}}$, can also be correlated with ${\displaystyle 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$).

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, ${\displaystyle s_i}$, these also contribute to variations in average growth rates. The variations in the ${\displaystyle s_i}$ within the community have substantial effects over a time $\sim 1/{\sigma }_s$ with ${\sigma }_s$ roughly the width of the ${\displaystyle 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-⟨\mathrm{{\rm Y}}⟩$, with angular brackets denoting a time average. As introduced earlier, the bias of strain $i$ is its average growth rate at low abundance in the absence of migration (Figure 1B). As we shall see, the size of the community is limited by the condition that the inter-strain variation in general fitness is no larger than variation in average drive from interactions. This means that the biases are of order $1/\sqrt{L}$.

The local population of each strain undergoes wild fluctuations over a logarithmic range $\mathcal{M}=\log (1/m)$ which is quite large. During blooms, the instantaneous growth and decay rates of local populations are substantially larger than the systematic biases (‘Appendix 1’) and change rapidly from growth to decay as seen in Figure 1. The time for abundances to fluctuate from large to small — the duration of blooms — is of order $\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.

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 ${\displaystyle 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’.

Figure 2

Download asset Open asset

The observations in Figure 2 illustrate one of the crucial findings of this work: spatiotemporally chaotic ecological dynamics can allow — but do not guarantee — gradual strain-level diversification up to arbitrarily high number of strains. The behavior depends on the symmetry parameter $\gamma$, which must be substantially negative for the STC to exist. Figure 2B shows that the average rate of diversification, $U\equiv ⟨dL/dT⟩$, 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 ${\displaystyle 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.

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.

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 ${\displaystyle s_i}$ between strains and show how these affect the evolutionary dynamics.

Exponential distribution of ${\displaystyle s_i}$

We first analyze the simplest case: exponentially distributed selective differences with scale $\mathrm{\Sigma }$ and probability density ${\displaystyle \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 ${\displaystyle s_i}$ of the community will change, and we are particularly interested in the dynamics of the population-weighted mean $\widehat{s}={\sum }_j{\overline{\nu }}_j{s}_j$.

The width of the distribution of ${\displaystyle 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 ${\displaystyle 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 ${\displaystyle s_i}$ would dominate over the drives and only the strains with the highest and quite similar ${\displaystyle s_i}$ would survive. Thus $L\gg 1/{\mathrm{\Sigma }}^2$ seems inconsistent. Even for the initial community with ${\displaystyle L_0}$ strains, we expect that ${\displaystyle 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, ${\displaystyle 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 $\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 ${\displaystyle L_0}$ when $K\gg 1/{\mathrm{\Sigma }}^2$, the $\mathrm{\Sigma }$-dependence of the steady-state $L$, and the linear increase of $\widehat{s}$ with number of successful invasions.

Figure 3

Download asset Open asset

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

When strain $i$ has very low abundance, its effects on the others are very small and the forces of the others on it, ${\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 ${\displaystyle \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 single-strain problem, with self consistency conditions on the bias correlations and response function.

With the dynamical mean field understanding of an assembled STC phase in hand, we can proceed to describe the evolutionary process in terms of the distributions of properties of the extant and newly invading strains — in particular their biases and consequent mean abundances. The distribution of biases is perturbed by the introduction of new strains and this can push some of the extant biases below ${\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 mean-drive part of the bias; $\mathcal{L}$ scales with the actual $L$ but discounts strains that are close to extinction.

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 ${\displaystyle \mathcal{N}(\xi )}$ and hence $\mathcal{L}$.

Strain $A$ can successfully invade the community and persist if ${\displaystyle {\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 one-by-one introduction of new strains, the bias of each extant strain undergoes some kind of random walk, and the strain goes extinct if its bias ventures below ${\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.

Figure 4

Download asset Open asset

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 larger-than-average (but still negative) bias.

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

It is convenient to work with the mean drives, ${\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: ${\displaystyle \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 coarse-grains over a substantial range of $T$ (but with range much smaller than $L$). From our analysis in ‘Appendix 8’, we see that $B$ and $D$ are order-unity coefficients which respectively characterize the average and mean-squared response of the bias to the invasion of a new strain. Both are proportional to $\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 mean-square response of the system to a random perturbation applied simultaneously on all the strains — precisely the effect of a successful invasion. The Langevin equation for the drives must be supplemented by a boundary condition that if $\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 eigenvalue-like condition which determines whether the diversification rate of the community is negative or positive, and that the coefficients $B$ and $D$ play a role in determining whether the community diversifies or not. This is consistent with our numerical results, which show that certain parameter regimes allow diversification and other regimes do not — even in the absence of general fitness differerences.

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.

Figure 5

Download asset Open asset

However, the density of biases near ${\xi }_{\mathrm{c}}$ determines the response of the community to evolutionary perturbations, since these low-bias strains are the ones most susceptible to extinction. In particular, the predicted linearly vanishing density of biases determines the distribution of the number, $\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 }$, ${\displaystyle \mathcal{P}(\ell )\sim e^{-\beta \ell }/\sqrt{\ell }}$ (with $\beta$ an order-unity coefficient) which is close to exponential as observed in Appendix 3—figure 1B.

A similar analysis for the initial randomly-assembled community shows that for fixed ${\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 ${\displaystyle 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 ${\displaystyle L_0}$ strains — with a relatively small coefficient — to go extinct before $L$ starts steadily increasing, and this will occur over of order ${\displaystyle 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 ${\displaystyle \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 close-to-marginal bias. In Figure 5A, we see that the mean abundances in an evolved community are more narrowly distributed than in an assembled community, with both fewer highly abundant and fewer rare strains. This is consistent with our picture of the bias distribution being smoothed out toward ${\xi }_{\mathrm{c}}$ due to invasion-triggered extinctions, resulting in the depletion of low-abundance strains. The depletion of abundant strains is likely due to the kill-the-winner dynamics which rewards invading strains that push the most abundant extant strains down.

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

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 ${\displaystyle \mathcal{P}(s)}$ decays faster or slower than exponentially. A heuristic understanding of how the dynamics of $L$ depend on ${\displaystyle \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 ${\displaystyle s_i}$ and cannot decrease further. Thereafter, the shape of ${\displaystyle \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 ${\displaystyle s_i}$ play little role and the population-weighted mean fitness, $\widehat{s}$, only increases gradually. But once $\widehat{s}$ is a few times $\mathrm{\Sigma }$, the tail of ${\displaystyle \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 ${\displaystyle 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 ${\displaystyle \mathcal{P}(s)\sim e^{-s/\hat{\mathrm{\Sigma }}}}$ with the effective width of the extant ${\displaystyle 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 ${\displaystyle 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 ${\displaystyle s_i}$ of both extant strains and of potentially-successful invaders. This suggests that understanding the general behavior at long evolutionary times can be built on understanding the case of the simple exponential distribution ($\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 ${\displaystyle 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 $\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 ${\displaystyle s_i}$ is pushed into the tail of ${\displaystyle \mathcal{P}(s)}$. To crudely take into account the effect of a large initial number of strains when $K\gtrsim 1/{\mathrm{\Sigma }}^2$, $\log (T)$ can be replaced by $\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 $\log [(T+K){\mathrm{\Sigma }}^2]$ prediction by only an ${\displaystyle \mathcal{O}(1)}$ constant, as expected.

Figure 6

Download asset Open asset

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}{(\log T)}^{3-2/\psi }$. Nevertheless, for ${\displaystyle \psi >1}$, the number of strains grows without bound albeit as a sub-linear power of the cumulative number of successful invasions. The average number of extinctions per successful invasion gradually decreases towards one as ${\displaystyle \mathcal{P}(s)}$ in the pertinent range, $s\approx \widehat{s}$, becomes closer and closer to exponential.

For longer-than-exponential tails, ${\displaystyle \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.

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 ${\displaystyle d}$-dimensional array of islands with nearest-neighbor migration; a spatial continuum with diffusive dispersal; or a mixture of long and short distance dispersal events as driven by wind, ocean currents, or hitchhiking on migrant animals (Hallatschek and Fisher, 2014). With real spatial structure, local sub-populations are much more prone to extinction and cannot be as readily rescued by migration from another location where the strain is blooming. Thus, in contrast to the regime we have worked in for this paper, recovery from local extinctions must play a crucial role. The dynamics of invasions, extinctions, and repopulation is very different than in the spatial mean field model: if the underlying dynamics is diffusive, invasion and repopulation will occur by propagating Fisher-Kolmogorov-Petrovsky-Piscunov (FKPP) fronts (Fisher, 1937). The properties of FKPP waves are known to be highly sensitive to dynamics at the wavefront, and the effects of demographic fluctuations have been investigated (Korolev et al., 2010). But the approximately multiplicative “noise” from the ecological interactions will surely change this, and even for a single wave understanding the impact of these larger fluctuations is still an open question (Rocco et al., 2002; Rocco et al., 2000).

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

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

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

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

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

The primary timescale that governs the interplay of evolution and ecology is the ecological relaxation time which is ${\displaystyle \mathcal{O}(\mathcal{M}L)}$, with $\mathcal{M}\equiv \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 ${\displaystyle \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 self-organized chaotic state.

The time to go extinct for a strain destined to do so depends logarithmically on the extinction threshold $1/N$, but as long as $Nm$ is very large, whether extinctions occur is not strongly dependent on $N$, as analyzed in PAF. In our simulations we choose for convenience ${\displaystyle Nm\gg 1}$ so that $\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 ${\displaystyle \sim N{\tau }_{\mathit{\text{gen}}}}$ with ${\displaystyle {\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 ${\displaystyle N\overline{\nu }m{\tau }_{\mathit{\text{gen}}}\gg 1}$, there are many migrants arriving and the dynamics is essentially deterministic. If the island-average $\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 ${\displaystyle Nm\gg L}$ for which the population dynamics of the persistent strains are not strongly influenced by the stochastic or migratory demographic fluctuations.

How do the statistics of the interaction matrix $V_{ij}$ change when conditioned on evolutionary history? How does this compare to conditioning on non-extinction in an assembled community? We ran a set of 10 replicates without general fitness differences over a period where the community grows from 50 initial strains to $\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 ${\displaystyle 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 cross-diagonal 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 $\tilde{\gamma }$, ${\tilde{\sigma }}_V$, ${\tilde{\mu }}_V$ and $\tilde{\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 high-abundance 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 ${\displaystyle \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 fixed-point phase where niche interactions — large negative diagonal terms in $V$ — stabilize the diversity instead of spatial structure. Bunin observed that, conditional on assembly, correlations between the effect of strains $k$ and $l$ both on strain $i$ become negative, that is ${\displaystyle \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 ${\displaystyle \mathcal{O}(1/K)}$ shift in the mean of the interaction matrix toward less competition, which is consistent with ${\displaystyle {\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.

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

Here, the correlation function $C(t,t^{\prime })$ of ${\zeta }_i(t)$, the response function $R(t,t^{\prime })$, and the island-averaged abundance ${\overline{\nu }}_i$, must be determined self-consistently. In order to use this description for the evolved (in addition to assembled) communities we have assumed that the statistics of the $V$ do not change substantially in evolved communities: we show in ‘Appendix 4’ that our simulation results are consistent with this Ansatz. Here, we have made the time-dependence in ${\zeta }_i(t)$ explicit — but elsewhere in the paper when we refer to the time-averaged 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 well-behaved time averages equal to island averages, denoted by overbars. The average growth rate of strain $i$ on a single island excluding migration is

where the static susceptibility, $X$, is the time-integrated total response function ${\displaystyle 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, $⟨\frac{d\log {\nu }_i}{dt}⟩=0$; therefore ${\displaystyle {\overline{r}}_i<0}$ must be exactly compensated by the average effect of the migration term $m(⟨{\overline{\nu }}_i/{\nu }_i⟩-1)$, which is always positive (by the Cauchy-Schwarz inequality) and can be much larger than $m$ when ${\nu }_i$ is small.