Beyond the adiabatic limit in systems with fast environments:
a $\tau$-leaping algorithm

We look at systems composed of discrete individuals. We will refer to this synonymously as the ‘system proper’, or ‘the population’. Each of the individuals is of one of $S$ species (or types), labelled $i=1,\dots,S$. We write $n_{i}$ for the number of individuals of species $i$ in the population, and $\mathbf{n}=(n_{1},\dots,n_{S})$. The system evolves in an external environment, whose state we write as $\sigma$. These states are time dependent, and can either take discrete values or be continuous.

The dynamics in the population proceeds through reactions $r=1,\dots,R$. Each of these reactions converts a number of individuals from one type into another. Time in the model is continuous, and we assume that the dynamics is Markovian. We then write $R_{r,\sigma}(\mathbf{n})$ for the rate of reaction $r$ if the environment is in state $\sigma$ and the population in state $\mathbf{n}$. The stoichiometric coefficient $\nu_{r,i}$ indicates how the number of individuals of type $i$ changes when a reaction of type $r$ occurs. Each $\nu_{r,i}$ is an integer, which can be negative, zero, or positive. We write ${\mbox{\boldmath$\nu$}}_{r}=(\nu_{r,1},\dots,\nu_{r,S})$. The rates $R_{r,\sigma}(\mathbf{n})$ and the stoichiometric coefficients fully specify the dynamics of the population when the environment is in state $\sigma$.

The state of the environment undergoes a Markovian stochastic process, governed by a master equation if states are discrete or by a stochastic differential equation in the case of continuous environmental states. These dynamics can depend on the state of the population $\mathbf{n}$. If the environmental states are discrete we write $q_{\sigma\to\sigma^{\prime}}(\tau)$ for the probability of finding the environment in state $\sigma^{\prime}$ at a particular point in time, given that $\tau$ units of time earlier it was in state $\sigma$. If the environment is continuos then $q_{\sigma\to\sigma^{\prime}}(\tau)$ is a probability density for $\sigma^{\prime}$ (at given $\sigma$). We call $q_{\sigma\to\sigma^{\prime}}(\tau)$ the transition kernel of the environmental process. We write ${\mbox{\boldmath$\rho$}}^{*}$ for the stationary distribution of the environmental dynamics. For discrete environmental states the entries $\rho^{*}_{\sigma}$ denote the probability of finding the system in state $\sigma$ in the stationary state. For continuous environments $\rho^{*}_{\sigma}$ is the stationary probability density for $\sigma$.

A conventional reaction system (without external environment) is governed by a chemical master equation of the form

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}P(\mathbf{n},t)=$
	$\displaystyle\sum_{r}\big{[}R_{r}(\mathbf{n}-{\mbox{\boldmath$\nu$}}_{r})P(\mathbf{n}-{\mbox{\boldmath$\nu$}}_{r},t)-R_{r}(\mathbf{n})P(\mathbf{n},t\big{)}].$		(1)

The notation is as in Sec. II.1, the only difference is that there is no subscript $\sigma$, as there is no environment. Sample paths entail events (reactions) which can occur at any point in continuous time, separated by exponentially distributed random waiting times. In each such event the state of the system $\mathbf{n}$ changes, and accordingly the reaction rates $R_{r}(\mathbf{n})$ can also change. Sample paths can be generated for example using the celebrated Gillespie algorithm (Gillespie, 1976, 1977).

The $\tau$-leaping algorithm for such conventional reaction systems is built around the idea of keeping reaction rates constant over finite time steps of length $\tau$ Gillespie (2001). That is to say, if the state of the population is $\mathbf{n}$ at time $t$, then the assumption is made that this state $\mathbf{n}$ and the rates $R_{r}(\mathbf{n})$ do not change until the end of the time step. The algorithm does not account for potential changes of the rates as individual reactions occur, and instead directly ‘leaps’ to time $t+\tau$. This is justified provided the so-called ‘leap condition’ is fulfilled Gillespie (2001): broadly speaking the time step $\tau$ must be sufficiently small so that the state $\mathbf{n}$ in the continuous-time system does not change significantly in a time interval of length $\tau$.

Making the approximation of constant $\mathbf{n}$ in the time interval from $t$ to $t+\tau$, the number of reactions of type $r$ that fire in this interval follows a Poissonian distribution with parameter $\tau R_{r}(\mathbf{n})$. Accordingly, realisations of Poissonian random variables $m_{1},\dots,m_{R}$ are drawn, and the corresponding numbers of each reactions are executed simultaneously. This generates a new state $\mathbf{n}^{\prime}$ at time $t+\tau$, with entries $n_{i}^{\prime}=n_{i}+\sum_{r=1}^{R}m_{r}\nu_{i,r}$. The process then repeats with updated rates $R_{r}(\mathbf{n}^{\prime})$.

The idea of the $\tau$-leaping algorithm we introduce for systems in external environments is similar. As in the conventional algorithm we discretise time, and keep the composition of the population $\mathbf{n}$ fixed during each iteration. It is only updated at the end of each step. From now on we use $\Delta t$ for the duration of a step instead of $\tau$.

The difference to the conventional case is the external environment. If the environmental state space is discrete then switches of the environment can in principle be simulated in continuous time along with the other reactions (using Gillespie algorithm (Gillespie, 1976, 1977)). They can also be dealt with by means of the conventional $\tau$-leaping algorithm, again along with the other reactions. These are natural simulation approaches when the environment operates on a similar time scale as the reactions in the population. Not much can then be gained by distinguishing between environmental processes and the dynamics in the system proper.

If the environment is infinitely fast compared to the reactions in the population, then the environment reaches stationarity on very short time scales. One can average over environmental states, see for example Bowen et al. (1963); Segel and Slemrod (1989); Newby and Bressloff (2010); Bressloff and Newby (2014); Hufton et al. (2019b, a). If the environment is discrete, for example, we can use average rates

$$R_{r}^{*}(\mathbf{n})\equiv\sum_{\sigma}\rho^{*}_{\sigma}R_{r,\sigma}(\mathbf{n}).$$

(2)

