# 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 allows for more efficient simulation and thus the consideration of larger systems.

Cao et al. improved the method to prevent the generation of negative populations.[3][4]

## Algorithm

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

the change is

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} _{j}}$ (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 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.[5] 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 been calculated values ${\displaystyle R_{j}(\mathbf {x} )}$. An important factor in this since ${\displaystyle X_{i}}$ is an integer value, then 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.