Efficient simulation of stochastic chemical kinetics with the Stochastic BulirschStoer extrapolation method
 Tamás SzékelyJr.^{1}Email author,
 Kevin Burrage^{1, 2},
 Konstantinos C Zygalakis^{3} and
 Manuel Barrio^{4}Email author
DOI: 10.1186/17520509871
© Székely et al.; licensee BioMed Central Ltd. 2014
Received: 22 January 2013
Accepted: 5 June 2014
Published: 18 June 2014
Abstract
Background
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 BulirschStoer method is an established method for solving ordinary differential equations that possesses both of these qualities.
Results
In this paper, we present the Stochastic BulirschStoer method, 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 BulirschStoer 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 BulirschStoer 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.
Keywords
Stochastic simulation Discrete stochastic methods BulirschStoer τleap Highorder methodsBackground
Microscopic processes with few interacting components can have considerable effects at the macroscopic scale [1–3]. Stochasticity is a defining property of these processes, which can have so few component particles that random fluctuations dominate their behaviour [4, 5]. Stochastic simulation methods take proper account of these fluctuations, as opposed to deterministic methods that assume a system does not deviate from its mean behaviour [6]; although deterministic methods can often be useful for an approximate description of the dynamics of a system, their results are not always representative [7, 8].
A common stochastic modelling approach is to consider the system as a continuoustime 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 a Poisson random variable at each step. The original stepsize selection procedure has since been modified to improve accuracy [13, 14]. To deal with issues of negative populations, Tian and Burrage [15] and Chatterjee et al.[16] introduced the binomial τleap, which samples a binomial random variable at each step. In addition, a newer binomial τleap [17] and a multinomial τleap [18] have since been proposed.
In his seminal paper, Gillespie also proposed the midpoint τleap, a higherorder method that allows for larger timesteps by reducing the inherent bias of the τleap [12]. More recently, other higherorder methods have been developed, such as the unbiased τleap [19], randomcorrected τleap [20], θtrapezoidal τleap [21], and extrapolated τleap [22]. Rather than focussing on adaptively optimising the timestep to reduce processing time, these methods instead improve their order of accuracy so that they find more accurate results for a given stepsize. Because of this, they can use larger timesteps for a desired error level, reducing processing time.
In this paper, we introduce a new adaptivestepsize method for simulating discrete Markov paths, which we call the Stochastic BulirshStoer (SBS) method. This is inspired by the deterministic method of the same name, a very accurate method for solving ODEs, based on Richardson extrapolation. Its high accuracy due to extrapolation, and its ability to adaptively maximise the timestep make the BulirschStoer method one of the most powerful ODE solvers. Possessing these same advantages, our SBS method is a very efficient and accurate new discrete stochastic numerical method.
The SBS calculates several approximations for the expected number of reactions occurring per timestep τ using stepsizes τ/2, τ/4, and so on, and extrapolates these to arrive at a very accurate estimate; the state of the system is then found by sampling a Poisson distribution with this parameter. The SBS is in some ways similar to the extrapolated τleap methods we proposed in a previous paper [22]. These involve running simulations with a τleap method of choice over the full time period of interest, and then taking moments and extrapolating them. The SBS is also based on extrapolation, but the extrapolation is carried out inside each timestep, rather than at the end of the simulation, allowing τ to be optimised at each step.
Overview of stochastic methods
We start with a chemical system of N species and M reactions, interacting in a fixed volume Ω at constant temperature, that is both wellstirred and homogeneous. Thus we assume that individual molecules undergo both reactive collisions and nonreactive ones, and the latter is more frequent than the former, mixing the molecules thoroughly [11]. Individual molecules are not tracked, rather it is their total numbers that we are interested in. These are stored in an N×1 state vector, x≡X(t)≡(X _{1},…,X _{ N })^{ T }, that contains the integer number of each type of molecule at some time t. Reactions are represented by an N×M matrix consisting of stoichiometric vectors ν _{ j }≡(ν _{1j },…,ν _{ N j })^{ T },j=1,…,M, which dictate how each reaction changes the system state, and an M×1 vector of propensity functions a _{ j }(x), where a _{ j }(x)d t gives the probabilities of each reaction occurring in an infinitesimal time interval dt. Together, these three variables fully characterise the chemical system as it evolves through time. In this paper, we adopt the following notation: a bold font variable refers to an N×1 vector, e.g. X(t), and unless otherwise specified the indices i=1,…,N and j=1,…,M.
The SSA is a statistically exact method for generating Monte Carlo paths. That is, a histogram built up from an infinite number of simulations of the SSA will be identical to the true histogram of the system. It is the stochastic method of choice for many researchers, but it has one main limitation: as with other stochastic methods, many realisations (usually starting at 10^{4} or 10^{5}) must be simulated to get a reasonable idea of the histogram shape, and SSA simulations can be slow depending on the problem.
This method is effective but a little crude: unless simplicity is the most important consideration or τ must be fixed, as is the case when used in the context of Richardson extrapolation [22], it is advisable to use more advanced τleap methods.
Adaptively changing τ at each step can give even greater gains in speed, and this is easily introduced into Algorithm 2. A successful approach in current implementations is to find τ such that the mean and variance of the change in propensities over [t,t+τ) are bounded by some fraction ε≪1 of a _{ j }(X _{ n }). The advantage of this is that τ is controlled to stick more closely to the leap condition, ensuring better accuracy, while at the same time maximimising τ for faster simulations. There have been several successive improvements for best selecting τ[12–14], and these methods can achieve a very high efficiency.
Methods
Extrapolation
where the e _{ k } are constant vectors and depend only on the final time and k _{1}<k _{2}<k _{3},…. Eq. (2) tells us that this method has order of accuracy k _{1}.
Neville table built from q initial approximations Y T τ _{ 1 } ,…, Y T τ _{ q } with order k _{ 1 } (first column) and extrapolated to find a solution of order k _{ q } , that is ${Y}_{T}^{{\tau}_{1},{\tau}_{q}}$
Order  k _{1}  k _{2}  k _{3}  …k _{ q } 

