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]
Algorithm
The algorithm is analogous to the Euler method for deterministic systems, but instead of making a fixed change
the change is
where is a Poisson distributed random variable with mean
.
Given a state with events
occurring at rate
and with state change vectors
(where
indexes the state variables, and
indexes the events), the method is as follows:
- Initialise the model with initial conditions
.
- Calculate the event rates
.
- Choose a time step
. This may be fixed, or by some algorithm dependent on the various event rates.
- For each event
generate Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.): K_j \sim \text{Poisson}(R_j\tau)
, which is the number of times each event occurs during the time interval .
- Update the state by
- where
is the change on state variable
due to event
. 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
).
- Repeat from Step 2 until some desired condition is met (e.g. a particular state variable reaches 0, or time
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 by a specified tolerance
(Cao et al. recommend
, although it may depend on model specifics). This is achieved by bounding the relative change in each state variable
by
, where
depends on the rate that changes the most for a given change in
.Typically
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 auxiliary values (where
is the number of state variables
), and should only require reusing previously calculated values
. An important factor in this since
is an integer value, then there is a minimum value by which it can change, preventing the relative change in
being bounded by 0, which would result in
also tending to 0.
- For each state variable
, calculate the auxiliary values
- For each state variable
, determine the highest order event in which it is involved, and obtain
- Calculate time step
as
This computed is then used in Step 3 of the
leaping algorithm.
References
- ↑ Lua error in package.lua at line 80: module 'strict' not found.
- ↑ Lua error in package.lua at line 80: module 'strict' not found.
- ↑ Lua error in package.lua at line 80: module 'strict' not found.
- ↑ 4.0 4.1 Lua error in package.lua at line 80: module 'strict' not found.
- ↑ Lua error in package.lua at line 80: module 'strict' not found.
- ↑ Lua error in package.lua at line 80: module 'strict' not found.