In the case of continuous environments the sum is to be replaced with an integral. These rates are functions of $\mathbf{n}$ only, the environmental process has been averaged out. Noise from the environmental process plays no role in the dynamics if these average rates are used. This corresponds to making a quasi-stationary state approximation for the fast-moving environment Bowen et al. (1963); Segel and Slemrod (1989).

The aim of this paper is to go beyond this adiabatic limit, and to construct a $\tau$-leaping algorithm which captures some elements of extrinsic noise. We focus on the limit of a fast, but not infinitely fast environmental dynamics.

Broadly speaking the $\tau$FE algorithm is constructed around the idea of treating the reaction rates $R_{r}(\mathbf{n})$ as stochastic variables in each discrete time step. These random variables represent the rates one obtains when averaging the environmental process over the time step $\Delta t$. Assuming that the rate of change of the environment is finite these average rates will remain stochastic. In the limit of infinitely fast environments the deterministic limit in Eq. (2) is recovered, and there is no stochasticity from the environment.

To construct the random reaction rates for each step, we make an approximation: we use a Gaussian distribution for the rates, with means as in Eq. (2) and with variances and correlations derived from the original combined process of the population and environment. We describe this in detail in the next section.

Here we assume the environmental states are discrete, $\sigma\in\{1,\dots,M\}$. The dynamics of the environment is governed by the rates $\lambda A_{\sigma\to\sigma^{\prime}}(\mathbf{n})$ for transitions from $\sigma$ to $\sigma^{\prime}$. The factor $\lambda$ is introduced to control the time-scale separation between reactions in the population and the switching of the environment. We use the notation $\mathbf{A}(\mathbf{n})$ for the $M\times M$ matrix with elements $A_{\sigma\to\sigma^{\prime}}(\mathbf{n})$. We also set $A_{\sigma\to\sigma}(\mathbf{n})=-\sum_{\sigma^{\prime}\neq\sigma}A_{\sigma\to\sigma^{\prime}}(\mathbf{n})$. The combined dynamics of population and environment are then described by the master equation

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}P(\mathbf{n},\sigma,t)=$
	$\displaystyle\sum_{r}\big{[}R_{r,\sigma}(\mathbf{n}-{\mbox{\boldmath$\nu$}}_{r})P(\mathbf{n}-{\mbox{\boldmath$\nu$}}_{r},\sigma,t)-R_{r,\sigma}(\mathbf{n})P(\mathbf{n},\sigma,t)\big{]}$
	$\displaystyle+\lambda\sum_{\sigma^{\prime}}\big{[}A_{\sigma^{\prime}\to\sigma}(\mathbf{n})P(\mathbf{n},\sigma^{\prime},t)-A_{\sigma\to\sigma^{\prime}}(\mathbf{n})P(\mathbf{n},\sigma,t)\big{]}.$		(3)

The rates $\lambda A_{\sigma\to\sigma^{\prime}}(\mathbf{n})$ can depend on the state of the population, $\mathbf{n}$. This means that $\mathbf{n}$ and $\sigma$ do not necessarily evolve in time independently. However, as mentioned above the state $\mathbf{n}$ of the system is kept constant during each $\tau$-leaping step. This in turn means that the transition rates for the environment also remain constant during each step.

We now focus on one such time step, starting at time $t$ and ending at $t+\Delta t$. We assume that $\mathbf{n}$ remains constant during this time interval. For the remainder of Sec. III.1 we suppress the potential dependence of $\mathbf{A}$ on $\mathbf{n}$, although it is always implied. We write $\rho_{\sigma}(t^{\prime})$ for the probability that the environment is in state $\sigma$ at time $t\in[t,t+\Delta t]$. We then have the master equation

$$\frac{\mathrm{d}{\mbox{\boldmath$\rho$}}}{\mathrm{d}t^{\prime}}=\lambda\mathbf{A}{\mbox{\boldmath$\rho$}}$$

(4)

for the environmental dynamics. The stationary distribution ${\mbox{\boldmath$\rho$}}^{*}$ for the environment is the solution of $\mathbf{A}{\mbox{\boldmath$\rho$}}^{*}=0$. If $\mathbf{A}$ depends on $\mathbf{n}$, then $\rho^{*}$ will also be a function of $\mathbf{n}$. Assuming that the environmental process is irreducible this stationary distribution is unique for any one $\mathbf{n}$.

The stochastic matrix $\mathbf{A}$ has one zero eigenvalue, which we write as $\mu_{1}=0$. The remaining eigenvalues are denoted by $\mu_{2},\dots,\mu_{M}$. We then have $\mu_{2},\dots,\mu_{M}<0$. The corresponding (right) eigenvalues are written as $\mathbf{v}_{1}={\mbox{\boldmath$\rho$}}^{*}$ (the eigenvector corresponding to eigenvalue $0$), and $\mathbf{v}_{2},\dots,\mathbf{v}_{M}$ respectively for the remaining eigenvectors. These are all understood to be column vectors of length $M$.

We note that the general solution of Eq. (4) can be written in the form

$${\mbox{\boldmath$\rho$}}(t^{\prime})={\mbox{\boldmath$\rho$}}^{*}+\sum_{\ell=2}^{M}c_{\ell}e^{\lambda\mu_{\ell}(t^{\prime}-t)}\mathbf{v}_{\ell},$$

(5)

with coefficients $c_{\ell}$ determined by the initial condition at the beginning of the time step $t^{\prime}=t$. More precisely these coefficients can be obtained from the linear system

$$\sum_{\ell=2}^{M}c_{\ell}\mathbf{v}_{\ell}={\mbox{\boldmath$\rho$}}(t)-{\mbox{\boldmath$\rho$}}^{*}.$$

(6)

We remark that there are $M-1$ coefficients, $c_{\ell}$ ($\ell=2,\dots,M$). The system in Eq. (6) technically consists of $M$ equations, but these are not independent due to normalisation of the probabilities on the right-hand side.