${\mathbf{Y}}_{T}^{{\tau}_{1}}$  
${\mathbf{Y}}_{T}^{{\tau}_{1},{\tau}_{2}}$  
Approximate  ${\mathbf{Y}}_{T}^{{\tau}_{2}}$  ${\mathbf{Y}}_{T}^{{\tau}_{1},{\tau}_{3}}$  
solutions  ${\mathbf{Y}}_{T}^{{\tau}_{2},{\tau}_{3}}$  ⋮  $\dots \phantom{\rule{1em}{0ex}}{\mathbf{Y}}_{T}^{{\tau}_{1},{\tau}_{q}}$  
${\mathbf{Y}}_{T}^{{\tau}_{3}}$  ⋮  ${\mathbf{Y}}_{T}^{{\tau}_{q2},{\tau}_{q}}$  
⋮  ${\mathbf{Y}}_{T}^{{\tau}_{q1},{\tau}_{q}}$  
${\mathbf{Y}}_{T}^{{\tau}_{q}}$ 
where p=τ _{ q−r }/τ _{ q−r+1} and k _{ q } is the order of the solution at column q (c.f. Eq. (2)), and r=1,…,q−1. This can be generalised to any order method with any appropriate error expansion. The only condition for extrapolation is existence of an error expansion of the form in Eq. (2).
BulirschStoer method
We give a brief overview of the deterministic BulirschStoer method here; Ref. [28] has an excellent description of the algorithm, as well as a guide to its implementation. At each step, a column of the Neville table, k, in which we expect the approximate solutions to have converged, as well as a stepsize τ are selected. The Neville table is then built up by running k MMPs, with stepsizes ${\widehat{\tau}}_{1}=\tau /2,\dots ,{\widehat{\tau}}_{q}=\tau /{n}_{q}$, where n _{ q }=2q,q=1,2,…,k and successively extrapolating the appropriate numerical approximations. The convergence of the solutions is evaluated based on the internal consistency of the Neville table, that is, the difference between the most accurate solution in column k and that in column k−1: from Table 1, this is $\Delta \mathbf{Y}(k,k1)={\mathbf{Y}}_{\tau}^{{\widehat{\tau}}_{1},{\widehat{\tau}}_{k}}{\mathbf{Y}}_{\tau}^{{\widehat{\tau}}_{2},{\widehat{\tau}}_{k}}$. As successive initial approximations ${\mathbf{Y}}_{\tau}^{{\widehat{\tau}}_{q}}$ are added to the first column, the extrapolated results in each new column converge to the true solution and Δ Y(k,k−1) shrinks. The final approximation at column k is acceptable if e r r _{ k }≤1, where e r r _{ k } is a scaled version of Δ Y(k,k−1) (see Appendix: Stochastic BulirschStoer full algorithm for more detail). If e r r _{ k }>1, the step is rejected and redone with $\tau =\frac{\tau}{2}$.
In a practical implementation, the initial step tests over q=1,…,k _{ m a x }, where k _{ m a x } is usually set as eight, in order to establish the k necessary to achieve the required accuracy and ensure the stepsize is reasonable; subsequent steps then test for convergence only in columns k−1,k and k+1 [28]. Because of its accuracy, the steps taken by the BulirschStoer method can be relatively large compared to other numerical solvers. τ is changed adaptively at each step, and is chosen to minimise the amount of work done (i.e. function evaluations $\widehat{n}+1$ of the MMP) per unit stepsize. In this way the BulirschStoer method adapts its order and stepsize to maximise both accuracy and computational efficiency.
Stochastic BulirschStoer method
The SBS method is based on its deterministic counterpart described in the previous section. There are some key issues that must be addressed in order to successfully adapt it into a stochastic method. The two most important ones are interlinked: first, what quantity should be calculated at each step, and second, how can stochasticity be introduced into the picture? The deterministic BulirschStoer method calculates X _{ n+1} from X _{ n } using the MMP to find the intermediate stages over [t,t+τ). However, stochasticity cannot simply be added to this scheme either inside or outside the MMP, as this would interfere with the extrapolation necessary for the Neville table. In order to update the state vector as in Eq. (1), we must find the number of reactions per step.
it is clear that the quantity we must calculate is $\underset{{t}_{n}}{\overset{{t}_{n}+\tau}{\int}}a\left(\mathbf{X}\right(t\left)\right)\mathit{\text{dt}}$, in order to then take a Poisson sample for the update (the τleap method approximates this as a(X(t _{ n }))τ). Thus, rather than calculating X _{ n+1} directly using the MMP, we need an accurate way to find the integral of the propensity functions over each step. Proceeding in a somewhat similar way to Algorithm 3, we arrive at Algorithm 4: the intermediate stages are found using the MMP, and the propensities calculated at each stage. These intermediate propensities are then fed into a composite trapezoidal method to give an accurate estimate of the integral.
This is our approximation to the underlying probability density function at each step. Combined with the extrapolation mechanism described previously and a way to adapt the stepsize, we have the full SBS method.
where τ is current timestep, τ _{ k } is the hypothetical next timestep for Romberg table column k and S _{1} and S _{2} are safety parameters, introduced in the next paragraph. Here e r r _{ k } is the local error relative to a mixed tolerance, with order $\mathcal{O}\left({\tau}^{2(k1)+1}\right)$ (see Appendix: Stochastic BulirschStoer full algorithm). Its ideal value is exactly one: if it is any smaller than this the step could have been made bigger, and if it is any larger it means our error bound is exceeded and the step must be redone using a smaller τ. At each step, a candidate next timestep is selected for each Neville table column k. As we know some measure of the work done for each k (the number of function evaluations of the MMP), we can calculate the efficiency of each of the k candidate timesteps as the work per unit τ. We then select the candidate τ _{ k } that gives the highest efficiency. For brevity, we have left the full description and stepbystep implementation of the SBS (Algorithm 5) until the Appendix.
The SBS uses several different parameters (see Algorithm 5), all of which have some effect on the results. S _{1} and S _{2} are both safety factors that resize the next timestep by some amount: the smaller they are, the smaller the timestep and the more accurate the solution. As always, however, there is a compromise between stepsize and speed, so one must be careful to optimise the parameters for maximum efficiency. The same is also true for the vectors a _{ t o l }, the absolute error tolerance, and r _{ t o l }, the relative error tolerance. These are used to scale the error that is calculated from the internal consistency of the Romberg table. They are usually set fairly low: around 10^{−6} is common. There is an additional consideration with the SBS, namely that of the column of convergence, k. Even when the safety factors are set high (meaning larger timesteps), the SBS can achieve very high accuracy by simply doing another extrapolation, and going to a higher column. For this reason, the relationship between the safety factors and accuracy is not a direct one, and it is advisable to check the timesteps and column of convergence for each new set of parameters.
Extension: SBSDA
which now replaces Eq. (6). Here ⌊ ⌉ denote rounding to the nearest integer and the value ten has been chosen heuristically as above this value a Poisson sample can be well represented by a Gaussian sample with the appropriate mean and variance.
This approach is somewhat similar to the unbiased τleap [19]. The key difference is that Ref.[19] uses this scheme in the context of a fixedstepsize τleap method, basing the entire method around this scheme. In contrast, the SBSDA is grounded in the BulirschStoer method, which it uses for its stepsize selection and combination of MMP and Richardson extrapolation to find the Poisson increment. The DA approach is only one part of the whole SBSDA, and is only used as an alternative to Eq. (6), in order to also find the variance of the number of reactions occurring over each step.
The SBS and SBSDA methods both calculate the parameter of the Poisson sample but they take different approaches to this. The two key differences are that (1) the SBSDA attempts to correct using the variance of the sampled distribution in order to better approximate the true Poisson parameter when the stepsize is large, but (2) it sacrifices some of its performance because of the inherent inaccuracies of Eqs. (10) and (11) (see Results, Higher order of accuracy and robustness section).
Results and discussion
To illustrate their effectiveness, we apply the SBS and SBSDA methods to four example problems of varying degrees of complexity. We compare them with the popular benchmark of the Euler τleap method (TL; most recent formulation)[14], and we also selected two newer methods that are intended to be representative of the most current, fastest and most accurate methods. These are the θtrapezoidal τleap (TTTL) [21], which has two stages and weak order two, and the unbiased τleap (UBTL) [19], which accurately estimates the mean and variance of the number of reactions that occur during one step. Although the authors of these methods have used fixed stepsizes in their works, we have implemented their methods using the same τadapting scheme as the Euler τleap. This actually makes them more advanced than originally described, but we believe this ensures a fairer comparison with the SBS.
Parameters varied for SBS, SBSDA, TL, TTTL and UBTL for each test system in order from fastest to slowest (left to right in efficiency plots)
Figure  System  Method  Parameters (r _{ t o l } for SBS methods, ε for TL methods)  

