Simulation methods with extended stability for stiff biochemical Kinetics
© Rué et al. 2010
Received: 19 November 2009
Accepted: 11 August 2010
Published: 11 August 2010
Skip to main content
© Rué et al. 2010
Received: 19 November 2009
Accepted: 11 August 2010
Published: 11 August 2010
With increasing computer power, simulating the dynamics of complex systems in chemistry and biology is becoming increasingly routine. The modelling of individual reactions in (bio)chemical systems involves a large number of random events that can be simulated by the stochastic simulation algorithm (SSA). The key quantity is the step size, or waiting time, τ, whose value inversely depends on the size of the propensities of the different channel reactions and which needs to be re-evaluated after every firing event. Such a discrete event simulation may be extremely expensive, in particular for stiff systems where τ can be very short due to the fast kinetics of some of the channel reactions. Several alternative methods have been put forward to increase the integration step size. The so-called τ-leap approach takes a larger step size by allowing all the reactions to fire, from a Poisson or Binomial distribution, within that step. Although the expected value for the different species in the reactive system is maintained with respect to more precise methods, the variance at steady state can suffer from large errors as τ grows.
In this paper we extend Poisson τ-leap methods to a general class of Runge-Kutta (RK) τ-leap methods. We show that with the proper selection of the coefficients, the variance of the extended τ-leap can be well-behaved, leading to significantly larger step sizes.
The benefit of adapting the extended method to the use of RK frameworks is clear in terms of speed of calculation, as the number of evaluations of the Poisson distribution is still one set per time step, as in the original τ-leap method. The approach paves the way to explore new multiscale methods to simulate (bio)chemical systems.
It is by now very well known that the biochemical kinetics involving small numbers of molecules can be very different to kinetics described by the law of mass action and differential equations [1–3]. This effect is a property of the intrinsic noise of the system and is associated with the uncertainty of knowing when a reaction occurs and what that reaction is. At the molecular level such intrinsic uncertainty is, in turn, a consequence of the stochastic nature of the fluctuations of the potential energy surface for any chemical reaction in the condensed phase . When considering a collection of molecules, the intrinsic noise is accentuated when some chemical species have small numbers, as is often the case in genetic regulatory models where there are small numbers of key transcription factors that can bind to a limited number of operator regions on DNA [5–15]. Kurtz  and Gillespie  realised this fact and developed discrete methods to deal with this situation. The stochastic simulation algorithm (SSA, see  for a review) describes the time evolution of the dynamics of the species in a well-stirred chemically reacting system as a discrete nonlinear Markov process, resulting in an exact method to sample from the probability density function described by the chemical master equation (CME). Gibson and Bruck proposed a more efficient implementation of the SSA called the next reaction method .
The basic idea of the SSA is that at each time point a waiting time to the next reaction and the most likely reaction to occur must be sampled from a joint probability density function leading to an appropriate update of the state vector. But if the rate constants and/or the numbers of molecules in the system are large then the waiting time (time step, τ) can be very small . Because of this Gillespie  introduced the Poisson τ-leap method, in which all reactions are allowed to fire in a given τ with a frequency extracted from a Poisson distribution. Since then many extensions of this idea have been developed. Cao et al.  have considered efficient mechanisms for selecting τ and have developed implicit methods suitable for simulating stiff systems. Tian and Burrage  introduced a modification of Poisson τ-leap methods known as Binomial τ-leap methods that avoids the issue of obtaining negative molecular numbers from which Poisson τ-leap methods can suffer. Chatterjee et al.  and Auger et al.  have considered modifications to Binomial τ-leap methods that improve some of the implementation aspects. On the other hand, Monk  and Mackey  noted the importance of representing delays, especially when representing processes such as transcription and translation. Accordingly, Bratsun et al.  and Barrio et al.  developed a delayed version of the Stochastic Simulation Algorithm. Leier et al.  and Anderson  extended these ideas to a τ-leap setting.
Although τ-leap methods can, in some cases, substantially improve computational efficiency compared with the SSA, when there is moderate stiffness in the system the efficiencies can be quite poor. One could resort to implicit τ-leap methods but then there are considerable implementation issues and subtleties. A different approach is to explore ideas from the numerical ODE (ordinary differential equations) and numerical SDE (stochastic differential equations) communities. Thus, with ODEs it is well known that stiffness leads to a step size restriction when using explicit methods and many classes of efficient implicit methods have been designed . However, in the case of moderately stiff systems explicit Runge-Kutta methods with extended stability regions along the negative real axis have proven to be especially effective [31, 32]. Runge-Kutta methods are a class of one step methods which gain their efficacy by computing intermediate approximations to the solution within a step. Explicit Runge-Kutta methods with extended stability regions are based on explicit Runge-Kutta methods whose stability function is a shifted and scaled Chebyshev polynomial or some variant thereof. In the stochastic setting, there are some subtleties designing fully implicit methods due to possible unboundedness of the solution as the Wiener increment can take positive or negative values with equal likelihood . Thus most methods are semi-implicit, that is implicit in the deterministic component. Abdulle and Cirilli  have, with some success, extended the ideas of explicit Chebyshev methods with extended stability regions to the SDE setting via their class of S-ROCK methods.
Here, we use the Runge-Kutta formulation to construct methods with large stability regions so that efficiencies are gained by allowing larger stepsizes. We note that this is exactly what Abdulle and Cirilli  do in the SDE setting, that is they use a Runge-Kutta formulation to construct methods with excellent stability properties and even though these methods are only weak order 1 they perform very well. It is noteworthy that in this work we are not using the Runge-Kutta formulation to get second order accuracy for τ-leap methods. This seems to be a difficult problem, just as it is the case for SDEs and will probably require double integrals of compensated processes to be simulated. In fact, Abdulle and Cirilli  also note that it is very difficult to construct weak order 2 methods with good stability properties and to our knowledge at the moment no such methods exist in the SDE setting. Note that in a stochastic setting we judge order of accuracy through two mechanisms: strong order (where trajectories are compared with the true solutions) and weak order (where moments are compared). Often a numerical method may have a higher weak order than its strong order. The Euler-Maruyama method is a case in point with weak order one and strong order a half.
Thus, in this paper, we explore a series of fully explicit multistage Runge-Kutta methods with extended stability for a fixed τ-leap stochastic simulation schema. Our methods involve the same number of Poisson evaluations per integration step as in the original τ-leap formulation but allow increasingly larger step sizes at the cost of an increasing series of deterministic evaluations in the internal stages. First we give some background on Runge-Kutta methods for ODEs and SDEs. In section Results we extend these ideas to the τ-leap methods and present a stability analysis for linear chemical kinetics, including its practical implementation. In section Numerical results we present numerical results for both the linear case and the classical stiff system described by the Schlögl reaction . Finally, in section Discussion we discuss further implications of this work and, in particular, possible extensions to multiscale modelling.
where b T = (β 1,...,β s ), w = Ae and e = (1,...,1)T. Here A is the matrix with entries α ij and w is the column vector w T = (w 1,...,w s)T. A Runge-Kutta method is said to be explicit if the s × s matrix A is strictly lower triangular. The method parameters are usually chosen so that a Runge-Kutta method has appropriate efficiency, order and stability characteristics. The Y i are considered to be approximations to the solution at the intermediate points t n + w i h for i = 1,...s.
where e is the unit vector, I s is the identity matrix of order s and ⊗ represents the Kronecker tensor product such that the (i, j) element of A ⊗ B is a ij B. Notice that, if Λ is a scalar value and taking z = h Λ, R(z) would be a scalar and take the form (4). Therefore we can refer to R seamlessly irrespective of whether the argument is a matrix or a scalar.
where r j = b T A j-1 e, j = 1,...,s. Hence, R(z) is a polynomial of at most degree s for any explicit method.
and non-overlapping Wiener increments are independent of one another. A sample of a Wiener increment W(t + h) - W(t) is simulated from a Normal random variable with mean 0 and variance h, N(0, h).
where ΔW n = (ΔW 1,....ΔW d )T and ΔW i := W i (t n + h) - W i (t n ), i = 1,...,d are normally distributed random numbers with mean 0 and variance h. This method is known as the Euler-Maruyama method and it is known to have strong order (pathwise order) and weak order (moment order) 1.
In the later case, the solution is mean square stable if Re [a] + Re [b 2] ≤ 0.
and in the (p, q) plane, with p, q ∈ ℝ, the stability region is a circle of radius 1 centered in (-1,0).
then X(t) = (A(t), B(t), C(t))T, ν 1 = (-1, -1, 1)T and a 1(X(t)) = cA(t)B(t).
and the algorithm repeats.
where L(τ, x) is given by (13) and represents the drift or expected stepchange. As our focus is explicit methods, the matrix A is strictly lower diagonal. We note that (17) requires the same number of samples of Poisson random variables per step as the Poisson τ-leap method.
Indeed any Runge-Kutta method for solving an ODE can be incorporated into this framework. We also note that other methods proposed in the literature can be put into this framework. For example, the midpoint method of Gillespie  can be represented with s = 2, b T = (0, 1), w = (0, 0.5)T and where the row-wise entries of A are 0, 0, 0.5, 0.
If the Runge-Kutta method for ODEs underlying a Runge-Kutta τ-leap method (17) has stability function given by (4), then when the latter is applied to (18) we show (Additional file 1) that
Note that with the constraint |R(z)| < 1, z = -τ(k 1 + k 2) then the spectral radius of R(τ W) is less than or equal to 1, and as there is only one eigenvalue equal to one hence we have boundedness of the mean.
We call this the relative variance at the stationary state associated to R.
Let us consider some particular cases of this result:
Poissonτ-leap For this method R(z) = 1 + z and . Thus, the equilibrium variance doubles at z = -1, it rises fourfold at z = -1.5 and is unbounded at z = -2.
This is the methodology we propose in the following section for the derivation of particular Runge-Kutta methods with s steps. Note that if we required the same limitation on the variance with the standard Poisson τ-leap method we could only take . Thus with the two stage method we can take a stepsize almost six times as large.
This was first shown by Cao et al. . In fact only those Runge-Kutta methods that have a stability function given by (25) can preserve the variance exactly for linear problems. These methods include the implicit midpoint and trapezoidal rules and have to be implicit.
and optimise the value of l s, ϵ such that the range for which this holds is (-l s, ϵ, 0].
Thus, similar to the case s = 2 in which we had one free parameter, γ, to optimise, if we assume r 1 = b T e = 1 then we have s - 1 parameters, r 2,...,r s , we can optimise. In this case, though, the search of the optimal set of parameters has to be performed with numerical optimisation methods rather than analytically. The problem of finding optimal sets of parameters can be stated as a nonlinear program, NLP, and thus its solution approximated numerically (see details in Additional file 1).
Stability regions for methods with bounded relative variance and optimal stability
Stability (l s, ϵ)
Factor vs. τ-leap
Norm. factor τ-leap
It is thus clear that these methods are computationally more efficient than the general case as they only require s-1 evaluations of the expected step-change f(·) instead of the s(s - 1)/2 required in the general framework (17). A collection of methods have been implemented in a branch of the ByoDyn package, v.5.0 .
Details of the Schlögl reaction system
k1 = 3·10-7
k1x(x - 1)A/2
k2 = 10-4
k2x(x - 1)(x - 2)/6
k3 = 3.5
k 3 x
k4 = 10-3
k 4 B
Finally, we have tested the performance of our methods on a larger system of chemical reactions with stiffness due to different reaction time scales and species amounts ranging over several orders of magnitude. For this purpose we considered the Huang and Ferrell model for the mitogen-activated protein kinase (MAPK) cascade . This model is available from the BioModels database  and consists of 22 species interacting through 30 reaction channels. The set of parameters used here (see Additional file 1 for details) renders the model stiff and with species amounts ranging from none up to 3·105 molecules. With the chosen initial conditions the system undergoes a transient change and finally settles down into a stationary state at around t = 150 minutes. We have simulated the model using SSA (Gillespie's Direct Method), the Poisson τ-leap and the RK methods presented here. To produce fair comparisons, all methods have been rewritten in ANSI C using the Mersenne twister  pseudorandom number generator from the GNU Scientific Library. The GNU C Compiler was used to compile the sources with the -O2 optimisation flag. The algorithms were run on an Intel(R) Core(TM)2 Duo Processor E8500 at 3.16 GHz and 6 MB cache. We have run the system to a final time T = 200. Simulations run with SSA took 61, 841 ± 74 seconds.
We have compared the methods in two distinct situations. First we have run them with the same time step τ = 5·10-5. In this case, the Poisson τ-leap method took 51.7 ± 0.4 seconds while Optimal RK τ-leap methods with s = 3 and s = 5 took 86.1 ± 0.4 seconds and 113.9 ± 0.3 seconds respectively. Hence, at the same time step the RK methods are approximately 66% and 120% slower than the Poisson τ-leap due to the multiple evaluations of the propensity functions per step. However, there is an important difference in the results. The relative variance at the steady state is 1.3 (see Additional file 1) for the Poisson τ-leap while for both RK τ-leap methods with s = 3 and s = 5 (ϵ = 0.1) it is less than 1.04.
Then we have compared these methods when run at their respective maximum time steps such that the relative variance at the stationary state is bounded to 1.1 (estimated from the simulations). The maximum time steps allowed with this constraint were: τ = 2·10-5 for the Poisson τ-leap, τ = 3.5·10-4 for the RK τ-leap (3, 0.1) and τ = 9.5·10-4 for the Optimal RK τ-leap (5, 0.1). With this setting, the runtimes obtained were: 111.9 ± 0.7 seconds for the Poisson τ-leap, 15.7 ± 0.06 seconds for the RK τ-leap (3, 0.1) and 7.8 ± 0.02 seconds for the Optimal RK τ-leap (5, 0.1). Thus, in this case the Poisson τ-leap approximately 7.1 and 14.3 times slower than the RK methods, respectively.
Biochemical kinetics typically deals with multiscale problems, in which several scales of time, space and concentrations, simultaneously affect the dynamical behaviour of the system. Thus, the systems biology community is deeply interested in the development of methods that lead to a multiscale view of biochemical systems. As a first step in this workflow, we have presented here a new set of methods that considerably expands the classical τ-leap implementation, from a stability perspective. The importance of the results shown here embraces not only the increase in computational speed for stochastic simulations, a key element for the understanding of the intrinsically noisy biological systems, but more importantly, a way to deal with fast reactions in multiscale settings. The methods developed here have been demonstrated for a first example of stiff system, the classical Schlögl autocatalytic reaction, and can be straightforwardly incorporated into hybrid SSASDE-ODE frameworks.
We see from Table 1 that if we require a bound on the equilibrium variance of 0.1 then the Poisson τ-leap method must take while for the RK methods the bounds on |z| are approximately 4 and 10, respectively with s = 3, 5. This is a very considerable improvement and all the more striking given that the same number of Poisson random variables are simulated per step in all cases.
Initially we had hoped that an approach via Chebyshev methods using ideas from ODEs and SDEs applied to the discrete cases would have been fruitful. It turns out that while such methods have good mean behaviour, the variance behaviour is poor. This is because the variance growth function satisfies (24) and an s-stage Chebyshev method would have s - 1 poles and zeros due to the oscillations in the stability function. Similar issues arise even in the damped forms of the Chebyshev formulation. This means that our optimisation approach is the only way of getting good bounds on ψ(z).
Our results on the nonlinear bimodal Schlögl problem show that the RK methods still behave appropriately even on nonlinear problems. For example, from Figure 3 we see that the Poisson τ-leap method is not very accurate with τ = 0.4 and quite poor in picking up the second peak with τ = 0.8. On the other hand the RK methods match the peak quite well, albeit with a slight shift in that peak. Furthermore, numerical results from the MAPK cascade simulations show that our methods can run an order of magnitude faster than the Poisson τ-leap and still give the same accuracy in the results.
Finally, we note that we could extend our RK methods to allow more than one set of Poisson random variables to be simulated per step. We imagine that this would allow even bigger stepsizes but at the cost of taking more simulation time in that the additional Poisson sampling is expensive. We emphasise that although our analysis of these new methods has been given for unimolecular reactions, the simulations of the nonlinear Schlögl reaction and the MAPK cascade indicate that these methods have a more general applicability and we will consider nonlinear analysis via Taylor series expansions in future work.
PR would like to thank Marta Dies for helpful discussions. PR acknowledges Obra Social "la Caixa" for funding through the Graduate Fellowship program. Support from Spanish MCINN grant CTQ2008-00755/BQU and from EC-funded projects BioBridge (FP6-2005-LIFESCIHEALTH-7 037909), QosCosGrid (FP6-2005-IST-5 033883) and VPH (FP7-2007-IST-223920) is highly appreciated. JVF participates in the COM-BIOMED network.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.