Tau-Leaping and Its Step-Size Constraints

Why Approximate Several Events at Once?

Gillespie’s direct method simulates one reaction event at a time, updating the state and propensities after each event. This is exact for a specified continuous-time Markov chain, but can be slow when many events occur in a short interval. Tau-leaping speeds up simulation by advancing a time interval \(\tau\) at once and approximating the number of times each reaction fires during that interval.

The approximation is useful only when propensities do not change too much during a leap. The state update must also remain physically possible: molecule or cell counts cannot be negative. These requirements compete. A longer \(\tau\) can reduce computational cost, but it generally makes both the frozen-propensity approximation and the negative-count risk less acceptable.

From Propensities to Poisson Reaction Counts

Consider a system with \(M\) reaction channels and \(d\) species. Let \(X(t)\) be the state vector, with nonnegative integer components, and let \(a_j(x)\) be the propensity of reaction \(j\) when the state is \(x\). In a short interval of length \(\tau\), the chance of reaction \(j\) firing is approximately \(a_j(x)\tau\), provided its propensity remains nearly constant. The tau-leap approximation makes that assumption for the whole interval.

The counts are independent under the frozen-propensity approximation. They do not generally equal the joint event counts of the original system over an interval in which propensities change as reactions occur. A leap is accurate when the state remains close enough to \(x\) that using \(a_j(x)\) throughout the interval is reasonable.

Here \(\Delta X=X(t+\tau)-x\), \(\mu(x)=\sum_j\nu_j a_j(x)\) is the mean state-change rate, and \(\Gamma(x)=\sum_j\nu_j\nu_j^{\mathsf T}a_j(x)\) is the covariance rate. To obtain these expressions, write \(\Delta X=\sum_j\nu_jK_j\). A Poisson variable with mean \(a_j(x)\tau\) has both mean and variance \(a_j(x)\tau\). Linearity of expectation gives the first formula; independence of the \(K_j\) makes cross-covariances zero and gives the second.

Choosing a Step to Limit Propensity Changes

A leap should be short enough that reaction propensities do not change substantially while counts are being accumulated. One practical local criterion approximates the change in each propensity by a first-order expansion around the current state. For propensity \(a_j\), this gives \(\Delta a_j\approx\nabla a_j(x)^{\mathsf T}\Delta X\). The mean and variance of this linearized change are, respectively, \(\nabla a_j(x)^{\mathsf T}\mu(x)\tau\) and \(\nabla a_j(x)^{\mathsf T}\Gamma(x)\nabla a_j(x)\tau\).

Choose a relative tolerance \(\epsilon\), such as \(0.03\) or \(0.05\), and require both the magnitude of the expected change and its standard deviation to be no more than \(\epsilon a_j(x)\). For positive \(a_j(x)\), the resulting local bounds are:

A zero denominator imposes no bound from that particular condition. Apply both conditions across the propensities being monitored and take the smallest finite bound. This is a local approximation, not a universal guarantee: a nonlinear propensity can change more than its first-order expansion predicts, especially after an unusually large leap. If a propensity is zero at the current state, a relative tolerance based on it is not useful; use an absolute change tolerance or treat the newly enabled reaction with a more cautious step.

Checking the Risk of Negative Populations

A Poisson count is unbounded, so a tau-leap can draw more consuming events than the current population can support. For species \(i\), define the total weighted consumption and production during the leap as \(D_i=\sum_{j:\nu_{ij}<0}(-\nu_{ij})K_j\) and \(Q_i=\sum_{j:\nu_{ij}>0}\nu_{ij}K_j\). The proposed updated count is exactly \(x_i-D_i+Q_i\). Consequently, the actual negative-count event is \(D_i>x_i+Q_i\).

A convenient conservative event is \(D_i>x_i\). It can be checked without knowing the production count, and it contains the actual negative-count event: if \(D_i>x_i+Q_i\), then \(Q_i\geq0\) implies \(D_i>x_i\). Thus \(P(D_i>x_i)\) is an upper bound on the probability of a negative updated count. It can be conservative because production during the same leap may offset consumption. In particular, the condition \(C_i>x_i\), when \(C_i\) denotes consumption, is not the actual negative-count condition if production is possible.

For unit-consumption reactions, the sum of their independent Poisson counts is Poisson, so the conservative tail probability is straightforward to calculate. With stoichiometric coefficients greater than one, \(D_i\) is a weighted sum of Poisson counts and is not generally Poisson; calculate its tail from the specified distribution or use a valid upper bound. If the bound exceeds an acceptable risk threshold, reduce \(\tau\), handle the at-risk reactions more carefully, or switch to exact simulation for the step. A common safeguard is to identify reactions that can rapidly deplete a low-count species and avoid treating their counts as unconstrained Poisson draws over a long leap.

A Practical Step-Selection Sequence

The propensity-change and depletion checks address different failure modes. Neither replaces the other. A disciplined procedure is to apply both at the current state, then re-evaluate after every accepted leap:

Diagnostics, Alternatives, and Limitations

Warning signs include frequent negative-count proposals, large changes in propensities across accepted steps, or results that change materially when the tolerances are tightened. Check that each propensity uses the correct current state, that consumption and production use the appropriate signs of \(\nu_{ij}\), and that reaction counts are drawn with means \(a_j(x)\tau\). Clipping a negative count to zero changes the model’s transition rules and can conceal an unsuitable step size.

The main alternative is Gillespie’s exact direct method, which simulates the next event time and reaction and updates propensities after every event. It avoids the frozen-propensity approximation and Poisson overdraw, but can require many steps when events are frequent. Tau-leaping is attractive when many reactions occur while propensities vary slowly. A hybrid strategy can use leaps for abundant, slowly varying reactions and exact event simulation when low counts or rapid propensity changes make a leap unsafe.

The bounds above are local: the derivative-based propensity criterion is based on a first-order approximation, and the depletion bound may be conservative. Meeting both checks is not a proof that the simulated path has negligible bias. Their tolerances should be chosen in light of the scientific question, and sensitivity to smaller steps should be examined when accuracy matters. Tau-leaping also preserves stochastic variation in reaction counts; it is not simply a deterministic update using expected counts.

Common Mistake