Calculating the probability $q_{\sigma\to\sigma^{\prime}}(\Delta t)$ to find the environment in state $\sigma^{\prime}$ at the end of the time step if it was in $\sigma$ at the beginning of the step is now mainly a matter of computing the coefficients $c_{\ell}$. We write $c_{\ell,\sigma}$ for the value the coefficient $c_{\ell}$ takes when $\rho_{\sigma^{\prime}}(t)=\delta_{\sigma^{\prime},\sigma}$ (i.e., when the system starts in state $\sigma$ at the beginning of the step).

We then have

$$q_{\sigma\to\sigma^{\prime}}(\Delta t)=\rho^{*}_{\sigma^{\prime}}+\sum_{\ell=2}^{M}c_{\ell,\sigma}e^{\lambda\mu_{\ell}\Delta t}v_{\ell,\sigma^{\prime}},$$

(7)

where $v_{\ell,\sigma^{\prime}}$ is the $\sigma^{\prime}$-entry of the eigenvector $\mathbf{v}_{\ell}$ of $\mathbf{A}$.

If the matrix $\mathbf{A}$ depends on the population state $\mathbf{n}$, the parameters $\mu_{\ell},v_{\ell,\sigma^{\prime}},$ and $c_{\ell,\sigma}$ can also depend on $\mathbf{n}$. For simplicity of notation we have not included this potential dependence in the above equations.

The $\tau$-leaping algorithm proceeds in discrete time intervals of length $\Delta t$. We continue to focus on one such interval $[t,t+\Delta t]$. The state of the population at the beginning of the step is $\mathbf{n}$ and we assume that this state does not change until the end of the interval. We do however take into account the fact that the state of environment $\sigma$ can undergo changes in the interval from $t$ to $t+\Delta t$. As a consequence, $R_{r,\sigma}(\mathbf{n})$ (at fixed $\mathbf{n})$ is also a function of time.

We then introduce the time-averaged quantities

$$\overline{R_{r}}(\mathbf{n})\equiv\frac{1}{\Delta t}\int_{t}^{t+\Delta t}\mathrm{d}t~{}R_{r,\sigma(t)}(\mathbf{n}),$$

(8)

noting that the time average is over the duration of the time step only as opposed to a long-term asymptotic time average. Given that the time step $\Delta t$ is finite ($\Delta t<\infty$) and assuming that the environment fluctuates with finite rates ($\lambda<\infty$), the quantity $\overline{R_{r}}(\mathbf{n})$ is a stochastic variable as it depends on the realisation of the environmental process. In one given time interval, the rates $\overline{R_{r}}(\mathbf{n})$ for different $r$ will be correlated as they all derive from the same path of the environment. As we will show below, the fluctuations of the random variables $\overline{R_{r}}(\mathbf{n})$ in any one time step are inversely proportional to $\lambda\Delta t$ to leading order. In the limit $\lambda\Delta t\to\infty$ the $\overline{R_{r}}(\mathbf{n})$ become deterministic.

We assume that the distribution for $\sigma$ at the beginning of the time step is the stationary distribution ${\mbox{\boldmath$\rho$}}^{*}$. This is the case for example, if then environmental state is drawn from the stationary distribution at the beginning of the simulation. The distribution for $\sigma(t^{\prime})$ is then also the stationary distribution at each time $t^{\prime}\in[t,t+\Delta t]$. Writing $\left\langle{\dots}\right\rangle$ for the average over realisations of the environmental process, we then have

$$\left\langle{\overline{R_{r}}(\mathbf{n})}\right\rangle=R_{r}^{*}(\mathbf{n}),$$

(9)

with $R_{r}^{*}(\mathbf{n})$ as in Eq. (2).

However, $\sigma(t^{\prime})$ ($t^{\prime}\in[t,t+\Delta t]$) will generally be correlated with $\sigma(t)$. Neglecting these means to operate in the adiabatic limit. We would like to retain some of these correlations. In order to compute second moments $\left\langle{\overline{R}_{r}(\mathbf{n})\overline{R}_{s}(\mathbf{n})}\right\rangle$ we first use Eq. (8). This leads to averages of the type $\left\langle{R_{r,\sigma(t_{1})}(\mathbf{n})R_{s,\sigma(t_{2})}(\mathbf{n})}\right\rangle$, where $t_{1}$ and $t_{2}$ are two times in the interval from $t$ to $t+\Delta t$. The second moments can then be expressed in terms of $q_{\sigma\to\sigma^{\prime}}(\cdot)$ as follows

	$\displaystyle\left\langle{\overline{R}_{r}(\mathbf{n})\overline{R}_{s}(\mathbf{n})}\right\rangle=$
	$\displaystyle\frac{1}{\Delta t^{2}}\sum_{\sigma\sigma^{\prime}}\int_{t}^{t+\Delta t}\mathrm{d}t_{1}\int_{t_{1}}^{t_{1}+\Delta t}\mathrm{d}t_{2}\Big{\{}\rho^{*}_{\sigma}q_{\sigma\to\sigma^{\prime}}(t_{2}-t_{1})$
	$\displaystyle\quad\quad\times\big{[}R_{r,\sigma}(\mathbf{n})R_{s,\sigma^{\prime}}(\mathbf{n})+R_{r,\sigma^{\prime}}(\mathbf{n})R_{s,\sigma}(\mathbf{n})\big{]}\Big{\}}.$		(10)

Further details are given in Appendix A. Using Eq. (7) we find

	$\displaystyle\left\langle{\overline{R}_{r}(\mathbf{n})\overline{R}_{s}(\mathbf{n})}\right\rangle-R_{r}^{*}(\mathbf{n})R_{s}^{*}(\mathbf{n})=$
	$\displaystyle\frac{1}{\Delta t^{2}}\sum_{\sigma\sigma^{\prime}}\sum_{\ell=2}^{M}\bigg{\{}\rho^{*}_{\sigma}c_{\ell,\sigma}v_{\ell,\sigma^{\prime}}$
	$\displaystyle\times\big{[}R_{r,\sigma}(\mathbf{n})R_{s,\sigma^{\prime}}(\mathbf{n})+R_{r,\sigma^{\prime}}(\mathbf{n})R_{s,\sigma}(\mathbf{n})\big{]}$
	$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\int_{t}^{t+\Delta t}\mathrm{d}t_{1}\int_{t_{1}}^{t_{1}+\Delta t}\mathrm{d}t_{2}~{}e^{\lambda\mu_{\ell}(t_{2}-t_{1})}\bigg{\}}.$		(11)

