Open Access Open Badges Methodology article

A higher-order numerical framework for stochastic simulation of chemical reaction systems

Tamás Székely1*, Kevin Burrage12, Radek Erban3 and Konstantinos C Zygalakis45

Author Affiliations

1 Department of Computer Science, University of Oxford, Oxford, OX1 3QD, UK

2 Department of Mathematics, Queensland University of Technology, Brisbane, Qld 4001, Australia

3 Mathematical Institute, University of Oxford, Oxford, OX1 3LB, UK

4 School of Mathematics, University of Southampton, Southampton, SO17 1BJ, UK

5 Mathematics Section, École Polytechnique Fédérale de Lausanne, Station 8, CH-1015 Lausanne, Switzerland

For all author emails, please log on.

BMC Systems Biology 2012, 6:85  doi:10.1186/1752-0509-6-85

Published: 15 July 2012



In this paper, we present a framework for improving the accuracy of fixed-step methods for Monte Carlo simulation of discrete stochastic chemical kinetics. Stochasticity is ubiquitous in many areas of cell biology, for example in gene regulation, biochemical cascades and cell-cell interaction. However most discrete stochastic simulation techniques are slow. We apply Richardson extrapolation to the moments of three fixed-step methods, the Euler, midpoint and θ-trapezoidal τ-leap methods, to demonstrate the power of stochastic extrapolation. The extrapolation framework can increase the order of convergence of any fixed-step discrete stochastic solver and is very easy to implement; the only condition for its use is knowledge of the appropriate terms of the global error expansion of the solver in terms of its stepsize. In practical terms, a higher-order method with a larger stepsize can achieve the same level of accuracy as a lower-order method with a smaller one, potentially reducing the computational time of the system.


By obtaining a global error expansion for a general weak first-order method, we prove that extrapolation can increase the weak order of convergence for the moments of the Euler and the midpoint τ-leap methods, from one to two. This is supported by numerical simulations of several chemical systems of biological importance using the Euler, midpoint and θ-trapezoidal τ-leap methods. In almost all cases, extrapolation results in an improvement of accuracy. As in the case of ordinary and stochastic differential equations, extrapolation can be repeated to obtain even higher-order approximations.


Extrapolation is a general framework for increasing the order of accuracy of any fixed-step stochastic solver. This enables the simulation of complicated systems in less time, allowing for more realistic biochemical problems to be solved.

Stochastic simulation algorithms; τ-leap; High-order methods; Monte Carlo error