Tau leaping

Algorithm

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

x ( t + τ ) = x ( t ) + τ x ′ ( t )

the change is

x ( t + τ ) = x ( t ) + P ( τ x ′ ( t ) )

where P ( τ x ′ ( t ) ) is a Poisson distributed random variable with mean τ x ′ ( t ) .

Given a state x ( t ) = { X i ( t ) } with events E j occurring at rate R j ( x ( t ) ) and with state change vectors v j (where i indexes the state variables, and j indexes the events), the method is as follows:

1. Initialise the model with initial conditions x ( t 0 ) = { X i ( t 0 ) } .
2. Calculate the event rates R j ( x ( t ) ) .
3. Choose a time step τ . This may be fixed, or by some algorithm dependent on the various event rates.
4. For each event E j generate K j ∼ Poisson ( R j τ ) , which is the number of times each event occurs during the time interval [ t , t + τ ) .
5. Update the state by
6. Repeat from Step 2 until some desired condition is met (e.g. a particular state variable reaches 0, or time t 1 is reached).

Algorithm for efficient step size selection

This algorithm is described by Cao et al. The idea is to bound the relative change in each event rate R j by a specified tolerance ϵ (Cao et al. recommend ϵ = 0.03 , although it may depend on model specifics). This is achieved by bounding the relative change in each state variable X i by ϵ / g i , where g i depends on the rate that changes the most for a given change in X i .Typically 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 2 N auxiliary values (where N is the number of state variables X i ), and should only require reusing previously calculated values R j ( x ) . An important factor in this since X i is an integer value, then there is a minimum value by which it can change, preventing the relative change in R j being bounded by 0, which would result in τ also tending to 0.

1. For each state variable X i , calculate the auxiliary values μ i ( x ) = ∑ j v i j R j ( x ) σ i 2 ( x ) = ∑ j v i j 2 R j ( x )
2. For each state variable X i , determine the highest order event in which it is involved, and obtain g i
3. Calculate time step τ as τ = min i { max { ϵ X i / g i , 1 } | μ i ( x ) | , max { ϵ X i / g i , 1 } 2 σ i 2 ( x ) }

This computed τ is then used in Step 3 of the τ leaping algorithm.