Neutral dynamics with environmental noise: Age-size statistics and species lifetimes

Neutral dynamics, where taxa are assumed to be demographically equivalent and their abundance is governed solely by the stochasticity of the underlying birth-death process, has proved itself as an important minimal model that accounts for many empirical datasets in genetics and ecology. However, the restriction of the model to demographic $\left[\mathcal{O}\left(\sqrt{N}\right)\right]$ noise yields relatively slow dynamics that appears to be in conflict with both short-term and long-term characteristics of the observed systems. Here we analyze two of these problems—age-size relationships and species extinction time—in the framework of a neutral theory with both demographic and environmental stochasticity. It turns out that environmentally induced variations of the demographic rates control the long-term dynamics and modify dramatically the predictions of the neutral theory with demographic noise only, yielding much better agreement with empirical data. We consider two prototypes of “zero mean” environmental noise, one which is balanced with regard to the arithmetic abundance, another balanced in the logarithmic (fitness) space, study their species lifetime statistics, and discuss their relevance to realistic models of community dynamics.

Figure 1

Variance through time plots for three empirical datasets: (a) Trees of the tropical forest in Barro Colorado Island (BCI), considering trees with diameter $>10\mathrm{mm}$ (data are from www.ctfs.si.edu) from 1990 to 2005 (every 5 y). (b) Fish samples collected from the cooling-water filter screens at Hinkley Point B Power Station in Bristol, from 1980 to 2008 [19]. We only consider estuarine fish species (not crustaceous organisms). (c) Herbaceous plant dataset [20] comprising a series of 51 quadrats of 1 ${\mathrm{m}}^2$ from mixed Kansas grass prairies where all individual plants were mapped every year from 1932 to 1972. In all cases we gathered $\text{ln}(N_{t+\mathrm{\Delta }t}/N_t)$ (natural logarithm was used) from every two same species abundance observations that are $\mathrm{\Delta }t$ apart, and plotted the variance of this quantity vs $\mathrm{\Delta }t$ (the details of this technique are presented in Ref. [15]). $m$ denotes the angular coefficient of the linear fit of data points $[m=4\times {10}^{-3}$ for (a), $5.5\times {10}^{-2}$ for (b), and $0.01$ for (c)], while $P_{val} [0.01$ for (a), $P<{10}^{-4}$ for (b) and (c)] and $R^2 [0.97, 0.72$ and $0.97$ for (a), (b), and (c) correspondingly] give the significance and the goodness of the fit, respectively. In all three cases the variance appears to grow linearly in time, as suggested for a system with logarithmically balanced environmental noise, where the slope $m$ indicates the value of the effective $D$ defined in Eq. (3). The intercept of the linear fit is above zero for the fish and the grassland, indicating the effect of sampling errors [15]. The negative intercept for trees indicates a delayed response of the system to the changing environment. Since the time scales are relatively small, the same analysis for arithmetically balanced noise, where the variance of $(N_{t+\mathrm{\Delta }t}-N_t)/N_t$ is plotted against $t$, yields almost the same results (not shown).

Figure 2

(a) A plot of ${\delta }_t$, the chance to survive at time $t$ for a certain levels of overall diffusion coefficient $D$ [log-log scale, the lines correspond to different $D$ values in the legends are ordered from bottom (low) to top (high $D$) in both panels]. The simulated process is the standard Wright-Fisher process with nonoverlapping generations. (b) Plotting ${\delta }_{t}t$ vs time in linear-log scale. (c) A plot of ${\delta }_{t}t$ against $tD$. In this case the data collapse as predicted by Eq. (B7). This collapse suggests that the proper scaling of the cumulative survival probability is ${\delta }_t=\frac{\mathcal{F}\left(tD\right)}{t}$, with $\mathcal{F}\left(x\right)\sim x^{1/2}$ for $x\to \infty$ and $\mathcal{F}\left(x\right)\sim \text{const}$ when $x\to 0$.

Figure 3

The survival probability of a species $\delta$ vs time in case B (arithmetically balanced noise). The logarithm of $t\delta \left(t\right)$ is plotted against $t$ (inset), showing an exponential decay with $D$ dependent slope. When plotted against $Dt$ (main panel) the data collapse, indicating the dependence of $\mathcal{F}$ on $Dt$.

Figure 4

The mean of the logarithmic (red) and arithmetic (black) abundance fluctuations as a function of the abundance, extracted from simulation of a neutral community model with both demographic and environmental noise. A Moran-type continuous time process with discrete individuals was simulated. At each small time step ($dt=0.1$) we generate the number of births and deaths for each species. The number of births is a Poisson distributed number with mean ${\alpha }_i\left(t\right)dt$, where ${\alpha }_i\left(t\right)$ is the instantaneous growth rate (fitness) of species $i$. The number of deaths in the species is binomially distributed, with the probability $N_{\mathrm{tot}}/K$ per individual, where $N_{\mathrm{tot}}$ is the instantaneous size of the community and $K={10}^5$. In addition, a Poisson number of new immigrants is drawn, with mean $\mu dt$, where $\mu =2$. Each immigrant founds a new species. The birth rates ${\alpha }_i$ are given by ${\alpha }_i=\text{exp}\left({\gamma }_i\right)$ where the ${\gamma }_i$ are generated by an Orenstein-Uhlenbeck process, ${\dot{\gamma }}_i=-{\gamma }_i/\tau +\eta$, so that they are Gaussian distributed with mean 0 and variance 0.001 and correlation time 2. This way the birth rates are guaranteed positive.

Figure 5

The Morse potential $V\left(y\right)=exp\left(2y\right)/4D-exp\left(y\right)/2$ for $D=0.1$.

Figure 6

Numerical simulation of the Fokker-Planck equation (B7), compared to the analytic approximation, Eq. (B12), for $D=0.1, Dt=10$, and $Dt=25$.

Figure 7

A graph of the ensemble-averaged width $\overline{\zeta \left(s\right)}$ calculated by averaging $\zeta$ as computed from the analytic generating function Eq. (B5) over ${10}^6$ realizations of binary noise $\gamma =\pm {\gamma }_0$, with ${\gamma }_0=0.1$. Also shown is the analytic formula, Eq. (B18).