SBS  10^{−5}  5×10^{−6}  10^{−6}  5×10^{−7}  2.5×10^{−7}  10^{−7}  
SBSDA  10^{−5}  5×10^{−6}  10^{−6}  5×10^{−7}  10^{−7}  10^{−8}  
Figure 2  Chain decay  TL  0.125  0.1  0.075  0.05  0.04  0.03 
TTTL  0.2  0.15  0.1  0.075  0.05  0.04  
UBTL  4  2  1  0.8  0.6  0.5  
SBS  10^{−2}  10^{−3}  10^{−4}  10^{−5}  10^{−6}  10^{−7}  
SBSDA  10^{−3}  10^{−4}  10^{−5}  10^{−6}  10^{−7}  10^{−8}  
Figure 3  MichaelisMenten  TL  0.2  0.15  0.1  0.07  0.05  0.04 
TTTL  0.3  0.25  0.2  0.1  0.075  0.06  
UBTL  7  5  3  2  1.5  1  
SBS  10^{−5}  8×10^{−6}  4×10^{−6}  10^{−6}  5×10^{−7}  10^{−7}  
SBSDA  2×10^{−4}  1.5×10^{−4}  10^{−4}  9×10^{−5}  8×10^{−5}  6×10^{−5}  
Figure 4  Schlögl  TL  0.046  0.044  0.042  0.04  0.0375  0.035 
TTTL  0.06  0.056  0.052  0.048  0.046  0.044  
UBTL  0.3  0.28  0.26  0.24  0.22  0.2  
SBS  10^{−5}  10^{−6}  10^{−7}  10^{−8}  10^{−9}  10^{−10}  
SBSDA  10^{−5}  10^{−6}  10^{−7}  10^{−8}  10^{−9}  10^{−10}  
Figure 5  Enzymes  TL  0.15  0.1  0.07  0.05  0.03  0.02 
TTTL  0.3  0.2  0.12  0.08  0.06  0.04  
UBTL  5  3  1.5  0.8  0.4  0.3 
The same number of simulations were run with all methods. We plotted probability density functions (PDFs) for each species and compared them to ones obtained from a reference set of 10^{6} SSA simulations, all generated from histograms using identical bins. We defined the histogram error as the L ^{1} distance between the probabilities of each method and the SSA in each bin. The runtime is the time taken to run a single simulation, obtained by dividing the total runtime by the number of simulations.
We show the probability distributions of all the simulation methods, as well as plots of histogram error versus (single) runtime. We refer to the latter as ‘efficiency’ plots, as they clearly indicate some measure of computational efficiency. If a method is both fast and has low error, it is efficient: its points are concentrated towards the origin. In contrast, points to the top right indicate low efficiency (i.e. a slow and inaccurate method).
Chain decay system
Efficiencies of each method as calculated according to Eq. ( 12) for each test system (higher is better)
System  SBS  SBSDA  TL  TTTL  UBTL  Bin sizes (for species 1,2,…N) 

