Spatiotemporal ecological chaos enables gradual evolutionary diversification without niches or tradeoffs

  1. Aditya Mahadevan
  2. Michael T Pearce
  3. Daniel S Fisher  Is a corresponding author
  1. Department of Physics, Stanford University, United States
  2. Department of Applied Physics, Stanford University, United States
16 figures, 5 tables and 1 additional file


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

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

Evolution of number of strains without general fitness differences.

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

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

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

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

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

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

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

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

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

Appendix 2—figure 1

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


Table 1
Definitions of commonly used quantities.
STCSpatiotemporally chaotic state
KNumber of strains put into the initial assembled community
VMatrix of pairwise strain interactions
γSymmetry parameter of the interaction matrix; E[VijVji]=δij+γ
INumber of islands
NPopulation size on each island, fixed to be constant
νi,αFractional abundance of strain i on island α
ν¯iTime (or space) average of strain i abundance
νfloorMigration floor mν¯L/M: ~ lower range of local abundances
ΥαLagrange multiplier maintaining iνi,α=1; Υα=iνi,α(si+jVijνj,α)
siGeneral fitness of strain i
P(s)Probability distribution of the si
ΣCharacteristic scale of the s distribution
ψExponent characterizing tail of of P(s)exp[(s/Σ)ψ/ψ]
mMigration rate between islands
MRange of fluctuations in logν; M=log(1/m)
ρCorrelation between parent’s and mutant’s interactions with other strains
TEvolutionary time in epochs, equal to number of attempted invasions
ZNumber of successful invasions
LNumber of extant strains at any point in the evolution
L0Number of strains surviving in initial assembled community
UAverage diversification rate; U=dL/dT
s^Mean general fitness of extant strains;  s^=jνjsj
Σ^The scale of P(s) of extant strains; Σ^=[ddslogP(s)]1|s=s^
ξiBias of strain i, its growth rate at low abundance without migration; ξi=ζi+si-Υ¯, scales as  1/L
N(ξ)Mean abundance of a strain as a function of its bias
ξcCritical bias (negative in the STC) below which strains go extinct, scales as 1/L
ζiMean drive on strain i by other strains in its absence; ζi=jVijν¯ji
LEffective number of extant strains; L=1/jν¯j2; scales with  L
χiStatic response of strain i to perturbations; χi=dν¯idξi
XTotal static response;  X=iχi
ΞFragility of the community to perturbations;  Ξ=jχj2/(1-jχj2)
Appendix 4—table 1
Definition of some statistics of the interaction matrix of a community. The standard deviation is calculated as Stdevi[{xi}]=1Li=1L(xi-1Lj=1Lxj)2 (the population standard deviation) and Corr[X,Y] is the Pearson correlation coefficient calculated for a sample (the covariance normalized by the product of standard deviations).
Appendix 4—table 2
Means (up to corrections due to diagonal entries) and standard deviations of the chosen statistics from 500 realizations of the matrix V from the original ensemble with γ=-0.8, and two values of L pertaining to the different evolutionary conditions. Here we do not condition on either non-extinction or evolution.
mean (L = 500)−0.8100.0446
standard deviation (L = 500)0.0010.0020.0010.0015
mean (L = 150)−0.8100.0816
standard deviation (L = 150)0.00360.0060.00310.0046
Appendix 4—table 3
Statistics of the interaction matrix for both evolved and assembled communities. Evolved communities grew to 500 strains from an initial K=50 and assembled communities had L0500 after a single epoch of duration 12ML, both with si=0. The last row is from a simulation with exponential P(s) and Σ=0.07, which results in L150 at steady state. The mean of each statistic is reported across 10 simulation replicates. In Appendix 4—table 4 we show these statistics appropriately normalized with respect to the original V ensemble. Bold entries correspond to those more than 3 standard deviations away from the original ensemble mean.
assembled (Σ=0)−0.80141.00120.00260.0337
diversifying (Σ=0, ρ=0)−0.80151.00280.00610.0289
diversifying (Σ=0, ρ=0.95)−0.80821.04270.14220.0328
steady state (Σ=0.07, ρ=0)−0.80741.02050.01650.0977
Appendix 4—table 4
Properties of the interaction matrix for evolved and assembled communities, displayed in terms of number of standard deviations from the mean in the original V ensemble. These are the same data as in Appendix 4—table 3, but normalized according to appropriate scale of deviations calculated from original V ensemble. Numbers in boldface have magnitude greater than 3, and are thus clearly statistically significant.
assembled (Σ=0)−−7.2
diversifying (Σ=0, ρ=0)−−10.4
diversifying (Σ=0, ρ=0.95)−8.221.4142.2−7.9
steady state (Σ=0.07, ρ=0)−

Additional files

Download links

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

Downloads (link to download the article as PDF)

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

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

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