For fixed $\ell\in\{2,\dots,M\}$ the integral in the last expression evaluates to

	$\displaystyle\int_{t}^{t+\Delta t}\mathrm{d}t_{1}\int_{t_{1}}^{t_{1}+\Delta t}\mathrm{d}t_{2}~{}e^{\lambda\mu_{\ell}(t_{2}-t_{1})}=$
	$\displaystyle-\dfrac{\Delta t}{\lambda\mu_{\ell}}+\frac{1}{(\lambda\mu_{\ell})^{2}}\left[e^{\lambda\mu_{\ell}\Delta t}-1\right].$		(12)

For $\lambda\Delta t\gg 1$ the first term dominates after inserting into Eq. (11), as also observed in Hufton et al. (2016). We are then left with

	$\displaystyle\left\langle{\overline{R}_{r}(\mathbf{n})\overline{R}_{s}(\mathbf{n})}\right\rangle-R_{r}^{*}(\mathbf{n})R_{s}^{*}(\mathbf{n})=$
	$\displaystyle-\frac{1}{\lambda\Delta t}\sum_{\sigma\sigma^{\prime}}\sum_{\ell=2}^{M}\bigg{\{}\frac{1}{\mu_{\ell}}R_{r,\sigma}(\mathbf{n})R_{s,\sigma^{\prime}}(\mathbf{n})$
	$\displaystyle~{}~{}~{}~{}\times\big{[}\rho^{*}_{\sigma}c_{\ell,\sigma}v_{\ell,\sigma^{\prime}}+\rho^{*}_{\sigma^{\prime}}c_{\ell,\sigma^{\prime}}v_{\ell,\sigma}\big{]}\bigg{\}}.$		(13)

The main challenge in implementing the $\tau$FE algorithm is then to find the average rates from Eq. (9) for all $r$, and the second moments from Eq. (13) for any pair $r,s$ of reactions affected by the environment.

Without loss of generality, we assume that only the rates for the reactions $r=1,\dots,L$ ($L\leq R$) depend on the environmental state $\sigma$.

The $\tau$FE algorithm with time step $\Delta t$ proceeds as follows:

1. 1.

   Initiate the population in state $\mathbf{n}(0)$. Set time to $t=0$.

2. 2.

   Compute $R^{*}_{r}(\mathbf{n})$ for $r=1,\dots,L$ using Eq. (2), and the covariances $\Xi_{rs}(\mathbf{n})\equiv\left\langle{\overline{R}_{r}(\mathbf{n})\overline{R}_{s}(\mathbf{n})}\right\rangle-R_{r}^{*}(\mathbf{n})R_{s}^{*}(\mathbf{n})$ using Eq. (13) for every pair $r,s\in\{1,\dots,L\}$.

3. 3.

   (i) First consider the reactions with rates dependent on the environment: Draw correlated Gaussian random numbers $\ell_{1},\dots,\ell_{L}$ such that $\left\langle{\ell_{r}}\right\rangle=R^{*}_{r}(\mathbf{n})$, and $\left\langle{\ell_{r}\ell_{s}}\right\rangle-\left\langle{\ell_{r}}\right\rangle\left\langle{\ell_{s}}\right\rangle=\Xi_{rs}(\mathbf{n})$. If $\ell_{r}<0$ for any $r\in\{1,\dots,L\}$ set $\ell_{r}=0$.
    
   (ii) For the remaining reactions $r\in\{L+1,\dots,R\}$ set $\ell_{r}=R_{r}(\mathbf{n})$. These are the reactions with rates independent of the environment.

4. 5.

   Draw independent Poissonians random numbers $m_{r}$, $r=1,\dots,R$, each with parameter $\ell_{r}\Delta t$.

5. 6.

   Update the state of the population, $\mathbf{n}(t+\Delta t)=\mathbf{n}(t)+\sum_{r}m_{r}{\mbox{\boldmath$\nu$}}_{r}$.

6. 7.

   Increment time by $\Delta t$ and go to 2.

We note that the mean of the $\ell_{r}$ in step 3(i) is of order $(\lambda\Delta t)^{0}$, and their variance of order $(\lambda\Delta t)^{-1}$. Truncation of the $\ell_{r}$ will therefore only be required very rarely when $\lambda\Delta t\gg 1$.

Evaluating the expressions in Eqs. (2) and (13) in step 2 requires eigenvalues $\mu_{\ell}$ of the transition matrix $\mathbf{A}(\mathbf{n})$ for the environment, the eigenvectors, $\mathbf{v}_{\ell}$ (including the stationary distribution $\mathbf{v}_{1}={\mbox{\boldmath$\rho$}}^{*}$), and the coefficients $c_{\ell,\sigma}$ for all $\sigma$. If the environmental process is independent of the state of the population (the $A_{\sigma\to\sigma^{\prime}}$ are not functions of $\mathbf{n}$), then these quantities do not depend on $\mathbf{n}$, and only need to be calculated once at the beginning.

In Sec. IV we first test the $\tau$FE algorithm on different models with discrete environments. However that the algorithm can also be extended to the case of environmental dynamics with continuous states. This will be discussed in Sec. V.

This system models the dynamics of two genes, which produce two different regulatory proteins: X (a transcription factor) and Y (an inhibitor that titrates $X$ into an inactive complex). Specifically, we use the activator-titration circuit described in Lin and Buchler (2018a). The reactions are as follows:

	$\displaystyle\emptyset\xrightarrow{\Omega\beta_{X}}X,\quad\emptyset\xrightarrow{\Omega\beta_{Y,\sigma}}Y,\quad(\sigma=0,1)$
	$\displaystyle X\xrightarrow{\delta_{X}}\emptyset,\quad Y\xrightarrow{\delta_{Y}}\emptyset,$
	$\displaystyle X+Y\xrightarrow{\alpha/\Omega}\emptyset,$
	$\displaystyle E_{0}\xrightarrow{\lambda n_{X}\kappa_{Y}/\Omega}E_{1},\quad\text{and,}\quad E_{1}\xrightarrow{\lambda\theta_{Y}}E_{0},$		(14)