Chain decay  14.7  22.3  0.4  4.4  14.1  2, 5, 5 
MichaelisMenten  7.1  11.7  2.3  7.7  1.2  5, 5, 5, 5 
Schlögl  19.1  7.1  18.8  18.8  0.3  10 
Enzymes  21.9  50.5  12.6  6.3  3.0  100, 100, 50, 50, 50, 50, 50, 50 
MichaelisMenten system
Schlögl system
and species A and B are held constant at 10^{5} and 2×10^{5} units respectively. We used the initial condition X(0)=250, which is an intermediate value between the two stable states. The simulation time was T=10, and we ran 10^{5} simulations for each method. The SBS safety factors were S _{1}=S _{2}=0.05, those for the SBSDA as S _{1}=S _{2}=0.125, and a _{ t o l }=10^{−6} for both.
Mutually inhibiting enzymes system
Further comparisons
All of our test systems have more than one species, and so far we have only presented results for X _{1}. This can often be unrepresentative of the full picture. The chain decay system is a clear example of this. Additional file 1: Figure A1 shows the efficiency plots for all three species. Only looking at X _{1} could lead one to think that the UBTL is the most efficient method for simulating this system. But including the other two species reveals that the SBSDA is, in fact, the most efficient overall. This is important, because it is clear that factors such as linear/nonlinear propensities, population size and stiffness all affect each reaction and species in a different way. Thus it is overall performance we are interested in.
This varies for each system, and is comparable only when the bins used to calculate errors are identical. In other words, the values are a direct comparison of the efficiency of each method for each test system, but should not be used across different test systems. Table 3 compares the efficiencies of each simulation method for every test system.
Additional file 1: Figures A1 to A4 contain a full picture of our computational results for all four example systems and all simulation methods. The overarching trend was the following: the SBS was very accurate, returning the lowest error in many cases. The TL was unexpectedly efficient for some systems. The TTTL also achieved good efficiency for longer runtimes, but the SBS had a flatter efficiency curve than the TTTL, with the TTTL quickly losing accuracy at lower runtimes even when it was more accurate than the SBS at higher runtimes. A clear trend emerges: the SBS is a very accurate method. Moreover, it is the most efficient method we tested, maintaining its accuracy better at low runtimes than the other methods.
SBS methods excel when we want a short runtime with high accuracy. In this case, we set the safety factors high, allowing large steps and a corresponding increase in extrapolations. This retains a high accuracy, whilst reducing runtime because of the large timesteps. In contrast, when we allow a longer runtime, we set the safety factors low, restricting the stepsize and removing the need for higher extrapolation. In many of our test examples, we have seen the SBS only using one extrapolation throughout the simulation. This is a waste of the extrapolation capability of the SBS, and it is no surprise that in these cases it is not the most efficient method, especially as the stepsize adaptation scheme adds some overhead to each step.
Higher order of accuracy and robustness
A thorough study of the order properties of τleaping methods was first given by Rathinam et al.[35], who showed that for linear reactions the Euler τleap method is weak order one in the moments under the scaling τ→0. This analysis was extended by Li [36] to nonlinear propensity functions by considering SDEs driven by Poisson random measures (see also [37]). Li showed that the Euler τleap method is precisely the Euler method applied to this SDE and hence inherits the properties of strong order half and weak order one. However, there are issues with using the scaling condition τ→0 as the τleap condition requires that $\sum {a}_{j}\left(\mathbf{X}\right)\tau \gg 1$. Anderson et al.[38] overcame this scaling condition by considering order under a large volume scaling Ω→∞. In this case, by letting $X/\Omega =\mathcal{O}\left(1\right)$ and with τ=Ω ^{−β },0<β<1, global strong and weak order convergence can be established. Hu et al.[39] investigate these issues in greater detail through the use of rooted tree expansions of the local truncation errors for the moments and covariance, thus generalising the approach first applied to SDEs by Burrage and Burrage [40]. This analysis shows that while some τleap methods may have higher order moments (for instance, the midpoint τleap has order two moments for linear systems), their covariance is invariably of unit order, unless this is specially taken into consideration (as with the TTTL method, which has order two moments and covariance). As Hu et al.[39] point out, these issues arise as a consequence of the differences between the infinitesimal generators for deterministic ODEs and jump processes.
It is wellknown that the BulirschStoer method has a high order of accuracy: this is the reason it is able to use large steps whilst still finding very accurate solutions. This is because of the Richardson extrapolation that is used at each step on the MMP solutions (which themselves have order two as well as an error expansion containing only even powers of $\widehat{\tau}$, resulting in very high order solutions with little work). In contrast, rather than the MMP solutions for X(t), the SBS instead extrapolates at each step the deterministic quantities $\Delta {a}^{{\widehat{\tau}}_{q}}({t}_{n},{t}_{n}+\tau )$ (or the mean and variance of K _{ j } given by Eqs. (10) and (11) in the case of the SBSDA), calculated using the composite trapezoidal rule, which also has a known error expansion. Thus the extrapolation is performed on a deterministic variable: the mean of the Poisson update.
We investigated the behaviour of the SBS and SBSDA methods on two simple systems: first, the linear system X→2X, X(0)=1000, and second, the nonlinear system X+Y→∅, X(0)=Y(0)=10000. As the SBS changes both stepsize and Romberg table column k (i.e. order of accuracy) adaptively, we used a restricted version, which had both a fixed stepsize and k. We ran simulations with both SBS and SBSDA, for k=1 and k=2, that is no extrapolation and one extrapolation, respectively. In addition, we also used the TL and TTTL methods for comparison, as they have known weak orders of accuracy (one and two, respectively). The gradients of the errors for the different methods are computed based on a linear leastsquares regression of the data points.
where ε>0 is a smallnoise term. It is wellknown that Langevin SDEs represent an intermediate regime between discrete stochastic chemical kinetics and the deterministic regime, arising as the number of molecules X in the system increases. In particular, ε behaves as $\frac{1}{\sqrt{X}}$[11]. For such systems, Milstein and Tretyakov [42] showed that the global weak order of numerical methods to solve the above SDE has the general form $\mathcal{O}({\tau}^{p}+{\tau}^{q}{\epsilon}^{r})$, where q<p. When noise is ignored (ε=0), the SDE becomes an ODE and the weak order of its approximate solution is just the deterministic order term $\mathcal{O}\left({\tau}^{p}\right)$. Milstein and Tretyakov [41] also performed an analysis in the strong sense, and again found the general form of the global strong error to be $\mathcal{O}({\tau}^{p}+{\tau}^{q}{\epsilon}^{r}),q<p$. In addition, Buckwar et al.[43] also examined smallnoise SDEs in a strong sense for some wellknown classes of RungeKutta methods. The implication of the extra term in the stochastic order is that although the underlying deterministic order of the method may be high, the stochastic order is restricted by the noise term. However, when the noise is small, this term will also become small, thus allowing the stochastic order to increase, possibly even up to the deterministic order $\mathcal{O}\left({\tau}^{p}\right)$. This is also the case if the stepsize is large. In fact, it occurs over a range of values of τ and ε: it is trivial to see that the condition for the deterministic term to dominate is $\tau \gg {\epsilon}^{\frac{r}{pq}}$.
Now, a standard way of mathematically investigating discrete stochastic methods is to analyse SDEs with jumps; this is the approach of Li [36]. In particular, Li [36] shows that the Euler discretisation of such an SDE is the Euler τleap method. Our hypothesis is that the analysis of Milstein and Tretyakov [42] is also applicable to SDEs with jumps; this then tells us something about the behaviour of discrete stochastic methods when noise levels are medium or small. A smallnoise analysis similar to that of Milstein and Tretyakov [42] for the SBS is beyond the scope of this paper but we postulate that there is a smallnoise error expansion in both τ and ε for the SBS method. The SBS is most useful when applied to systems with relatively larger biochemical populations (thus small noise), as it is in these cases that the SSA is prohibitively slow and an approximate method is necessary. Combined with the fact that the SBS often uses large stepsizes, this implies that in many systems of interest, the global weak order of the SBS may, in fact, not be far from its deterministic (high) order. This also explains the behaviour of the SBS in our numerical tests, where the timesteps were large and the populations moderate (implying moderate noise), thus making it likely that the condition for the deterministic order term to dominate was met.
Implementation issues
There are two distinct approaches to determining τ(0): first, as described previously, we could set τ(0) to an arbitrary value and run the initial step through as many columns as necessary (up to k _{ m a x }) until it finds the required accuracy. Should τ(0) be so large that it drives the populations negative, it would also be reduced here until it reaches a more suitable size for the given problem. In addition, if τ(0) is still larger than its optimum value, it is reduced over the next several steps until it has reached this optimum value (and vice versa if it is too small). This is the standard approach for the deterministic BulirschStoer method, and it is the one we have taken in our simulations. However, in the stochastic regime there is another approach: we could set τ(0) as some multiple of 1/a _{0}(x _{0}) (the expected size of an SSA step in state x _{0}[12]), along with an initial guess of the Romberg table column to aim for. As 1/a _{0}(x _{0}) is very small, this is a more conservative approach, but τ is increased to its optimum value over the first few steps. It could be useful for systems that are very stiff, or that oscillate, whose timestep must be very small at certain parts of the solution domain and larger timesteps could result in large errors. There seems to be no substantial difference in accuracy between the two approaches, and we believe both are equally valid.
Conclusions
Our results have shown that the SBS is generally a very accurate method, at least comparable to or, in most cases, better than its competitors. However, the real strength of the SBS is this accuracy combined with the fact that its efficiency curve has a relatively low gradient; in other words, it is an accurate method that loses little of its accuracy as it is speeded up, allowing for fast, robust and accurate simulations. This is because as runtime is shortened, the SBS uses more and more extrapolations to maintain its accuracy. At the same time, the use of larger timesteps means less overhead overall, allowing the SBS to be very efficient. It is in such parameter regimes that the SBS can achieve its full potential. In addition to this, we believe the SBS is also able to achieve high weak order (in the moments; the variance remains one) in the smalltomoderate noise regime, that is when the number of molecules in the system is moderate to large, and also when the timesteps are large compared to the noise level. Its performance in this regime is accelerated as more and more extrapolations are performed, giving it exceptional accuracy.
As the SBS is an explicit method, it is not necessarily suited for solving especially stiff problems. In such cases, RungeKutta methods with larger regions of stability, such as the stochastic RungeKutta method [44], are more ideal, as well as implicit or multiscale methods [45–47]. The initial stepsize of the SBS should be chosen appropriately, as it may be possible for the SBS to settle in a higherstepsize regime, which could affect accuracy, or a lowstepsize regime, which could affect runtime. In addition, τ(0) should be chosen such that it is within the stability region of the modified midpoint method. Running a few preliminary simulations can help choose τ(0).
In previous work, we have extended Richardson extrapolation into the discrete stochastic regime [22]. In this framework, full simulations with fixed stepsize are run over t=[0,T], and their moments are extrapolated to find accurate approximations to the moments at time T. In contrast, the SBS uses extrapolation within each timestep and varies τ to optimise efficiency. Thus the SBS is a complementary approach to extrapolated τleap methods that has two advantages: first, the stepsize can be adapted to lower runtime and eliminate the need for finding a suitable range of fixed stepsizes; second, the SBS returns an entire histogram, rather than just the moments. This can be desirable in many cases, especially if the solutions do not follow a simple distribution such as a Gaussian or Poisson, or have multiple stable states.
In this paper we have introduced a new efficient and robust simulation method, the Stochastic BulirschStoer method, which can also achieve higher weak order in the moments for certain systems. This is inspired by the deterministic method of the same name, and as such it also boasts the two main advantages of that method: its speed and its high accuracy. We have shown using numerical simulations that for a range of example problems, it is generally the most efficient and robust out of all recent τleap methods that we tested, which are the current stateoftheart in fast stochastic simulation. Thus the SBS is a promising new method to address the need for fast and efficient discrete stochastic methods.
Appendix: Stochastic BulirschStoer full algorithm
where τ _{ q } is a set of hypothetical new stepsizes adjusted from the current stepsize τ. S _{1} and S _{2} are safety factors 0<S _{1},S _{2}<1, that ensure τ is not set too large because of errors in the MMP and composite trapezoidal rule approximations.
where n _{ q }=2q. The optimal column k for the next timestep is given by the lowest W _{ q }, and the optimal stepsize by the corresponding τ _{ q }. In reality, after the initial step only columns k−1, k and k+1 are tested for convergence, as otherwise the convergence is likely to be an artifact or the timestep is far off its optimal size. This helps reduce the runtime but makes the implementation more complicated.
Abbreviations
 SSA:

