Many-Species Ecological Fluctuations as a Jump Process from the Brink of Extinction

Highly diverse ecosystems exhibit a broad distribution of population sizes and species turnover, where species at high and low abundances are exchanged over time. We show that these two features generically emerge in the fluctuating phase of many-variable model ecosystems with disordered species interactions, when species are supported by migration from outside the system at a small rate. We show that these and other phenomena can be understood through the existence of a scaling regime in the limit of small migration, in which large fluctuations and long timescales emerge. We construct an exact analytical theory for this asymptotic regime that provides scaling predictions on timescales and abundance distributions that are verified exactly in simulations. In this regime, a clear separation emerges between rare and abundant species at any given time, despite species moving back and forth between the rare and abundant subsets. The number of abundant species is found to lie strictly below a well-known stability bound, maintaining the system away from marginality. At the same time, other measures of diversity, which also include some of the rare species, go above this bound. In the asymptotic limit where the migration rate goes to zero, trajectories of individual species abundances are described by non-Markovian jump-diffusion processes, which proceeds as follows: A rare species remains so for some time, then experiences a jump in population sizes after which it becomes abundant (a species turnover event) and later sees its population size gradually decreasing again until rare, due to the competition with other species. The asymmetry of abundance trajectories under time reversal is maintained at a small but finite migration rate. These features may serve as fingerprints of endogenous fluctuations in highly diverse ecosystems.

Popular Summary

Ecosystems often harbor many coexisting species whose population sizes may undergo complex high-dimensional dynamics driven by the interactions between them. Strikingly, such dynamics feature large fluctuations in population sizes, a phenomenon observed in nature and computational models. The understanding and fundamental properties of these dynamics have remained elusive, including the distribution of population sizes, the timescales involved, and the number of coexisting species. Here, we present an analytical framework for this dynamical state that allows us to address these different questions.

Our approach builds on dynamical mean-field theory. This powerful description exactly maps the collective deterministic dynamics to independent stochastic evolutions for the dynamics of the population sizes. When we couple the ecosystem to its surroundings through migration of individuals at a small rate, we find that long timescales and large fluctuations reinforce each other. We characterize this behavior by deriving an exact effective theory over long timescales. Notably, we find that the population size of any species alternates between periods of rarity and abundance, with characteristic sawtooth trajectory profiles.

This work makes experimentally testable predictions and may be instrumental in understanding ecosystem dynamics in broader contexts, for instance, for spatially structured ecosystems or under ecoevolutionary dynamics.

Figure 1

Species richness and species abundance distributions at fixed points (left column) and a persistently fluctuating state (right column). (a) Dynamics reaching a fixed point. (b) Persistent dynamics, where the abundances fluctuate indefinitely in the range $\lambda \lesssim N_i\lesssim 1$. (c)–(f) The corresponding distributions for $N$ and $\ln N$ denoted by $P(N)$ and $\widehat{P}(\ln N)$, respectively. The distributions contain three parts [see (e) and (f)]: a top part, which remains $O({\lambda }^0)$ when $\lambda \to 0^{+}$, containing $S_{\mathrm{top}}^{*}$ species, an intermediate part at $\lambda \ll N\ll 1$, and a low part at $O(\lambda )$. The intermediate part contains a finite fraction of the species in the dynamical phase but not in the fixed-point phase. (g) The three definitions of species richness normalized by $\mathrm{var}(\alpha )$ so that 1 is the fixed-point stability bound as a function of the interaction variability $\sigma$. For $\sigma <{\sigma }_c$, a fixed point is reached, and all three definitions coincide and lie below the stability bound and agree with known theory (solid line). At $\sigma ={\sigma }_c$, the stability bound is reached. Beyond it, there are persistent fluctuations, and the three definitions no longer coincide: $S_{\mathrm{top}}^{*}$ is lower than the stability bound, while $S_{\text{inter}}^{*}$ and $S_{\text{growth}}^{*}$ are above it. The diversities are obtained by solving the rescaled dynamics defined in Sec. 4 and agree quantitatively with careful analysis of the full equations of motion (see Sec. 7). Simulation parameters for all figures are given in Appendix pp1-s7.

Figure 2