where the $E_{\sigma}$ denote states of the environment $(\sigma=0,1)$. These two environmental states represent situations in which a transcription factor $X$ is bound to the promoter of gene $Y$ (state $E_{1}$), or no transcription factor is bound ($E_{0}$), respectively. The first two reactions in Eq. (14) describe the production of the two proteins ($X$ and $Y$). The production rates are $\beta_{X}$ and $\beta_{Y,\sigma}$. The former is independent of the environmental state, the latter explicitly depends on $\sigma$ (i.e., on the presence or absence of a bound transcription factor). The reactions in the second line of Eq. (14) describe degradation of $X$ and $Y$, and the reaction in the third line captures titration. The binding and unbinding processes of the transcription factor are described by the reactions in the last line. The parameter $\Omega$ in the reaction rates determines the typical number of particles in the system, for further details see Lin and Buchler (2018a). We write $n_{X}$ for the number of $X$-particles in the system, and similarly $n_{Y}$ is the number of $Y$-particles. One finds $n_{X},n_{Y}={\cal O}(\Omega)$ in the stationary state.

Mathematically, the model consists of two species in the population (with numbers of particles $n_{X},n_{Y}$), and two environmental states, $\sigma=0,1$. We therefore have $S=2,M=2$. Eqs. (9) and (13) can be evaluated explicitly for this case, see also Hufton et al. (2019b). The only process affected by the state of the environment is the production of $Y$, with rate $\beta_{Y,\sigma}$. This rate becomes a (clipped) Gaussian random variable in the $\tau$FE algorithm, with first moment

$$\left\langle{\overline{\beta}_{Y}}\right\rangle=\beta_{Y}^{*}=\dfrac{\theta_{Y}\beta_{Y,0}+n_{X}\kappa_{Y}\beta_{Y,1}/\Omega}{\theta_{Y}+n_{X}\kappa_{Y}/\Omega},$$

(15)

and with variance

	$\displaystyle\sigma_{\beta_{Y}\beta_{Y}}^{2}$	$\displaystyle\equiv\langle\bar{\beta}_{Y}^{2}\rangle-{\beta_{Y}^{*}}^{2}$
		$\displaystyle=\dfrac{2n_{X}\kappa_{Y}\theta_{Y}/\Omega}{\lambda\Delta t(n_{X}\kappa_{Y}/\Omega+\theta_{Y})}\left(\beta_{Y,0}-\beta_{Y,1}\right)^{2}.$		(16)

Further details of the derivation can be found in Appendix B.

Simulation results for this model are shown in Fig. 1. In panels (a) and (b) we illustrate typical sample paths obtained from Gillespie simulations of the full model (population and environment), and from the $\tau$FE algorithm, respectively. We also show the stationary distributions for the quantities $n_{X}+n_{Y}$ and $n_{X}-n_{Y}$ as obtained from both simulation algorithms. The distributions in panels (c) and (d) are for $\lambda=10^{3}$ (i.e., moderately fast environmental dynamics), there are then remaining discrepancies between the $\tau$FE algorithm and simulations of the full model. In panels (e) and (f) the time-scale separation is larger ($\lambda=10^{4}$). The agreement improves as indicated by the Jensen-Shannon divergence (JSD) Fuglede and Topsoe (2004); Lin (1991) given in the figure.

We note at this point that the average CPU time to run a sample path up to time $t=10^{3}$ with parameters as in Fig. 1 (e) and (f) is $2.94$ seconds for the Gillespie algorithm, and $0.03$ seconds for the $\tau$FE algorithm (with a time step $\Delta t=0.1$). These average simulation times are obtained from ten runs. They indicate that the $\tau$FE algorithm can significantly increase efficiency while producing results of the quality shown in Fig. 1. We stress that our primary interest is the relative comparison of computing times, and not on absolute simulation times ¹¹1For completeness, we add that simulations were performed on a MacBook Pro (Mid 2014), with processor 2.6 GHz Dual-Core Inter Core i5, and memory 8 GB 1600 MHz DDR3..

Next, we consider a two-species birth-death process subject to an external environment which can be in one of three different states. This is a toy model chosen for illustration and does not represent any specific natural system. However, it captures elements of models of population dynamics.

The species in the population are labeled $A$ and $B$, and the environmental states $\sigma=0,1,2$. Particles of type $A$ are produced with rate $\Omega\alpha_{\sigma}$, and particles of type $B$ with rate $\Omega\beta_{\sigma}$. The subscript $\sigma$ indicates explicit dependence on the state of the environment. Particles are removed with constant per capita rates $d_{A}$ and $d_{B}$ respectively. The parameter $\Omega$ again sets the typical size of the population. We write $n_{A}$ and $n_{B}$ for the number of individuals of either species. The environmental states cycle stochastically through the sequence $\sigma=0,1,2,0,\dots$, with rate constants $\lambda k_{1}$, $\lambda k_{2}$ and $\lambda k_{0}$ for the three transitions. Mathematically, the reactions in this model are

	$\displaystyle\emptyset\xrightarrow{\Omega\alpha_{\sigma}}A,\quad\emptyset\xrightarrow{\Omega\beta_{\sigma}}B,\quad(\sigma=0,1,2)$
	$\displaystyle A\xrightarrow{d_{A}}\emptyset,\quad B\xrightarrow{d_{B}}\emptyset,$
	$\displaystyle E_{0}\xrightarrow{\lambda k_{1}}E_{1},\quad E_{1}\xrightarrow{\lambda k_{2}}E_{2},\quad E_{2}\xrightarrow{\lambda k_{0}}E_{0},$		(17)

where as before $E_{\sigma}$ denotes the environment. The rates $k_{0},k_{1},$ and $k_{2}$ are constant parameters, independent of the population state.

