Theory of time-averaged neutral dynamics with environmental stochasticity

Competition is the main driver of population dynamics, which shapes the genetic composition of populations and the assembly of ecological communities. Neutral models assume that all the individuals are equivalent and that the dynamics is governed by demographic (shot) noise, with a steady state species abundance distribution (SAD) that reflects a mutation-extinction equilibrium. Recently, many empirical and theoretical studies emphasized the importance of environmental variations that affect coherently the relative fitness of entire populations. Here we consider two generic time-averaged neutral models; in both the relative fitness of each species fluctuates independently in time but its mean is zero. The first (model A) describes a system with local competition and linear fitness dependence of the birth-death rates, while in the second (model B) the competition is global and the fitness dependence is nonlinear. Due to this nonlinearity, model B admits a noise-induced stabilization mechanism that facilitates the invasion of new mutants. A self-consistent mean-field approach is used to reduce the multispecies problem to two-species dynamics, and the large-$N$ asymptotics of the emerging set of Fokker-Planck equations is presented and solved. Our analytic expressions are shown to fit the SADs obtained from extensive Monte Carlo simulations and from numerical solutions of the corresponding master equations.

Figure 1

$P\left(x\right)$, the chance of finding the $A$-type at relative abundance $x$, is plotted for a two competing species system with one-sided mutation, environmental stochasticity, and demographic noise. In both figures $N=1000, \nu =0.01$, and the main panels are plotted using a double logarithmic scale. Results shown include those obtained from a Monte Carlo simulation (green circles), numeric solutions for the steady state of the master equations (B1, B2, B3) (red diamonds), and the analytic prediction of Eq. (8) (black line). In panel (a) the results are depicted for $\delta =0.5$ and $\gamma =0.2$, such that $\nu =g$. In panel (b) $\delta =1.25$ and $\gamma =0.4$ so $\nu =0.1g$. As discussed in the main text, when $\nu <g/2$ there is a peak at high values of $x$. To emphasize this peak we have added an inset where the same results are shown using a linear scale. The fit between the three curves is quite good.

Figure 2

$P_{\mathrm{modelB}}\left(x\right)$, the chance of finding the $A$ type at relative abundance $x$, is plotted for a system with two competing species with one-sided mutation, environmental stochasticity, and demographic noise. In both figures $N=1000, \nu =0.005$. Results shown include those obtained from a Monte Carlo simulation (filled circles), numeric solution for the steady state of the master equations, and the analytic prediction of Eq. (18) (lines with different colors). In panel (a) (plotted using a double logarithmic scale) the results are depicted for large value of $\delta , \delta =2$, so the outcomes imitate those obtained for model A, in particular the two power laws when $\gamma =0.2$ (green circles) and the peak close to $x=1$ when $\gamma =0.4$ (blue circles). In panel (b) (where the scale we have used is semilogarithmic) $\gamma =0.4$ while $\delta =0.4$ for the blue circles and $\delta =0.1$ for the greens. Since $\delta$ is small, the peak at $x=1/2$ is pronounced, and it becomes even steeper as $\delta$ decreases. Here, and in all other figures, the marker sizes were chosen to allow one to distinguish between the three data sets, but the actual width of the lines is smaller. We did not track the standard deviation associated with the outcomes of our Monte Carlo simulations, but the almost perfect agreement between the lines obtained using three different techniques implies that the corresponding error bars are tiny.

Figure 3

The species abundance distribution, $P\left(x\right)$ as a function of $x$, for a time-averaged neutral model (model A) with environmental stochasticity and demographic noise. In both figures $N={10}^4$ and $\nu =0.01$, and the results are plotted using a double logarithmic scale. The outcomes of a Monte Carlo simulation (green circles), numeric solution for the steady state of the master equations (see Appendix pp3) (red diamonds) and the analytic predictions of Eq. (21) (black line) are compared, and the fit is very good. In panel (a) the results are depicted for $\delta =0.5$ and $\gamma =0.5$, such that $\nu /g=0.16$. In panel (b) $\delta =0.25$ and $\gamma =0.2$ so $\nu /g=2$. In both panels the small $x$ behavior is $P\sim 1/x$, but in panel (a) this regime is very narrow since it requires $x\ll 1/G=1.6\times {10}^{-3}$. The $Gx\gg 1$ behavior obeys a power law in panel (a), where the environmental stochasticity dominates ($\nu /g<1/2$) and is dominated by an exponential decay in panel (b), where the mutation losses are stronger. In these parameters, $N$ is not big enough to justify the use of the asymptotic value $\zeta =1$. Instead, the value of ${\overline{f}}_{+}$ used in Eq. (21) was obtained by measuring the long-term average fraction of individuals in the plus state through the MC simulations.

Figure 4

The SAD, $P\left(x\right)$ vs $x$, for model B. In both figures $N={10}^5$ and $\nu =0.001$, and the results are plotted using a double logarithmic scale. The outcomes of a Monte Carlo simulation (blue circles), numeric solution for the steady state of the master equations (red line), and the analytic predictions of Eq. (32) (yellow line) are compared. In panel (a) the results are depicted for $\delta =0.5$ and $\gamma =\sqrt{0.025}$, such that $\nu /g=0.16$. In panel (b) $\delta =0.1$ and $\gamma =0.05$ so $\nu /g=8$. The values of $f_{+}$ and $\sigma$ were taken from the Monte Carlo simulations.