Stochastic simulation algorithm
 SBS:

Stochastic BulirschStoer
 ODE:

Ordinary differential equation
 MMP:

Modified midpoint method
 DA:

Degree of advancement
 PDF:

Probability density function
 TL:

(Euler) τleap
 TTTL:

θtrapezoidal τleap
 UBTL:

Unbiased τleap.
Declarations
Acknowledgements
TSz was supported by the Engineering and Physical Sciences Research Council through the Systems Biology Doctoral Training Centre, University of Oxford.
Authors’ Affiliations
References
 Arkin A, Ross J, McAdams HH: Stochastic kinetic analysis of developmental pathway bifurcation in phage λinfected Escherichia coli cells. Genetics. 1998, 149: 16331648.PubMed CentralPubMedGoogle Scholar
 Kaern M, Elston T, Blake W, Collins J: Stochasticity in gene expression: from theories to phenotypes. Nat Rev Genet. 2005, 6: 451464.View ArticlePubMedGoogle Scholar
 Munsky B, Neuert G, van Oudenaarden A: Using gene expression noise to understand gene regulation. Science. 2012, 336: 183PubMed CentralView ArticlePubMedGoogle Scholar
 Fedoroff N, Fontana W: Small numbers of big molecules. Science. 2002, 297: 11291131.View ArticlePubMedGoogle Scholar
 Elowitz M, Levine A, Siggia E, Swain P: Stochastic gene expression in a single cell. Science. 2002, 297: 11831186.View ArticlePubMedGoogle Scholar
 Gillespie DT, Mangel M: Conditioned averages in chemical kinetics. J Chem Phys. 1981, 75: 704709.View ArticleGoogle Scholar
 Goutsias J: Classical versus stochastic kinetics modeling of biochemical reaction systems. Biophys J. 2007, 92: 23502365.PubMed CentralView ArticlePubMedGoogle Scholar
 Wilkinson DJ: Stochastic modelling for quantitative description of heterogeneous biological systems. Nat Rev Genet. 2009, 10: 122133.View ArticlePubMedGoogle Scholar
 Kurtz TG: Limit theorems for sequences of jump Markov processes approximating ordinary differential processes. J Appl Probab. 1971, 8: 344356.View ArticleGoogle Scholar
 Gillespie DT: Exact stochastic simulation of coupled chemical reactions. J Phys Chem. 1977, 81: 23402361.View ArticleGoogle Scholar
 Gillespie DT: Stochastic simulation of chemical kinetics. Annu Rev Phys Chem. 2007, 58: 3555.View ArticlePubMedGoogle Scholar
 Gillespie DT: Approximate accelerated stochastic simulation of chemically reacting systems. J Chem Phys. 2001, 115: 17161733.View ArticleGoogle Scholar
 Gillespie DT, Petzold LR: Improved leapsize selection for accelerated stochastic simulation. J Chem Phys. 2003, 119: 82298234.View ArticleGoogle Scholar
 Cao Y, Gillespie DT, Petzold LR: Efficient step size selection for the tauleaping simulation method. J Chem Phys. 2006, 124: 044109View ArticlePubMedGoogle Scholar
 Tian TH, Burrage K: Binomial leap methods for simulating stochastic chemical kinetics. J Chem Phys. 2004, 121: 1035610364.View ArticlePubMedGoogle Scholar
 Chatterjee A, Vlachos DG, Katsoulakis MA: Binomial distribution based tauleap accelerated stochastic simulation. J Chem Phys. 2005, 122: 024112View ArticlePubMedGoogle Scholar
 Peng X, Zhou W, Wang Y: Efficient binomial leap method for simulating chemical kinetics. J Chem Phys. 2007, 126: 224109View ArticlePubMedGoogle Scholar
 Pettigrew MF, Resat H: Multinomial tauleaping method for stochastic kinetic simulations. J Chem Phys. 2007, 126: 084101View ArticlePubMedGoogle Scholar
 Xu Z, Cai X: Unbiased tauleap methods for stochastic simulation of chemically reacting systems. J Chem Phys. 2008, 128: 154112View ArticlePubMedGoogle Scholar
 Hu Y, Li T: Highly accurate tauleaping methods with random corrections. J Chem Phys. 2009, 130: 124109View ArticlePubMedGoogle Scholar
 Hu Y, Li T, Min B: A weak second order tauleaping method for chemical kinetic systems. J Chem Phys. 2011, 135: 024113View ArticlePubMedGoogle Scholar
 Székely T, Burrage K, Erban R, Zygalakis KC: A higherorder numerical framework for stochastic simulation of chemical reaction systems. BMC Syst Biol. 2012, 6: 85PubMed CentralView ArticlePubMedGoogle Scholar
 Richardson LF: The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam. Phil Trans Roy Soc Lond. 1910, 210: 307357.View ArticleGoogle Scholar
 Hairer E, Nørsett SP, Wanner G: Solving Ordinary Differential Equations: Nonstiff Problems. 2nd edition. 1993, Berlin: SpringerVerlagGoogle Scholar
 Bulirsch R, Stoer J: Numerical treatment of ordinary differential equations by extrapolation methods. Numerische Mathematik. 1966, 8: 113.View ArticleGoogle Scholar
 Deuflhard P: Recent progress in extrapolation methods for ordinary differential equations. SIAM Rev. 1985, 27 (4): 505535.View ArticleGoogle Scholar
 Gragg WB: On extrapolation algorithms for ordinary initial value problems. SIAM J Numer Anal. 1965, 2: 384403.Google Scholar
 Press WH, Teukolsky SA, Vetterling WT, Flannery BP: Numerical Recipes in C: The Art of Scientific Computing. 2nd edition. 1992, Cambridge: Cambridge University PressGoogle Scholar
 Kurtz TG: Strong approximation theorems for density dependent Markov chains. Stochastic Processes Appl. 1978, 6: 223240.View ArticleGoogle Scholar
 van Kampen NG: Stochastic Processes in Physics and Chemistry. 3rd edition. 2007, Amsterdam: ElsevierGoogle Scholar
 Goutsias J, Jenkinson G: Markovian dynamics on complex reaction networks. Phys Rep. 2013, 529: 199264.View ArticleGoogle Scholar
 Kurtz TG: Representations of Markov processes as multiparameter time changes. Ann Probab. 1980, 8: 682715.View ArticleGoogle Scholar
 MarquezLago T, Burrage K: Binomial tauleap spatial stochastic simulation algorithm for applications in chemical kinetics. J Chem Phys. 2007, 127: 104101View ArticlePubMedGoogle Scholar
 Elf J, Ehrenberg M: Spontaneous separation of bistable biochemical systems into spatial domains of opposite phases. Syst Biol. 2004, 1: 230236.View ArticleGoogle Scholar
 Rathinam M, Petzold LR, Cao Y, Gillespie DT: Consistency and stability of tauleaping schemes for chemical reaction systems. Multiscale Model Simul. 2005, 4 (3): 867895.View ArticleGoogle Scholar
 Li T: Analysis of explicit tauleaping schemes for simulating chemically reacting systems. Multiscale Model Simul. 2007, 6 (2): 417436.View ArticleGoogle Scholar
 Burrage K, Tian T: Poisson RungeKutta methods for chemical reaction systems. Advances in Scientific Computing and Applications. Edited by: Sun YLW, Tang T. 2004, Beijing/New York: Science Press, 8296.Google Scholar
 Anderson DF, Ganguly A, Kurtz TG: Error analysis of tauleap simulation methods. Ann Appl Probab. 2011, 21 (6): 22262262.View ArticleGoogle Scholar
 Hu Y, Li T, Min B: The weak convergence analysis of tauleaping methods: revisited. Comm Math Sci. 2011, 9: 965996.View ArticleGoogle Scholar
 Burrage K, Burrage PM: Order conditions of stochastic Runge–Kutta methods by BSeries. SIAM J Numer Anal. 2000, 38 (5): 16261646.View ArticleGoogle Scholar
 Milstein GN, Tretyakov MV: Numerical methods in the weak sense for stochastic differential equations with small noise. SIAM J Numer Anal. 1997, 34: 21422167.View ArticleGoogle Scholar
 Milstein GN, Tretyakov MV: Meansquare numerical methods for stochastic differential equations with small noises. SIAM J Sci Comput. 1997, 18: 10671087.View ArticleGoogle Scholar
 Buckwar E, Rößler A, Winkler R: Stochastic RungeKutta methods for Ito SODEs with small noise. SIAM J Sci Comput. 2010, 32: 17891808.View ArticleGoogle Scholar
 Rué P, VillaFreixà J, Burrage K: Simulation methods with extended stability for stiff biochemical kinetics. BMC Syst Biol. 2010, 4: 110123.PubMed CentralView ArticlePubMedGoogle Scholar
 Cao Y, Gillespie DT, Petzold LR: The adaptive explicitimplicit tauleaping method with automatic tau selection. J Chem Phys. 2007, 126: 224101View ArticlePubMedGoogle Scholar
 Goutsias J: Quasiequilibrium approximation of fast reaction kinetics in stochastic biochemical systems. J Chem Phys. 2005, 122: 184102View ArticlePubMedGoogle Scholar
 MacNamara S, Burrage K, Sidje R: Multiscale modeling of chemical kinetics via the master equation. Multiscale Model Simul. 2008, 6 (4): 11461168.View ArticleGoogle Scholar
Copyright
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 credited.