Details of the calculation of the average rates and their second moments can be found in Appendix C. The average production rates for the two types of particles are

	$\displaystyle\alpha^{*}$	$\displaystyle=$	$\displaystyle\dfrac{k_{0}\alpha_{0}+k_{1}\alpha_{1}+k_{2}\alpha_{2}}{k_{0}+k_{1}+k_{2}},$
	$\displaystyle\beta^{*}$	$\displaystyle=$	$\displaystyle\dfrac{k_{0}\beta_{0}+k_{1}\beta_{1}+k_{2}\beta_{2}}{k_{0}+k_{1}+k_{2}},$		(18)

while the covariance $\sigma_{\alpha\beta}=\sigma_{\beta\alpha}\equiv\langle\bar{\alpha}\bar{\beta}\rangle-\alpha^{*}\beta^{*}$ takes the form

	$\displaystyle\sigma_{\alpha\beta}$	$\displaystyle=\dfrac{\theta^{2}}{\lambda\Delta t}\Big{\{}(\alpha_{0}-\alpha_{1})(\beta_{0}-\beta_{1})\left(3k_{0}^{2}-k_{0,1}k_{0,2}\right)$
		$\displaystyle\quad+(\alpha_{1}-\alpha_{2})(\beta_{1}-\beta_{2})\left(3k_{1}^{2}-k_{1,0}k_{1,2}\right)$
		$\displaystyle\quad+(\alpha_{0}-\alpha_{2})(\beta_{0}-\beta_{2})\left(3k_{2}^{2}-k_{2,0}k_{2,1}\right)\Big{\}},$		(19)

with

$$\theta^{2}\equiv\dfrac{k_{0}k_{1}k_{2}}{(k_{0}k_{1}+k_{1}k_{2}+k_{2}k_{0})^{3}},$$

(20)

and $k_{\sigma,\sigma^{\prime}}=k_{\sigma}-k_{\sigma^{\prime}}$, for $\sigma,\sigma^{\prime}\in\{0,1,2\}$. The variance $\sigma_{\alpha\alpha}\equiv\langle\bar{\alpha}\bar{\alpha}\rangle-(\alpha^{*})^{2}$, is obtained by replacing all instances of $\beta_{\sigma}$ on the right-hand side of Eq. (19) with $\alpha_{\sigma}$. The analog $\sigma_{\beta\beta}$ is obtained similarly by replacing $\alpha_{\sigma}$ with $\beta_{\sigma}$.

In Fig. 2 we report results from numerical simulations for this model, both from conventional Gillespie algorithm of the full systems of environment and population, and using the $\tau$FE algorithm. Panels (a)–(f) show the stationary distributions of $n_{A}+n_{B}$ and $n_{A}-n_{B}$. As seen from the data for example in panel (a) the $\tau$FE algorithm displays deviations from Gillespie simulations of the full model when the environmental process is not sufficiently fast. We quantify these deviations again through the Jensen–Shannon divergence between the two distributions. The deviations reduce as the time scale separation $\lambda$ is increased, i.e., when the environmental process becomes faster relative to the dynamics within the population.

In order to examine if the $\tau$FE algorithm accurately reproduces dynamical features (i.e., properties of the system beyond the stationary distribution), we show spectral densities of the time series for $n_{A}$ and $n_{B}$ in Fig. 2(g) and (h). The spectral densities are defined as

	$\displaystyle S_{AA}(\omega)=\langle|\hat{n}_{A}(\omega)|^{2}\rangle,$
	$\displaystyle S_{AB}(\omega)=\langle\hat{n}_{A}^{\dagger}(\omega)\hat{n}_{B}(\omega)\rangle,$		(21)

where $\hat{n}_{A}(\omega)$ and $\hat{n}_{B}(\omega)$ are the Fourier transforms of $n_{A}(t)$ and $n_{B}(t)$, respectively. The dagger denotes complex conjugation. The data from the $\tau$FE algorithm (open symbols) in Fig. 2(g) and (h) compares well with spectra obtained from direct Gillespie simulations of the full model (solid lines). This shows the $\tau$FE method indeed captures the dynamics of $n_{A}$ and $n_{B}$. We also provide a comparison against the spectral densities obtained from conventional $\tau-$leaping simulations in the adiabatic limit, i.e., simulations with constant rates $\alpha^{*}$ and $\beta^{*}$ for the production events [Eq. (18)]. These are shown as full markers in Fig. 2(g) and (h). One then finds more substantial systematic deviations. This is because environmental fluctuations are discarded in the adiabatic limit. The $\tau$FE algorithm on the other hand captures the stochasticity of the environment to sub-leading order in $\lambda^{-1}$ in each iteration step.

We now consider a model studied in Lin et al. (2018); Hufton et al. (2019a), describing a single gene $G$ with a promoter site which can bind to a total of up to $N$ molecules of protein. The number of protein molecules bound, $\sigma$, plays the role of the environment in this setting. The rate for transitions from $\sigma$ to $\sigma+1$ depends on the number of protein molecules. The reactions in this model can be summarised as follows,

	$\displaystyle G_{\sigma}+P\xrightleftharpoons[\lambda k_{-}]{\lambda k_{+}/\Omega}G_{\sigma+1},$	$\displaystyle\quad\text{for}\quad\sigma<N,$
	$\displaystyle G_{\sigma}\xrightarrow{\Omega b_{\sigma}}G_{\sigma}+M,$
	$\displaystyle M\xrightarrow{d}\emptyset,\quad M\xrightarrow{\beta}M+P,$	$\displaystyle\quad P\xrightarrow{\delta}\emptyset,$		(22)

where $M$ and $P$ refer to molecules of mRNA and protein, respectively. The production rate $b_{\sigma}$ for mRNA depends on the number of protein molecules bound to the promoter. We refer to Lin et al. (2018); Hufton et al. (2019a) for further details. In the following we write $n_{M}$ and $n_{P}$ for the numbers of particles of either type. One interesting feature of this model is that the distribution of the protein and mRNA populations can become bimodal, as illustrated in Fig. 3. This leads to bistability, with trajectories transitioning between the two modes of the joint distribution of $n_{P}$ and $n_{M}$. Hence, the model describes a genetic switch.

In this model only the production rate of mRNA molecules is affected by the state of the environment. The average mRNA-production rate is found as

