Efficient simulation of stochastic chemical kinetics with the Stochastic Bulirsch-Stoer extrapolation method
© Székely et al.; licensee BioMed Central Ltd. 2014
Received: 22 January 2013
Accepted: 5 June 2014
Published: 18 June 2014
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.
In this paper, we present the Stochastic Bulirsch-Stoer 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.
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.
KeywordsStochastic simulation Discrete stochastic methods Bulirsch-Stoer τ-leap High-order methods
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 ; 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 continuous-time Markov jump process . The stochastic simulation algorithm (SSA) of Gillespie  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 ).
The first of these is commonly called the Euler or Poisson τ-leap ; 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  and Chatterjee et al. introduced the binomial τ-leap, which samples a binomial random variable at each step. In addition, a newer binomial τ-leap  and a multinomial τ-leap  have since been proposed.
In his seminal paper, Gillespie also proposed the midpoint τ-leap, a higher-order method that allows for larger timesteps by reducing the inherent bias of the τ-leap . More recently, other higher-order methods have been developed, such as the unbiased τ-leap , random-corrected τ-leap , θ-trapezoidal τ-leap , and extrapolated τ-leap . 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 adaptive-stepsize method for simulating discrete Markov paths, which we call the Stochastic Bulirsh-Stoer (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 Bulirsch-Stoer 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 . 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 well-stirred and homogeneous. Thus we assume that individual molecules undergo both reactive collisions and non-reactive ones, and the latter is more frequent than the former, mixing the molecules thoroughly . 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 104 or 105) 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 , 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.
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
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).
We give a brief overview of the deterministic Bulirsch-Stoer method here; Ref.  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 , 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 . As successive initial approximations 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 Bulirsch-Stoer full algorithm for more detail). If e r r k >1, the step is rejected and redone with .
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 . Because of its accuracy, the steps taken by the Bulirsch-Stoer 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 of the MMP) per unit stepsize. In this way the Bulirsch-Stoer method adapts its order and stepsize to maximise both accuracy and computational efficiency.
Stochastic Bulirsch-Stoer 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 Bulirsch-Stoer 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 , 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 (see Appendix: Stochastic Bulirsch-Stoer 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 step-by-step 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.
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 . The key difference is that Ref. uses this scheme in the context of a fixed-stepsize τ-leap method, basing the entire method around this scheme. In contrast, the SBS-DA is grounded in the Bulirsch-Stoer 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 SBS-DA, 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 SBS-DA methods both calculate the parameter of the Poisson sample but they take different approaches to this. The two key differences are that (1) the SBS-DA 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 SBS-DA 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), 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) , which has two stages and weak order two, and the unbiased τ-leap (UBTL) , 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, SBS-DA, TL, TTTL and UBTL for each test system in order from fastest to slowest (left to right in efficiency plots)
Parameters (r t o l for SBS methods, ε for TL methods)
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 106 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)
Bin sizes (for species 1,2,…N)
2, 5, 5
5, 5, 5, 5
100, 100, 50, 50, 50, 50, 50, 50
and species A and B are held constant at 105 and 2×105 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 105 simulations for each method. The SBS safety factors were S 1=S 2=0.05, those for the SBS-DA as S 1=S 2=0.125, and a t o l =10−6 for both.
Mutually inhibiting enzymes system
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 SBS-DA is, in fact, the most efficient overall. This is important, because it is clear that factors such as linear/non-linear 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., 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  to non-linear propensity functions by considering SDEs driven by Poisson random measures (see also ). 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 . Anderson et al. overcame this scaling condition by considering order under a large volume scaling Ω→∞. In this case, by letting and with τ=Ω −β ,0<β<1, global strong and weak order convergence can be established. Hu et al. 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 . 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. point out, these issues arise as a consequence of the differences between the infinitesimal generators for deterministic ODEs and jump processes.
It is well-known that the Bulirsch-Stoer 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 , 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 (or the mean and variance of K j given by Eqs. (10) and (11) in the case of the SBS-DA), 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 SBS-DA methods on two simple systems: first, the linear system X→2X, X(0)=1000, and second, the non-linear 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 SBS-DA, 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 least-squares regression of the data points.
where ε>0 is a small-noise term. It is well-known 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 . For such systems, Milstein and Tretyakov  showed that the global weak order of numerical methods to solve the above SDE has the general form , 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 . Milstein and Tretyakov  also performed an analysis in the strong sense, and again found the general form of the global strong error to be . In addition, Buckwar et al. also examined small-noise SDEs in a strong sense for some well-known classes of Runge-Kutta 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 . 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 .
Now, a standard way of mathematically investigating discrete stochastic methods is to analyse SDEs with jumps; this is the approach of Li . In particular, Li  shows that the Euler discretisation of such an SDE is the Euler τ-leap method. Our hypothesis is that the analysis of Milstein and Tretyakov  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 small-noise analysis similar to that of Milstein and Tretyakov  for the SBS is beyond the scope of this paper but we postulate that there is a small-noise 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.
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 Bulirsch-Stoer 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), 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.
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 small-to-moderate 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, Runge-Kutta methods with larger regions of stability, such as the stochastic Runge-Kutta method , 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 higher-stepsize regime, which could affect accuracy, or a low-stepsize 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 . 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 Bulirsch-Stoer 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 state-of-the-art 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 Bulirsch-Stoer 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.
Stochastic simulation algorithm
Ordinary differential equation
Modified midpoint method
Degree of advancement
Probability density function
TSz was supported by the Engineering and Physical Sciences Research Council through the Systems Biology Doctoral Training Centre, University of Oxford.
- Arkin A, Ross J, McAdams HH: Stochastic kinetic analysis of developmental pathway bifurcation in phage λ-infected Escherichia coli cells. Genetics. 1998, 149: 1633-1648.PubMed CentralPubMedGoogle Scholar
- Kaern M, Elston T, Blake W, Collins J: Stochasticity in gene expression: from theories to phenotypes. Nat Rev Genet. 2005, 6: 451-464.View ArticlePubMedGoogle Scholar
- Munsky B, Neuert G, van Oudenaarden A: Using gene expression noise to understand gene regulation. Science. 2012, 336: 183-PubMed CentralView ArticlePubMedGoogle Scholar
- Fedoroff N, Fontana W: Small numbers of big molecules. Science. 2002, 297: 1129-1131.View ArticlePubMedGoogle Scholar
- Elowitz M, Levine A, Siggia E, Swain P: Stochastic gene expression in a single cell. Science. 2002, 297: 1183-1186.View ArticlePubMedGoogle Scholar
- Gillespie DT, Mangel M: Conditioned averages in chemical kinetics. J Chem Phys. 1981, 75: 704-709.View ArticleGoogle Scholar
- Goutsias J: Classical versus stochastic kinetics modeling of biochemical reaction systems. Biophys J. 2007, 92: 2350-2365.PubMed CentralView ArticlePubMedGoogle Scholar
- Wilkinson DJ: Stochastic modelling for quantitative description of heterogeneous biological systems. Nat Rev Genet. 2009, 10: 122-133.View ArticlePubMedGoogle Scholar
- Kurtz TG: Limit theorems for sequences of jump Markov processes approximating ordinary differential processes. J Appl Probab. 1971, 8: 344-356.View ArticleGoogle Scholar
- Gillespie DT: Exact stochastic simulation of coupled chemical reactions. J Phys Chem. 1977, 81: 2340-2361.View ArticleGoogle Scholar
- Gillespie DT: Stochastic simulation of chemical kinetics. Annu Rev Phys Chem. 2007, 58: 35-55.View ArticlePubMedGoogle Scholar
- Gillespie DT: Approximate accelerated stochastic simulation of chemically reacting systems. J Chem Phys. 2001, 115: 1716-1733.View ArticleGoogle Scholar
- Gillespie DT, Petzold LR: Improved leap-size selection for accelerated stochastic simulation. J Chem Phys. 2003, 119: 8229-8234.View ArticleGoogle Scholar
- Cao Y, Gillespie DT, Petzold LR: Efficient step size selection for the tau-leaping simulation method. J Chem Phys. 2006, 124: 044109-View ArticlePubMedGoogle Scholar
- Tian TH, Burrage K: Binomial leap methods for simulating stochastic chemical kinetics. J Chem Phys. 2004, 121: 10356-10364.View ArticlePubMedGoogle Scholar
- Chatterjee A, Vlachos DG, Katsoulakis MA: Binomial distribution based tau-leap accelerated stochastic simulation. J Chem Phys. 2005, 122: 024112-View ArticlePubMedGoogle Scholar
- Peng X, Zhou W, Wang Y: Efficient binomial leap method for simulating chemical kinetics. J Chem Phys. 2007, 126: 224109-View ArticlePubMedGoogle Scholar
- Pettigrew MF, Resat H: Multinomial tau-leaping method for stochastic kinetic simulations. J Chem Phys. 2007, 126: 084101-View ArticlePubMedGoogle Scholar
- Xu Z, Cai X: Unbiased tau-leap methods for stochastic simulation of chemically reacting systems. J Chem Phys. 2008, 128: 154112-View ArticlePubMedGoogle Scholar
- Hu Y, Li T: Highly accurate tau-leaping methods with random corrections. J Chem Phys. 2009, 130: 124109-View ArticlePubMedGoogle Scholar
- Hu Y, Li T, Min B: A weak second order tau-leaping method for chemical kinetic systems. J Chem Phys. 2011, 135: 024113-View ArticlePubMedGoogle Scholar
- Székely T, Burrage K, Erban R, Zygalakis KC: A higher-order numerical framework for stochastic simulation of chemical reaction systems. BMC Syst Biol. 2012, 6: 85-PubMed 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: 307-357.View ArticleGoogle Scholar
- Hairer E, Nørsett SP, Wanner G: Solving Ordinary Differential Equations: Nonstiff Problems. 2nd edition. 1993, Berlin: Springer-VerlagGoogle Scholar
- Bulirsch R, Stoer J: Numerical treatment of ordinary differential equations by extrapolation methods. Numerische Mathematik. 1966, 8: 1-13.View ArticleGoogle Scholar
- Deuflhard P: Recent progress in extrapolation methods for ordinary differential equations. SIAM Rev. 1985, 27 (4): 505-535.View ArticleGoogle Scholar
- Gragg WB: On extrapolation algorithms for ordinary initial value problems. SIAM J Numer Anal. 1965, 2: 384-403.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: 223-240.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: 199-264.View ArticleGoogle Scholar
- Kurtz TG: Representations of Markov processes as multiparameter time changes. Ann Probab. 1980, 8: 682-715.View ArticleGoogle Scholar
- Marquez-Lago T, Burrage K: Binomial tau-leap spatial stochastic simulation algorithm for applications in chemical kinetics. J Chem Phys. 2007, 127: 104101-View ArticlePubMedGoogle Scholar
- Elf J, Ehrenberg M: Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases. Syst Biol. 2004, 1: 230-236.View ArticleGoogle Scholar
- Rathinam M, Petzold LR, Cao Y, Gillespie DT: Consistency and stability of tau-leaping schemes for chemical reaction systems. Multiscale Model Simul. 2005, 4 (3): 867-895.View ArticleGoogle Scholar
- Li T: Analysis of explicit tau-leaping schemes for simulating chemically reacting systems. Multiscale Model Simul. 2007, 6 (2): 417-436.View ArticleGoogle Scholar
- Burrage K, Tian T: Poisson Runge-Kutta methods for chemical reaction systems. Advances in Scientific Computing and Applications. Edited by: Sun YLW, Tang T. 2004, Beijing/New York: Science Press, 82-96.Google Scholar
- Anderson DF, Ganguly A, Kurtz TG: Error analysis of tau-leap simulation methods. Ann Appl Probab. 2011, 21 (6): 2226-2262.View ArticleGoogle Scholar
- Hu Y, Li T, Min B: The weak convergence analysis of tau-leaping methods: revisited. Comm Math Sci. 2011, 9: 965-996.View ArticleGoogle Scholar
- Burrage K, Burrage PM: Order conditions of stochastic Runge–Kutta methods by B-Series. SIAM J Numer Anal. 2000, 38 (5): 1626-1646.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: 2142-2167.View ArticleGoogle Scholar
- Milstein GN, Tretyakov MV: Mean-square numerical methods for stochastic differential equations with small noises. SIAM J Sci Comput. 1997, 18: 1067-1087.View ArticleGoogle Scholar
- Buckwar E, Rößler A, Winkler R: Stochastic Runge-Kutta methods for Ito SODEs with small noise. SIAM J Sci Comput. 2010, 32: 1789-1808.View ArticleGoogle Scholar
- Rué P, Villa-Freixà J, Burrage K: Simulation methods with extended stability for stiff biochemical kinetics. BMC Syst Biol. 2010, 4: 110-123.PubMed CentralView ArticlePubMedGoogle Scholar
- Cao Y, Gillespie DT, Petzold LR: The adaptive explicit-implicit tau-leaping method with automatic tau selection. J Chem Phys. 2007, 126: 224101-View ArticlePubMedGoogle Scholar
- Goutsias J: Quasiequilibrium approximation of fast reaction kinetics in stochastic biochemical systems. J Chem Phys. 2005, 122: 184102-View ArticlePubMedGoogle Scholar
- MacNamara S, Burrage K, Sidje R: Multiscale modeling of chemical kinetics via the master equation. Multiscale Model Simul. 2008, 6 (4): 1146-1168.View ArticleGoogle Scholar
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.