Species dynamics and correlations. (a) Species are initialized with order-one values. We follow $\ln N_i(t)$. In the transient regime, the changes in $\ln N_i(t)$ are comprised of three elements: downward motion, roughly in a straight line (maintaining its slope for a large part of the decay between largest and smallest values, corresponding to exponential decay of $N_i$), upward motion, roughly in a straight line, and slow changes at high population values $N_i(t)\simeq O(1)$. The dynamics slow down, with the time spent in each of these elements increasing with the time since the start of the dynamics; see Sec. 6. As the excursions grow longer, the average value of $\ln N_i(t)$ decreases linearly in $t$. This transient regime ends when a finite fraction of species reach the migration floor $N_i(t)\simeq O(\lambda )$ after a time of order $|\ln \lambda |$. At long times, after the transient regime is over, $\ln N_i(t)$ performs the three dynamical elements described above, and a fourth one, slowly changing around the migration floor. The vertical red dashed line marks the crossover between these two regimes. (b) At long times, timescales do not grow in time anymore, as seen in the correlation function that depends only on time differences $C(t,t^{\prime })=C(t-t^{\prime })$. (c),(d) This timescale is proportional to $|\ln \lambda |$. This is demonstrated by the data collapse when plotting $C$ against $t/|\ln \lambda |$ presented in (c), compared with the data without rescaling shown in (d).

Figure 3

The rescaled dynamics. Trajectories of $N$ and $z\equiv \ln (N)/|\ln \lambda |$ as a function of the rescaled time $s\equiv t/|\ln \lambda |$ in an example run of the limiting rescaled dynamics when $\lambda \to 0^{+}$.

Figure 4

Collapse of species abundance distributions. Numerically measured distributions of $N$ (a) and $z=\ln N/|\ln \lambda |$ (b) converge to the distributions predicted by the rescaled process as $\lambda \to 0^{+}$.

Figure 5

Growth of timescales in dynamics without migration. (a),(b) When $\lambda =0$, the collapse of $C(t,t+\tau )$ as a function of $\tau /t$ demonstrates the linear growth of timescales with the elapsed time. Data without rescaling, as a function only of $\tau$, are shown for comparison in (b). (c),(d)Crossover at finite $\lambda$ from a transient regime exhibiting the $\lambda =0$ phenomenology to the long-time behavior, as identified by the time constant ${\tau }_2$ of the correlation function. After an initial transient of finite time duration, ${\tau }_2$ grows linearly with the elapsed time and independently of $\lambda$, a trademark of the $\lambda =0$ dynamics. The crossover to the long-time behavior happens around times proportional to $|\ln \lambda |$, after which ${\tau }_2$ stabilizes to a value proportional to $|\ln \lambda |$. This scaling of the crossover is manifest in the data collapse in (c), with the data without rescaling displayed in (d) for comparison. ${\tau }_2(t,\lambda )$ is defined as the time when $[C(t,t+{\tau }_2)-C(t,\infty )]/[C(t,t)-C(t,\infty )]=e^{-1}$. The value of $C(t,\infty )$ is estimated through an exponential fit of the function $C(t,t+\tau )$ as a function of $\tau$.

Figure 6

Average top diversity measures, ${\varphi }_{\mathrm{top}}^N(S,\lambda ,\epsilon ),{\varphi }_{\mathrm{top}}^g(S,\lambda ,\epsilon )$, as a function of $S$ and $\lambda$. ${\varphi }_{\mathrm{top}}^N$ is the fraction of species with $N_i>\epsilon ={10}^{-3}$; ${\varphi }_{\mathrm{top}}^g$ also requires $g_i>0$. (a) ${\varphi }_{\mathrm{top}}^N(S,\lambda ,\epsilon )$ as a function of $S$ converges as ${\varphi }_{\mathrm{top}}^N(S\to \infty ,\lambda ,\epsilon )+\#/S$, at fixed $\lambda$. Data points at $S=500$ are excluded from the linear fit. (b) Large $S$ value of the top diversity ${\varphi }_{\mathrm{top}}^N(S\to \infty ,\lambda ,\epsilon )$ converges as ${|\ln \lambda |}^{-1/2}$ to its asymptotic value ${\varphi }_{\mathrm{top}}^N(S\to \infty ,\lambda \to 0^{+},\epsilon )$. The asymptotic value ${\varphi }_{\mathrm{top}}^N(S\to \infty ,\lambda \to 0^{+},\epsilon )$ is within 5% of the value predicted by the rescaled DMFT and well below the stability bound. For all reasonable values of the migration rate (here already when $\lambda \gtrsim {10}^{-40}$), the measured top diversity is significantly above the stability bound. ${\varphi }_{\mathrm{top}}^N(S\to \infty ,\lambda ,\epsilon )$ is inferred from the finite-$S$ scaling, panel (a). (c),(d) Same as panels (a) and (b) but for ${\varphi }_{\mathrm{top}}^g(S,\lambda ,\epsilon )$. The convergence to the asymptotic ${\varphi }_{\mathrm{top}}$ is more rapid; in particular, ${\varphi }_{\mathrm{top}}^g$ can be either below or above the stability bound.

