Model reduction methods for population dynamics with fast-switching environments: Reduced master equations, stochastic differential equations, and applications

We study stochastic population dynamics coupled to fast external environments and combine expansions in the inverse switching rate of the environment and a Kramers–Moyal expansion in the inverse size of the population. This leads to a series of approximation schemes, capturing both intrinsic and environmental noise. These methods provide a means of efficient simulation and we show how they can be used to obtain analytical results for the fluctuations of population dynamics in switching environments. We place the approximations in relation to existing work on piecewise-deterministic and piecewise-diffusive Markov processes. Finally, we demonstrate the accuracy and efficiency of these model-reduction methods in different research fields, including systems in biology and a model of crack propagation.

Figure 1

Schematic overview of the different model-reduction schemes for populations coupled to external environments with discrete states. Each column and row corresponds to a successive layer of approximation.

Figure 2

Spectra of fluctuations from direct simulations of the full model (closed symbols) and from the linear-noise approximation of the effective SDE capturing intrinsic and extrinsic noise ($k_0=k_1=1$, $\mathrm{\Omega }=20$, $\lambda =20$, ${\alpha }_0=0$, ${\alpha }_1=1$, and $\mathrm{\Delta }t=0.1$).

Figure 3

Stationary probability distribution of the populations of mRNA and protein molecules for the full model in Sec. 4a and the eight levels of model reduction in Fig. 1: (a) full model, (b) piecewise-diffusive process, (c) piecewise-deterministic Markov process, (d) reduced master equation, (e) SDE with switching noise and demographic noise, (f) SDE with switching noise, (g) master equation with average rates, (h) SDE with demographic noise, and (i) rate equation ($▴$ represents a $\delta$ peak). The parameters are $N=2$, $\mathrm{\Omega }=50$, $b_0=b_1=1$, $b_2=20$, $d=9.2$, $\beta =50$, $\delta =1$, $k_{-}=0.025$, $k_{+}=1$, and $\lambda =1250$.

Figure 4

Stationary distribution for the genetic circuit with exclusive binding (Sec. 4b) from numerical integration of (a) the full master equation with explicit environment, (b) the reduced master equation (C4), and (c) the adiabatic approximation which considers average rates. Also shown is (d) the marginal distribution of $n_A-n_B$ in order to compare the three distributions. The parameters are ${\alpha }_0={\beta }_0=0$, ${\alpha }_1={\beta }_1=1$, $\mathrm{\Omega }=20$, $\lambda {\kappa }_0=\lambda {\kappa }_1=20$, and $\gamma =\delta =1$.

Figure 5

(a) Sample path of the model of crack growth (Sec. 4c). Background shading indicates the state of the environment, with states 0, 1, and 2 shown progressively darker. (b) Survival probability as a function of time. The black line shows results from Monte Carlo simulations and the dashed line is the prediction of Eq. (43). (c) and (d) Same quantities as in (a) and (b), respectively, for tenfold increased switching rates. The model parameters are given in the text.

Figure 6

(a) Schematic illustrating an environment with two distinct conditions, with $N$ and $M$ identical stages, respectively. (b) Stationary distribution for different values of $N=M$. Histograms show the results of simulations and lines are from the theory described in the text. (c) Variance of the stationary distribution as a function of $N$ (again for the case $N=M$). The black line shows the results of the theory and orange circles are from simulations of the piecewise-deterministic Markov process. The parameters $N$ and $M$ have been generalized to include noninteger values, by considering $\gamma$-distributed waiting-time distributions in the two conditions (see the text). The dashed line shows the variance of the limit cycle obtained in the limit of a periodic environment. The parameters are $b_0=100/3$, $b_1=500/3$, and $\lambda k_0=\lambda k_1=20$.