Efficient simulation of stochastic chemical kinetics with the Stochastic Bulirsch-Stoer extrapolation method

Tamás Székely, Kevin Burrage, Konstantinos C Zygalakis, Manuel Barrio
2014 BMC Systems Biology  
Biochemical systems with relatively low numbers of components must be simulated stochastically in order to capture their inherent noise. Although there has recently been considerable work on discrete stochastic solvers, there is still a need for numerical methods that are both fast and accurate. The Bulirsch-Stoer method is an established method for solving ordinary differential equations that possesses both of these qualities. Results: In this paper, we present the Stochastic Bulirsch-Stoer
more » ... hod, a new numerical method for simulating discrete chemical reaction systems, inspired by its deterministic counterpart. It is able to achieve an excellent efficiency due to the fact that it is based on an approach with high deterministic order, allowing for larger stepsizes and leading to fast simulations. We compare it to the Euler τ -leap, as well as two more recent τ -leap methods, on a number of example problems, and find that as well as being very accurate, our method is the most robust, in terms of efficiency, of all the methods considered in this paper. The problems it is most suited for are those with increased populations that would be too slow to simulate using Gillespie's stochastic simulation algorithm. For such problems, it is likely to achieve higher weak order in the moments. Conclusions: The Stochastic Bulirsch-Stoer method is a novel stochastic solver that can be used for fast and accurate simulations. Crucially, compared to other similar methods, it better retains its high accuracy when the timesteps are increased. Thus the Stochastic Bulirsch-Stoer method is both computationally efficient and robust. These are key properties for any stochastic numerical method, as they must typically run many thousands of simulations. Spain Full list of author information is available at the end of the article of a system, their results are not always representative [7, 8] . A common stochastic modelling approach is to consider the system as a continuous-time Markov jump process [9] . The stochastic simulation algorithm (SSA) of Gillespie [10] is a simple and exact method for generating Markov paths. However, because it keeps track of each reaction, it can be too computationally costly for more complex systems or those with frequent reactions. Many approximate methods have since been developed, which use similar principles as the SSA but group many reactions into a single calculation, reducing computational time (for a recent review, see [11] ). The first of these is commonly called the Euler or Poisson τ -leap [12] ; it corresponds to the Euler method for ordinary differential equations (ODEs), and samples
doi:10.1186/1752-0509-8-71 pmid:24939084 pmcid:PMC4085235 fatcat:jpp7spfciradnkmenx2fbtvj2m