$$b^{*}=\dfrac{k_{-}^{2}b_{0}+k_{-}\tilde{k}_{+}b_{1}+\tilde{k}_{+}^{2}b_{2}}{k_{-}^{2}+k_{-}\tilde{k}_{+}+\tilde{k}_{+}^{2}},$$

(23)

with $\tilde{k}_{+}=k_{+}n_{P}/\Omega$. The second moment of the production rate takes the form

	$\displaystyle\sigma_{bb}^{2}$	$\displaystyle\equiv\langle\bar{b}^{2}\rangle-{b^{*}}^{2}$
		$\displaystyle=\dfrac{\theta^{2}}{\lambda\Delta t}\Big{\{}(b_{0}-b_{1})^{2}k_{-}\left(k_{-}^{2}+\tilde{k}_{+}k_{-}-\tilde{k}_{+}^{2}\right)$
		$\displaystyle\quad+(b_{0}-b_{2})^{2}2k_{-}\tilde{k}_{+}\left(k_{-}+\tilde{k}_{+}\right)$
		$\displaystyle\quad+(b_{1}-b_{2})^{2}\tilde{k}_{+}\left(\tilde{k}_{+}^{2}+\tilde{k}_{+}k_{-}-k_{-}^{2}\right)\Big{\}},$		(24)

with

$$\theta^{2}=\dfrac{2k_{-}\tilde{k}_{+}}{\left(k_{-}^{2}+k_{-}\tilde{k}_{+}+\tilde{k}_{+}^{2}\right)^{3}}.$$

(25)

Details of the calculation leading to Eqs. (23) and (24) can be found in Appendix D.

Figure 3 shows the stationary joint distribution of the number of mRNA and protein molecules for different values of the time-scale separation parameter $\lambda$. The figure shows data from Gillespie simulations of the full model [panels (a)–(c)], and data from the $\tau$FE algorithm [panels (d)–(f)]. The $\tau$FE algorithm captures the distribution profile with two local maxima. For low values of $\lambda$ (i.e., a relatively slow environmental process) the distribution obtained from $\tau$FE tends to be wider than those from the Gillespie algorithm. The agreement improves for faster environments, as indicated again by the Jensen–Shannon distances in Fig. 3.

In Fig. 4, we show the distribution and means of the sojourn times $t_{\ell}$ and $t_{\rm h}$ near the lower and higher modes of the stationary distribution. More precisely this is the time between entering and leaving a designated region around each of the modes. The lower maximum of the stationary distribution is sharper than the upper maximum (Fig. 3). Accordingly, we have chosen a smaller region at the lower mode than at the upper mode. For the lower mode, we use the region $0\leq n_{M}\leq 20$, $0\leq n_{P}\leq 1100$ which encloses the mode at $(n_{M},n_{P})=(10,500)$. For the higher mode we use the region $20\leq n_{M}\leq 80$, $1100\leq n_{P}\leq 2700$ enclosing the mode at $(n_{M},n_{P})=(30,1800)$.

The data shown in the figure is constructed from one long sample path (run until $t=10^{6}$), recording the points in time at which the system enters or leaves either region. Gillespie simulations operate in continuous time and the $\tau$FE algorithm in discrete time. In order to remove any artefacts resulting from this difference, the same time resolution ($0.05$) is used in both algorithms for the measurement of arrival and departure times. Because the lower mode is sharper than the upper maximum and because the sizes of the two detection regions are different the sojourn time $t_{\ell}$ at the lower mode is found to be smaller than that at the higher mode, $t_{\rm h}$.

The distributions of sojourn times in Figs. 4 (a) and (b) indicate that the $\tau$FE algorithm captures this dynamic quantity, provided the environmental process is sufficiently fast. This is confirmed in panels (c) and (d), where we show the mean sojourn times as a function of the relative speed $\lambda$ of the environment compared to the population dynamics. As seen in both panels, the $\tau$FE algorithm generates accurate measurements of the mean sojourn times $\left\langle{t_{\ell}}\right\rangle$ and $\left\langle{t_{\rm h}}\right\rangle$ in the limit $\lambda\gg 1$.

At the same time, stochastic effects due to the random environmental process are captured for large but finite $\lambda$. This can be seen in Fig. 4 (d): the mean sojourn time $\left\langle{t_{\rm h}}\right\rangle$ drops significantly as the environmental process becomes slower, and hence additional noise is injected into the population (there is no environmental noise in the adiabatic limit). While there are quantitative differences compared to exact simulations, the $\tau$FE algorithm captures this reduction of $\left\langle{t_{h}}\right\rangle$. Panel (c) reveals that there are also limitations to the precision of the $\tau$FE algorithm. The mean sojourn time $\left\langle{t_{\ell}}\right\rangle$ near the lower mode is affected much less by a reduction of the time-scale separation parameter $\lambda$ than the mean sojourn time at the upper mode. This indicates that the escape from this region is driven mostly by intrinsic noise rather than by environmental stochasticity. While the data from the two algorithms remains within approximately $10\%$ for sufficiently fast environmental dynamics ($\lambda\gtrsim 10^{4}$) the $\tau$FE algorithm is unable to capture the small rise of $\left\langle{t_{\ell}}\right\rangle$ observed in Gillespie simulations for intermediate values of $\lambda$.

In Table 1 we compare the the computing time required for both the Gillespie algorithm and the $\tau$FE method for different values of $\lambda$. The data in the table is the CPU time required to generate one sample path up to time $t=10^{3}$, averaged over ten runs. The model parameters are as in Figs. 3 and 4.

The full model comprises the reactions in the population and the environmental switching. The rates for the former reactions are independent of $\lambda$, the rates for the latter scale linearly in $\lambda$. Accordingly, one expects the computing time for Gillespie simulations of the full model to be linear in $\lambda$, with a non-zero intercept. The data in the table is consistent with this. We note that Gillespie algorithm does not require any time discretisation.