Figure 7

Asymmetry of blooms under time reversal. Trajectories of $N(t)$, while at high abundance display a clear asymmetry in time, with a rapid initial increase and more gradual decrease. This is conspicuous even for migration rates that are not very small (here, $\lambda ={10}^{-6}$). It is the finite-$\lambda$ counterpart of the asymptotic behavior at $\lambda \to 0^{+}$ featuring sharp jumps from zero to positive $N$; see Fig. 3. Shown are trajectories that go above $N={10}^{-2}$ at some time $t_{\mathrm{in}}$ and stay above it until time $t_{\mathrm{in}}+\mathrm{\Delta }t$, and reach $N\ge 0.1$ at some intermediate time. Here, $\mathrm{\Delta }t\simeq 5.8|\ln \lambda |$ that is the most common length of such trajectories. In light gray are 22 example trajectories, and the thick blue line shows the average over many such trajectories.

Figure 8

Enlargement of the small-$s$ behavior of $\widehat{C}(s)$. The correlation function is obtained by numerically solving the DMFT equations (5), (7), and (8). The derivative is finite at $s=0^{+}$, thereby confirming the absence of fast timescale in the Lotka-Volterra dynamics.

Figure 9

Collapse of species abundance distributions. Numerically measured distributions of $N$ (a) and $z=\ln N/t$ (b) converge to the distributions predicted by the rescaled process as $t\to \infty$.

Figure 10

Model behavior with partially symmetric or antisymmetric interaction matrices. Even when the interaction matrix possesses some degree of symmetry or asymmetry, ${\tau }_{\lambda }=|\ln \lambda |$ is the only timescale controlling the dynamics when $\lambda \to 0^{+}$. We plot the collapse of the correlation function $C_{\lambda }(t,t+\tau )$ as a function of $\tau /|\ln \lambda |$ for partially antisymmetric [$\gamma =-0.5$, (a)] and partially symmetric [$\gamma =0.5$, (b)] interaction matrices. The top diversity measured at finite $\lambda$ asymptotically goes below the May bound as $\lambda \to 0^{+}$ for partially antisymmetric [$\gamma =-0.5$, (c)] and partially symmetric [$\gamma =0.5$, (d)] interaction matrices. The diversity is evaluated by counting the number of species with $N_i>{10}^{-3}$. The form of the extrapolation in powers of $|\ln \lambda |$ is explained in Appendix pp1-s4. Parameters: (a) $S=5000,\mu =10,\sigma =4$, average over 40 realizations; (b) $S=5000,\mu =50,\sigma =1.3$, average over 40 realizations.

Figure 11

Behavior of the first and second moments of the population size. $⟨N⟩$ and $⟨N^2⟩$ are robust predictions that weakly depend on the number of species $S$ and the migration rate $\lambda$. In the inspected range of parameters, the variations in $⟨N⟩$ are of the order of 1% and those in $⟨N^2⟩$ of the order of 10%. Parameters: $\sigma =1.8$, $\mu =10$, average over 80 realizations.

Figure 12

Model behavior at higher values of $\sigma$, $\mu$. (a) Collapse of the correlation function $C_{\lambda }(t,t+\tau )$ as a function of $\tau /|\ln \lambda |$, showing that the only timescale in the problem scales as $|\ln \lambda |$. Compare with Fig. 2. (b) ${\varphi }_{\mathrm{top}}^N$ as a function of $\lambda$, showing the same dependence on $\lambda$, converging as ${|\ln \lambda |}^{-1/2}$. Compare with Fig. 6. The extrapolated asymptotic value at $\lambda \to 0^{+}$ is well below the May bound. Here, $\sigma =5$, $\mu =50$, $S=5000$. ${\varphi }_{\mathrm{top}}^N$ is defined with $N_i>\epsilon ={10}^{-3}$.

Figure 13

Identifying the top species by combined selection and dynamics. (a) Top species are identified by stopping the dynamics (solid lines), and removing species that are about to grow and disrupt the system, here the orange trajectory, and those with a negative invasion growth rate (by removing species $i$ if $N_i/g_i<1/2$), and then continuing the dynamics (dashed lines). The outcome is an equilibrium. (b) The distribution of abundances thus obtained is close to the asymptotic distribution obtained from the rescaled dynamics. It is much closer than the bare abundance distribution $P(N)$, and is also very different from the distribution in the equilibrium phase. $S=5000$, $\sigma =1.8$, $\mu =10$, $\lambda ={10}^{-10}$.