Tau-leaping

In probability theory, tau-leaping, or τ-leaping, is an approximate method for the simulation of a stochastic system.[1] It is based on the Gillespie algorithm, performing all reactions for an interval of length tau before updating the propensity functions.[2] By updating the rates less often this sometimes allows for more efficient simulation and thus the consideration of larger systems.

Many variants of the basic algorithm have been considered.[3][4][5][6][7]

Algorithm

The algorithm is analogous to the Euler method for deterministic systems, but instead of making a fixed change

${\displaystyle x(t+\tau )=x(t)+\tau x'(t)}$

the change is

${\displaystyle x(t+\tau )=x(t)+P(\tau x'(t))}$

where ${\displaystyle P(\tau x'(t))}$ is a Poisson distributed random variable with mean ${\displaystyle \tau x'(t)}$.

Given a state ${\displaystyle \mathbf {x} (t)=\{X_{i}(t)\}}$ with events ${\displaystyle E_{j}}$ occurring at rate ${\displaystyle R_{j}(\mathbf {x} (t))}$ and with state change vectors ${\displaystyle \mathbf {v} _{ij}}$ (where ${\displaystyle i}$ indexes the state variables, and ${\displaystyle j}$ indexes the events), the method is as follows:

1. Initialise the model with initial conditions ${\displaystyle \mathbf {x} (t_{0})=\{X_{i}(t_{0})\}}$.
2. Calculate the event rates ${\displaystyle R_{j}(\mathbf {x} (t))}$.
3. Choose a time step ${\displaystyle \tau }$. This may be fixed, or by some algorithm dependent on the various event rates.
4. For each event ${\displaystyle E_{j}}$ generate ${\displaystyle K_{j}\sim {\text{Poisson}}(R_{j}\tau )}$, which is the number of times each event occurs during the time interval ${\displaystyle [t,t+\tau )}$.
5. Update the state by

   ${\displaystyle \mathbf {x} (t+\tau )=\mathbf {x} (t)+\sum _{j}K_{j}v_{ij}}$

   where ${\displaystyle v_{ij}}$ is the change on state variable ${\displaystyle X_{i}}$ due to event ${\displaystyle E_{j}}$. At this point it may be necessary to check that no populations have reached unrealistic values (such as a population becoming negative due to the unbounded nature of the Poisson variable ${\displaystyle K_{j}}$).

6. Repeat from Step 2 onwards until some desired condition is met (e.g. a particular state variable reaches 0, or time ${\displaystyle t_{1}}$ is reached).

Algorithm for efficient step size selection

This algorithm is described by Cao et al.[4] The idea is to bound the relative change in each event rate ${\displaystyle R_{j}}$ by a specified tolerance ${\displaystyle \epsilon }$ (Cao et al. recommend ${\displaystyle \epsilon =0.03}$, although it may depend on model specifics). This is achieved by bounding the relative change in each state variable ${\displaystyle X_{i}}$ by ${\displaystyle \epsilon /g_{i}}$, where ${\displaystyle g_{i}}$ depends on the rate that changes the most for a given change in ${\displaystyle X_{i}}$. Typically ${\displaystyle g_{i}}$ is equal the highest order event rate, but this may be more complex in different situations (especially epidemiological models with non-linear event rates).

This algorithm typically requires computing ${\displaystyle 2N}$ auxiliary values (where ${\displaystyle N}$ is the number of state variables ${\displaystyle X_{i}}$), and should only require reusing previously calculated values ${\displaystyle R_{j}(\mathbf {x} )}$. An important factor in this is that since ${\displaystyle X_{i}}$ is an integer value, there is a minimum value by which it can change, preventing the relative change in ${\displaystyle R_{j}}$ being bounded by 0, which would result in ${\displaystyle \tau }$ also tending to 0.

1. For each state variable ${\displaystyle X_{i}}$, calculate the auxiliary values

   ${\displaystyle \mu _{i}(\mathbf {x} )=\sum _{j}v_{ij}R_{j}(\mathbf {x} )}$

   ${\displaystyle \sigma _{i}^{2}(\mathbf {x} )=\sum _{j}v_{ij}^{2}R_{j}(\mathbf {x} )}$

2. For each state variable ${\displaystyle X_{i}}$, determine the highest order event in which it is involved, and obtain ${\displaystyle g_{i}}$
3. Calculate time step ${\displaystyle \tau }$ as

   ${\displaystyle \tau =\min _{i}{\left\{{\frac {\max {\{\epsilon X_{i}/g_{i},1\}}}{|\mu _{i}(\mathbf {x} )|}},{\frac {\max {\{\epsilon X_{i}/g_{i},1\}}^{2}}{\sigma _{i}^{2}(\mathbf {x} )}}\right\}}}$

This computed ${\displaystyle \tau }$ is then used in Step 3 of the ${\displaystyle \tau }$ leaping algorithm.

References

1. Gillespie, D. T. (2001). "Approximate accelerated stochastic simulation of chemically reacting systems" (PDF). The Journal of Chemical Physics. 115 (4): 1716–1733. Bibcode:2001JChPh.115.1716G. doi:10.1063/1.1378322.
2. Erhard, F.; Friedel, C. C.; Zimmer, R. (2010). "FERN – Stochastic Simulation and Evaluation of Reaction Networks". Systems Biology for Signaling Networks. p. 751. doi:10.1007/978-1-4419-5797-9_30. ISBN 978-1-4419-5796-2.
3. Cao, Y.; Gillespie, D. T.; Petzold, L. R. (2005). "Avoiding negative populations in explicit Poisson tau-leaping". The Journal of Chemical Physics. 123 (5): 054104. Bibcode:2005JChPh.123e4104C. CiteSeerX 10.1.1.123.3650. doi:10.1063/1.1992473. PMID 16108628. S2CID 1652735.
4. Cao, Y.; Gillespie, D. T.; Petzold, L. R. (2006). "Efficient step size selection for the tau-leaping simulation method" (PDF). The Journal of Chemical Physics. 124 (4): 044109. Bibcode:2006JChPh.124d4109C. doi:10.1063/1.2159468. PMID 16460151.
5. Anderson, David F. (2008-02-07). "Incorporating postleap checks in tau-leaping". The Journal of Chemical Physics. 128 (5): 054103. arXiv:0708.0377. Bibcode:2008JChPh.128e4103A. doi:10.1063/1.2819665. ISSN 0021-9606. PMID 18266441. S2CID 1166923.
6. Chatterjee, Abhijit; Vlachos, Dionisios G.; Katsoulakis, Markos A. (2005-01-08). "Binomial distribution based τ-leap accelerated stochastic simulation". The Journal of Chemical Physics. 122 (2): 024112. Bibcode:2005JChPh.122b4112C. doi:10.1063/1.1833357. ISSN 0021-9606. PMID 15638577.
7. Moraes, Alvaro; Tempone, Raul; Vilanova, Pedro (2014-04-24). "Hybrid Chernoff Tau-Leap". Multiscale Modeling & Simulation. 12 (2): 581–615. CiteSeerX 10.1.1.756.9799. doi:10.1137/130925657. ISSN 1540-3467.