The running time for the $\tau$FE algorithm depends on the choice of the time step. The time step in turn affects the accuracy of the outcome. If $\Delta t$ is large, then $\tau$FE simulations are fast, but the approximation to the continuous-time full model becomes less good. On the other hand the time step must not be too small, as the construction of the algorithm requires sufficient averaging of the environmental process in each step [Eqs. (11)–(13)]. The time step for the $\tau$FE algorithm in Figs. 3 and 4, and in Table 1 is chosen inversely proportional to $\lambda$. This is to ensure that each time step captures a sufficient number of switches of the environmental state. Accordingly, we expect the computing time for the $\tau$FE algorithm to scale linearly in $\lambda$, with no intercept. Again, the running times we measured in our simulations are consistent with this expectation. Overall, Table  1 shows that the $\tau$FE algorithm is able to generate data of the accuracy as in Figs. 3 and 4 while reducing the computing effort approximately ten fold compared to full Gillespie simulations.

We turn now to systems which are subject to an environment with continuous states. Specifically, we follow Assaf et al. (2013a) and assume that the environmental state $\sigma$ follows an Ornstein–Uhlenbeck process (see also Roberts et al. (2015); Assaf et al. (2013b)),

$$\frac{\mathrm{d}\sigma}{\mathrm{d}t}=\lambda(m-\sigma)+\sqrt{2\lambda v^{2}}\,\eta(t),$$

(26)

where $\eta(t)$ is Gaussian white noise of unit amplitude, in particular $\left\langle{\eta(t)\eta(t^{\prime})}\right\rangle=\delta(t-t^{\prime})$. The parameter $m$ is the average value of $\sigma$ in the long run, whilst $v$ controls the magnitude of noise. As before, the parameter $\lambda>0$ indicates how quickly the environment changes relative to the dynamics in the population; $\lambda$ is the equivalent of $1/\tau_{c}$ in the notation of Assaf et al. (2013a).

The probability distribution of finding the environment in state $\sigma$ at time $t$, given that was in state $\sigma^{\prime}$ at time $t^{\prime}$, can be obtained from the Fokker–Planck equation for the Ornstein–Uhlenbeck process, and is given by (see e.g. Klebaner (2005); Risken (1996))

	$\displaystyle q_{\sigma^{\prime}\to\sigma}(t-t^{\prime})=\sqrt{\dfrac{1}{2\pi v^{2}(1-e^{-2\lambda(t-t^{\prime})})}}$
	$\displaystyle\times\exp\left[-\frac{\left(\sigma-\sigma^{\prime}e^{-\lambda(t-t^{\prime})}-m\left(1-e^{-\lambda(t-t^{\prime})}\right)\right)^{2}}{2v^{2}\left(1-e^{-2\lambda(t-t^{\prime})}\right)}\right].$		(27)

For $t\to\infty$ (and $t^{\prime}$ fixed) this quantity tends to the stationary distribution

$$\rho^{*}_{\sigma}=\sqrt{\dfrac{1}{2\pi v^{2}}}\exp\left[-\dfrac{\left(\sigma-m\right)^{2}}{2v^{2}}\right].$$

(28)

We note that it is not a requirement for the $\tau$FE algorithm that the environment follows an Ornstein–Uhlenbeck process. However, both functions $q_{\sigma^{\prime}\to\sigma}(t-t^{\prime})$ and $\rho^{*}_{\sigma}$ are required, as discussed in more detail below.

We proceed to describe how the $\tau$FE algorithm can be implemented for models with continuous environments (Sec. V.2).

In the case of discrete environments, continuous-time sample paths of the full model can be generated using the conventional Gillespie algorithm. This is an exact procedure: the ensemble of these sample paths faithfully describes the statistics of the full model. In Sec. IV we have used this as a benchmark to test the $\tau$FE algorithm. We are not aware of any analogous exact simulation method for models of discrete populations in a stochastic environment with continuous states. In order to test the $\tau$FE algorithm we therefore compare outcomes against those from approximation methods to generate paths of the combined set of the population and the environment. Several such methods exist, we describe these in Sec. V.3. The tests of the $\tau$FE algorithm against the baseline of these methods are described in Sec. VI.

We proceed similar to discrete case in Sec. III, replacing the sums over $\sigma$ in Eqs. (2) and (10) with integrals. We then have

$$R_{r}^{*}(\mathbf{n})=\int_{-\infty}^{\infty}\mathrm{d}\sigma\rho^{*}_{\sigma}R_{r,\sigma}(\mathbf{n}),$$

(29)

and the relation for the second moments turns into

	$\displaystyle\left\langle{\overline{R}_{r}(\mathbf{n})\overline{R}_{s}(\mathbf{n})}\right\rangle=$
	$\displaystyle\frac{1}{\Delta t^{2}}\int_{-\infty}^{\infty}\mathrm{d}\sigma\int_{-\infty}^{\infty}\mathrm{d}\sigma^{\prime}\int_{t}^{t+\Delta t}\mathrm{d}t_{1}\int_{t_{1}}^{t+\Delta t}\mathrm{d}t_{2}$
	$\displaystyle\times\Big{\{}\rho^{*}_{\sigma}q_{\sigma\to\sigma^{\prime}}(t_{2}-t_{1})$
	$\displaystyle\big{[}R_{r,\sigma}(\mathbf{n})R_{s,\sigma^{\prime}}(\mathbf{n})+R_{r,\sigma^{\prime}}(\mathbf{n})R_{s,\sigma}(\mathbf{n})\big{]}\Big{\}}.$		(30)

Depending on the form of the stationary distribution $\rho^{*}_{\sigma}$, the kernel $q_{\sigma\to\sigma^{\prime}}(t_{2}-t_{1})$ and the rates $R_{r,\sigma}(\mathbf{n})$ the integrals in Eqs. (29) and (30) can be carried out, and closed-form analytical expressions can be obtained. In Sec. VI we explore a number of different examples, further scenarios are also discussed Appendix F. Once the average rates and the second moments are calculated, the $\tau$FE algorithm is implemented as described in Sec. III.3.

In this section we summarise ‘conventional’ approaches to simulating discrete Markovian systems subject to environmental dynamics with continuous states. By ‘conventional’ we mean methods which produce explicit (approximate) sample paths of the environmental process. This is in contrast to the $\tau$FE algorithm, which generates paths only of the system proper.