Tau-leaping

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).

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.