Tau-Leaping as an Approximation to Gillespie Simulation

From Event-by-Event Simulation to Tau-Leaping

In Gillespie’s direct method, a simulation generates one reaction event at a time. After each event, it updates the state and recalculates the reaction propensities. This is exact for the specified continuous-time Markov chain, but can be computationally expensive when many reactions occur in a short interval. Tau-leaping speeds up the simulation by advancing a time step \(\tau\) at once: it approximates the number of times each reaction fires during that interval by a Poisson random variable, using propensities evaluated at the interval’s initial state.

Let \((\Omega,\mathcal{F},\mathbb{P})\) be a probability space, and let \((X_t)_{t\ge0}\) be a nonexplosive continuous-time Markov chain with state space \(S\subseteq\mathbb{Z}_{\ge0}^d\). There are \(R\) reaction channels. Reaction \(r\) changes the state by a stoichiometric vector \(\nu_r\in\mathbb{Z}^d\), and has propensity \(a_r(x)\ge0\) in state \(x\): conditional on the current state \(x\), the probability that reaction \(r\) fires in a short interval of length \(h\) is \(a_r(x)h+o(h)\). Write \(a_0(x)=\sum_{r=1}^R a_r(x)\). The process is considered with its natural filtration \((\mathcal{F}_t)_{t\ge0}\), generated by its history up to time \(t\).

The approximation treats reaction counts as if their rates remained fixed over the interval. In the actual chain, every reaction can change the state and therefore change all the propensities; Gillespie simulation accounts for those changes after each event. Tau-leaping is most plausible when a step is short enough that the propensities do not change much, and when the reactions expected to occur are numerous enough for batching to save work.

What the Frozen Poisson Counts Imply

For a Poisson random variable with mean \(m\), both its expectation and variance equal \(m\). Conditional independence therefore gives the following useful calculations for a tau-leap.

Proof. Write \(\Delta X=\sum_r\nu_rK_r\). Since \(\mathbb{E}[K_r\mid\mathcal{F}_t]=a_r(x)\tau\), linearity of conditional expectation gives \(\mathbb{E}[\Delta X\mid\mathcal{F}_t]=\sum_r\nu_ra_r(x)\tau\). For the covariance, independence conditional on \(\mathcal{F}_t\) makes every cross-covariance between distinct counts zero, and \(\operatorname{Var}(K_r\mid\mathcal{F}_t)=a_r(x)\tau\). Thus \(\operatorname{Cov}(\Delta X\mid\mathcal{F}_t)=\sum_r\nu_r\nu_r^{\mathsf T}a_r(x)\tau\). If a propensity is zero, its count is identically zero and contributes zero to both formulas, so the result also covers that case. \(\square\)

Why the Approximation Is Locally Accurate

A small-step comparison explains what tau-leaping captures and what it neglects. The next result concerns reaction counts within a single step, not a claim that a large-step leap reproduces the full trajectory distribution. Assume the propensities are finite at the starting state and bounded by some finite \(\Lambda\) in the states considered during a sufficiently short interval. This bound ensures that the probability of two or more events is of order \(\tau^2\).

Proof. For the exact chain, before its first event the state remains \(x\), so the probability of no event by time \(\tau\) is \(e^{-a_0(x)\tau}=1-a_0(x)\tau+O(\tau^2)\). The density that the first event occurs at time \(s\) and is of type \(r\) is \(e^{-a_0(x)s}a_r(x)\). If that event occurs, the probability of a further event before time \(\tau\) is at most \(\Lambda(\tau-s)\): a rate bounded by \(\Lambda\) gives probability at most \(\Lambda u\) of an event in an interval of length \(u\). Consequently, the probability of exactly one event of type \(r\) is \(a_r(x)\tau+O(\tau^2)\): replacing \(e^{-a_0(x)s}\) by \(1\) in the first-event integral makes an error of order \(\tau^2\), and the chance of a further event contributes at most \(a_r(x)\Lambda\int_0^\tau(\tau-s)\,ds= a_r(x)\Lambda\tau^2/2\). The probability of at least two events is \(O(\tau^2)\), since after a first event the rate is bounded by \(\Lambda\), and integrating the bound on a second event over the possible first-event times gives at most \(\Lambda^2\tau^2/2\).

For the tau-leap, the probability all counts are zero is \(e^{-a_0(x)\tau}\), since the counts are independent and their means sum to \(a_0(x)\tau\). The probability that \(K_r=1\) and all other counts are zero is \(a_r(x)\tau e^{-a_0(x)\tau}=a_r(x)\tau+O(\tau^2)\). Finally, two or more total firings have probability \(O(\tau^2)\): the sum of independent Poisson counts is Poisson with mean \(a_0(x)\tau\), and its probability of at least two is \(1-e^{-a_0(x)\tau}(1+a_0(x)\tau)=O(\tau^2)\). These calculations prove the stated first-order agreement. They do not remove the error caused by freezing rates over a step of non-negligible size. \(\square\)

Step Size, Changing Rates, and Invalid States

A smaller \(\tau\) generally reduces the time during which frozen propensities can diverge from the actual rates. But it also requires more updates to cover a fixed time interval, reducing the computational advantage. A larger \(\tau\) batches more reactions, but can miss important changes in rates and can produce an invalid state. There is no universal step size that guarantees accuracy for every reaction network.

Common pitfall: clipping negative counts. A tempting correction is to replace a Poisson count \(K\) by \(\min(K,2)\) in the pair-reaction example, so that \(A=4-2\min(K,2)\) cannot be negative. At \(\tau=1/2\), the original count has mean \(1\), but the clipped count has expectation \(\mathbb{E}[\min(K,2)]=\mathbb{P}(K\ge1)+\mathbb{P}(K\ge2)=(1-e^{-1})+(1-2e^{-1})=2-3/e\), not \(1\). Clipping therefore changes the count distribution and its mean; it is a different approximation, not a harmless implementation of the stated tau-leap. Better responses include reducing the step, using a method that respects population constraints, or using event-by-event simulation when the state is too close to a boundary for a leap.

In practice, step-size choices should consider both rate changes and the sizes of the populations involved. When a reaction is expected to change a critical, low-count species by a substantial fraction of its available population, a large leap is risky. When propensities vary substantially with state, freezing them for a long interval is also risky. Tau-leaping is valuable when batching meaningfully reduces computation while the step remains short relative to these changes; its Poisson counts are a modeling approximation, not a guarantee of valid states or exact trajectories.