Approximate Numerical Integration of the Chemical Master Equation for Stochastic Reaction Networks

07/24/2019 ∙ by Linar Mikeev, et al. ∙ Universität Saarland 0

Numerical solution of the chemical master equation for stochastic reaction networks typically suffers from the state space explosion problem due to the curse of dimensionality and from stiffness due to multiple time scales. The dimension of the state space equals the number of molecular species involved in the reaction network and the size of the system of differential equations equals the number of states in the corresponding continuous-time Markov chain, which is usually enormously huge and often even infinite. Thus, efficient numerical solution approaches must be able to handle huge, possibly infinite and stiff systems of differential equations efficiently. We present an approximate numerical integration approach that combines a dynamical state space truncation procedure with efficient numerical integration schemes for systems of ordinary differential equations including adaptive step size selection based on local error estimates. The efficiency and accuracy is demonstrated by numerical examples.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The chemical master equation (CME) describes the random dynamics of many stochastic reaction networks in physics, biology, and chemistry, amongst other sciences. More specifically, many biochemical reaction networks can be appropriately modeled by multivariate continuous-time Markov chains (CTMCs), also referred to as Markov jump processes, where at any time the system state is represented by a vector of the numbers of each molecular species present in the reaction network and the CME is a system of differential equations whose solution, given an initial probability distribution, provides the transient (time-dependent) state probabilities.

Exact analytical solutions of the CME are only available for very small reaction networks or special cases such as, e.g., monomolecular reactions or reversible bimolecular reactions [26, 29]. Therefore, in general computational approaches are required. However, because reaction rates typically differ by several orders of magnitude, the system dynamics possess multiple time scales and the corresponding equations are stiff. Furthermore, the size of the state space typically increases exponentially with the model dimensionality, that is, with the number of molecular species in the reaction network. This effect is often referred to as the curse of dimensionality or state space explosion. Stochastic simulation of the reaction network and the numerical solution of the CME are two common complementary approaches to analyze stochastic reaction networks governed by the CME, both of which must be properly designed to cope with huge, potentially infinite multidimensional state spaces, and stiffness.

Stochastic simulation does not solve the CME directly but imitates the reaction network dynamics by generating trajectories (sample paths) of the underlying CTMC, where a stochastically exact imitation is in principle straightforward [3, 16, 17]. Several mathematically equivalent implementations have been developed [10, 15, 30, 34, 48]. Stochastic simulation does not suffer from state space explosion, because the state space needs not be explicitly enumerated, but in particular stochastic simulation of stiff systems becomes exceedingly slow and inefficient, because simulating all successive reactions in order to generate trajectories of the underlying CTMC advances only in extremely small time steps. In order to tackle this problem approximate stochastic simulation techniques for accelerated trajectory generation have been developed, in particular various tau-leaping methods [1, 8, 9, 19, 43, 44, 50, 53, 56] where, instead of simulating every reaction, at any time a time step size is determined by which the simulation is advanced.

A general problem of stochastic simulation, however, remains also with approximate techniques: stochastic simulation constitutes an algorithmic statistical estimation procedure that tends to be computationally expensive and only provides estimates whose reliability and statistical accuracy in terms of relative errors or confidence interval half widths depend on the variance of the corresponding simulation estimator. Estimating the whole probability distribution by stochastic simulation is enormously time consuming. Therefore, often only statistical estimates of the expected numbers of molecules are considered, which requires less effort, but in any case, when estimating expectations, probabilities, or any other relevant system property, many stochastically independent and identically distributed trajectories must be generated in order to achieve a reasonable statistical accuracy

[49]. Hence, stochastic simulation is inherently costly. In many application domains it is sometimes even referred to as a method of last resort.

Several hybrid approximation approaches combine stochastic simulation with deterministic numerical computations. One way to do this is by distinguishing between low molecular counts, where the evolution is described by the CME and handled by stochastic simulation, and high molecular counts, where an approximation of the CME by the continuous-state Fokker-Planck equation is viable [24, 46]

. Another way is to consider time scale separations where parts of the system (states or reactions) are classified as either slow or fast and the different parts are handled by different stochastic simulation approaches and numerical solution techniques

[5, 14, 22, 23, 31, 32, 42, 47, 55], or even without resorting to stochastic simulation for any part [6, 51]. As a major drawback, however, for time scale separation methods it is in general hard to define what is slow or fast. In fact, many systems possess multiple time scales rather than only two and a clear separation of time scales is impossible.

Some numerical approximation approaches tackle the state space explosion problem by restricting the analysis of the model to certain subsets of states where the truncated part of the state space has only a sufficiently small probability. For instance, finite state projection (FSP) algorithms [4, 39, 40] consider finite parts of the state space that can be reached either during the whole time period of interest or during multiple time intervals into which the time period of interest is split. Then the computation of transient probabilities is conducted based on the representation of the transient probability distribution as the product of the initial probability distribution times a matrix exponential involving the generator matrix of the underlying CTMC restricted to the finite projection. However, computing matrix exponentials is well known to be an intricate issue in itself [37, 38] and with FSP algorithms it must be repeated multiple times, corresponding to repeated expansions of the finite state projection. Alternative state space truncation methods are based on adaptive wavelet methods [25, 27], or on a conversion to discrete time where it is dynamically decided which states to consider at a certain time step in a uniformized discrete-time Markov chain [2, 33].

While examples show that these methods can handle the state space explosion in some cases, for many systems they are still not feasible, because even the considered finite part of the state space is too large, or because the specific system dynamics hamper the computational methods involved. In particular, even if the state space is relatively small or can be reduced to a manageable size, the stiffness problem must be handled satisfactorily and it is not clear whether these methods are suitable for stiff systems.

In this paper, we present a novel numerical approximation method that tackles both the state space explosion problem and the stiffness problem where we consider efficient approximate numerical integration based on the fact that the CME can be cast in the form of a system of ordinary differential equations (ODEs) and the computation of transient probabilities, given an initial probability distribution, is an initial value problem (IVP). We combine a dynamical state space truncation procedure, efficient numerical integration schemes, and an adaptive step size selection based on local error estimates. The dynamical state space truncation keeps the number of considered states manageable while incurring only a small approximation error. It is much more flexible and can reduce the state space much more than the aforementioned methods. The use of efficient ODE solvers with adaptive step size control ensures that the method is fast and numerically stable by taking as large as possible steps without degrading the method’s convergence order.

Our method approximates the solution of the CME by truncating large, possibly infinite state spaces dynamically in an iterative fashion. At a particular time instant , we consider an approximation of the transient probability distribution and temporarily neglect states with a probability smaller than a threshold , that is, their probability at time is set to zero. The CME is then solved for an (adaptively chosen) time step during which the truncated state space is adapted to the probability distribution at time . More precisely, certain states that do not belong to the truncated part of the state space at time are added at time , when in the meantime they receive a significant amount of probability which exceeds . Other states whose probabilities drop below are temporarily neglected. The smaller the significance threshold is chosen the more accurate the approximation becomes.

The dynamical state space truncation approach combined with a 4th-order 4-stage explicit Runge-Kutta integration scheme was applied in [35, 36]

to the computation of certain transient rare event probabilities in nonstiff models. The step size at each step was chosen heuristically according to the smallest expected sojourn time in any of the significant states, without any local error control. This works well for nonstiff models, but stiff models require more sophisticated adaptive step size selection strategies based on local error estimates and often even implicit integration schemes

[7, 20, 21, 52].

In the present paper we substantially generalize and improve our previous work. We extend the dynamical state space truncation to work for the whole class of Runge-Kutta methods, explicit and implicit, each of which can be equipped with an adaptive step size selection based on local error estimates. This yields an accurate, numerically stable and computationally efficient framework for the approximate solution of the CME for stochastic reaction networks with extremely huge, possibly infinite, multidimensional state spaces.

In the next section we introduce the notation for stochastic reaction networks, the underlying CTMC and in particular the CME that describes the transient probability distribution. Section 3 describes the dynamical state space truncation, its efficient implementation for general Runge-Kutta methods, and the adaptive step size selection based on local error estimates. Examples are provided in Section 4. Finally, Section 5 concludes the paper and outlines further research directions.

2 Chemical Master Equation

Consider a well-stirred mixture of molecular species in a thermally equilibrated system of fixed volume, interacting through different types of chemical reactions, also referred to as chemical reaction channels, . At any time

a discrete random variable

describes the number of molecules of species and the system state is given by the random vector .

The system changes its state due to one of the possible reactions, where each reaction channel , is defined by a stoichiometric equation


with an associated stochastic reaction rate constant reactants products and the corresponding stoichiometric coefficients where indexes those species that are involved in the reaction. Mathematically, the stoichiometry is described by the state change vector where is the change of molecules of species due to That is, if a reaction of type occurs when the system is in state then the next state is or, equivalently, state is reached, if a reaction of type occurs when the system is in state

For each reaction channel the reaction rate is given by a state-dependent propensity function where is the conditional probability that a reaction of type occurs in the time interval given that the system is in state at time That is


The propensity function is given by times the number of possible combinations of the required reactants and thus computes as


where is the number of molecules of species present in state , and is the stoichiometric coefficient of according to (1). Because at any time the system’s future evolution only depends on the current state, is a time-homogeneous continuous-time Markov chain (CTMC), also referred to as a Markov jump process, with -dimensional state space . Since the state space is countable it is always possible to map it to , which yields a numbering of the states.

The conditional transient (time-dependent) probability that the system is in state at time given that the system starts in an initial state at time is denoted by


and the transient probability distribution at time is the collection of all transient state probabilities at that time, represented by the row vector . The system dynamics in terms of the state probabilities’ time derivatives are described by the chemical master equation (CME)


which is also well known as the system of Kolmogorov forward differential equations for Markov processes [54]. These stochastic reaction kinetics are physically well justified since they are evidently in accordance with the theory of thermodynamics [18, 54]. In the thermodynamic limit the stochastic description converges to classical deterministic mass action kinetics [28].

Note that (5

) is the most common way to write the CME, namely as a partial differential equation (PDE), where

as well as are variables. However, for any fixed state the only free parameter is the time parameter such that (5) with fixed is an ordinary differential equation (ODE) with variable . In particular, when solving for the transient state probabilities numerical ODE solvers can be applied. We shall therefore use the notation in the following.

3 Approximate Numerical Integration of the CME

Numerical integration methods for solving the CME, given , discretize the integration interval and successively compute approximations , where are the mesh points and is the step size at the -th step, . With single-step methods each approximation is computed in terms of the previous approximation only, that is, without using approximations . For advanced methods the step sizes (and thus , the number of steps) are not determined in advance, but variable step sizes are determined in the course of the iteration.

The system of ODEs described by the CME (5) is typically large or even infinite, because there is one ODE for each state in the underlying CTMC, that is, the size of the system of ODEs equals the size of the CTMC’s state space. Thus, its solution with standard numerical integration methods becomes computationally infeasible. However, one can exploit that at any time only a tractable number of states have “significant” probability, that is, only relatively few states have a probability that is greater than a small threshold.

The main idea of our dynamical state space truncation for numerical integration methods is to integrate only those differential equations in the CME (5) that correspond to significant states. All other state probabilities are (temporarily) set to zero. This reduces the computational effort significantly since in each iteration step only a comparatively small subset of states is considered. Based on the fixed probability threshold , we dynamically decide which states to drop or add, respectively. Due to the regular structure of the CTMC the approximation error of the algorithm remains small since probability mass is usually concentrated at certain parts of the state space. The farther away a state is from a “significant set” the smaller is its probability. Thus, in most cases the total error of the approximation remains small. Since in each iteration step probability mass may be “lost” the approximation error at step is the sum of all probability mass lost (provided that the numerical integration could be performed without any errors), that is,


It is important to note that other than static state space truncation approaches our dynamical approach allows that in the course of the computation states can “come and go”, that is, states join the significant set if and only if their current probability is above the threshold and states in the significant set are dropped immediately when their current probability falls below . Furthermore, states that have previously been dropped may come back, that is they are re-considered as significant as soon as they receive a probability that exceeds the threshold . This is substantially different from state space exploration techniques where only the most probable states are generated but states are never dropped as time progresses like for instance in [12] with regard to approximating stationary distributions. Our dynamical state space truncation approach is also much more flexible than finite state projection (FSP) algorithms [4, 39, 40] which work over pre-defined time intervals with the same subset of states, where in particular for stiff systems many reactions can occur during any time interval, so that in order to safely meet reasonable accuracy requirements the resulting subset of states is often still extremely large. In contrast to that we update our set of significant states in each adaptively chosen time step, without much overhead. Furthermore, by numerically integrating the ODEs we avoid the intricate computation of matrix exponentials required in FSP algorithms and by using an efficient data structure we do not even need to generate any matrices.

In order to avoid the explicit construction of a matrix and in order to work with a dynamic set of significant states that changes in each step, we use for a state a data structure with the following components:

  • fields , for the approximated probabilities and , respectively,

  • for all with a pointer to the successor state as well as a field with the rate .

2: 0, 0
4:while  do
5:     compute
6:     compute ,
7:     if step successful then
8:         update
9:          , i i+1
10:     end if
11:end while
Algorithm 1 A numerical adaptive step size integration scheme.

The workflow of the numerical integration scheme is given in pseudocode in Algorithm 1. We start at time and initialize the set as the set of all states that have initially a probability greater than . We compute the initial time step . In each iteration step we compute the approximation using an explicit or implicit Runge-Kutta method (see Sections 3.1, 3.2). We check whether the iteration step was successful computing the local error estimate and ensuring that error tolerance conditions are met. If so, then for each state we update the field with , and remove the state from if its probability becomes less than . Based on the local error estimate, we choose a time step for the next iteration (or the repetition of the iteration in case it was not successful). This and the computation of the initial time step is detailed in Section 3.3.

3.1 Runge-Kutta methods

We consider the whole family of Runge-Kutta methods, which proceed each time step of given step size in multiple stages. More precisely, a general -stage Runge-Kutta method proceeds according to the iteration scheme


which is uniquely defined by weights with and coefficients , . Thus, is an approximation to , where , and is the component of the vector that corresponds to state . Hence, is the probability of at stage . If for all , then the sum in Equation (8) effectively runs only from to , which means that for each the computation of includes only previous stage terms . Therefore, can be computed sequentially, that is, for all yields explicit integration schemes. If there is at least one with , then the integration scheme is implicit, which implies that the solution of at least one linear system of equations is required per iteration step.

1:for  1 to  do
2:     for all  do
4:         for all  :                 do
5:              if  then  end if
8:         end for
9:     end for
10:end for
11:for all  do
13:      0, …, 0
14:end for
Algorithm 2 A single iteration step of a general explicit -stage Runge-Kutta scheme, defined by , and for , with dynamical state space truncation, which approximates the solution of the CME.

A single iteration step for general explicit -stage Runge-Kutta schemes is given in pseudocode in Algorithm 2. Note that for each state besides , , we additionally store fields for the stage terms and initialize them with zero. We compute the approximation of based on Equation (7) by traversing the set times. In the first rounds (lines 1-10) we compute according to Equation (8) and in the final round (lines 11-14) we accumulate the summands and zero . While processing state in round , for each reaction channel , we transfer probability mass from state to its successor , by subtracting a term from and adding the same term to (lines 6-7). We do so after checking (line 4) whether is already in , and if not, whether it is worthwhile to add to , that is, we guarantee that will receive enough probability mass and that will not be removed in the same iteration. Thus, we add to (line 5) only if the inflow to is greater or equal than a certain threshold . Obviously, may receive more probability mass from other states and the total inflow may be greater than . Thus, if a state is not a member of and if for each incoming transition the inflow probability is less than , then this state will not be added to even if the total inflow is greater or equal than . This small modification yields a significant speed-up since otherwise all states that are reachable within at most transitions will always be added to , but many of the newly added states will be removed in the same iteration.

3.2 Implicit Methods

The advantage of implicit methods is that they can usually take larger (and thus fewer) steps, which comes at the price of an increased computational effort per step, but paying that price can lead to large speed-up of the overall integration over the time interval . It is common sense that in general the efficient solution of stiff ordinary differential equations requires implicit integration schemes [21]. In the present paper, with regard to implicit numerical integration we restrict ourselves to the implicit Euler method, hence the special case of a one-stage Runge-Kutta method with and , which yields


and requires to solve the linear system


for in each step.

Of course, when considering a standard approach to the numerical integration of the CME, where no state space truncation is considered, then this linear system is huge, possibly infinite, and its solution is often impossible. In conjunction with our dynamical state truncation procedure, the linear system is reduced similarly to the reduction of the state space and the number of differential equations to be integrated per step, respectively. However, there a subtleties that must be properly taken into account.

Firstly, since we do not need to maintain huge matrices but we use the previously described dynamical data structure the solution of the linear system must be accordingly implemented with this data structure. Secondly, some states that are not significant at time may receive a significant probability at time and must be included in the linear system. Thus, a dynamical implementation of the solution of the linear system is required. Therefore, iterative solution techniques for linear systems that are usually simply defined by a fixed matrix must be properly adapted to the dynamical data structure and the dynamical state space truncation.

In fact, when using implicit numerical integration schemes in conjunction with the dynamical state space truncation procedure, in principle the solution of a linear system in each integration step is a challenging potential bottleneck. It is therefore a key point and a key contribution to implement it efficiently.

1:while convergence not reached do
2:     for all  do
3:         for all  :                 do
4:              if  then  end if
6:         end for
7:     end for
8:     for all  do
10:         check convergence for state
12:          0
13:     end for
14:end while
Algorithm 3 A single iteration step of an implicit Euler scheme using the Jacobi method with dynamical state space truncation, which approximates the solution of the CME.

We illustrate the solution of the linear system (10) using the Jacobi method, which yields the following iterative scheme


where . The pseudocode is given in Algorithm 3. In the -th iteration of the adaptive numerical integration scheme (Algorithm 1), we store the “old” approximation of the state probability in the field and the “new” approximation in the field . We initialize with the state probability from the -th iteration . In lines 2-7 for each state we compute the sum and store it in a field . While processing state , for each reaction channel , we transfer probability mass from state to its successor . Similarly to Algorithm 2, we only add a new state to if it receives enough probability mass. In lines 9-14 for each state we compute the “new” approximation of according to (11) and check whether the convergence criterion


is fulfilled for some relative and absolute tolerances and . Algorithm 3 terminates if (12) holds for all states . After the convergence of the Jacobi method, the field contains the approximation .

For our numerical experiments we use the Gauss-Seidel method, which is known to converge faster than the Jacobi method. The iterative solution is given as


where is an indicator (or flag) that takes the value if the state has been already processed, and otherwise. Thus, in (13) in the summation, we use the “new” approximations of the processed states and the “old” probability if the approximation was not yet updated in the current iteration. We modify Algorithm 3 as follows. We compute the sum as before, but after processing a state in line 9, we mark it as processed and propagate to the successor states which are not marked as processed. After this, the field contains the sum required for (13).

3.3 Local Error Control and Adaptive Step Size Selection

The accuracy as well as the computing time of numerical integration methods depend on the order of the method and the step size. The error in a single step with step size is approximately with a factor that varies over the integration interval. Hence, one crucial point for the efficiency of numerical integration methods is the step size selection. It is well known that methods with constant step size perform poorly if the solution varies rapidly in some parts of the integration interval and slowly in other parts [7, 20, 21, 52]. Therefore, adaptive step size selection so that the accuracy and the computing time are well balanced is highly desirable for explicit and implicit integration schemes. For both classes of schemes we base our step size selection strategy on local error estimates.

Our goal is to control the local error and, accordingly, to choose the step size so that at each step for all states ,


where RTOL and ATOL are user-specified relative and absolute error tolerances. In particular, note that we use a mixed error control, that is, a criterion that accounts for both the relative and the absolute error via corresponding relative and absolute error tolerances, because in practice using either a pure relative error control or a pure absolute error control can cause serious problems, see, e.g., Section 1.4 of [52]. Of course, the true local errors are not available and we must estimate them along with each integration step.

For the explicit and the implicit Euler method we compute a local error estimate similarly to the step doubling approach, that is, we approximate by taking the time step and independently taking two consecutive time steps of length . The local error estimate is then the vector


with components , where and denote the approximations computed with time step and with two consecutive time steps of length , respectively.

The embedded Runge-Kutta methods provide an alternative way for the step size control. Along with the approximation of order , they deliver the approximation of order computed as


where with . Then the local error estimate is the vector


with components .

Now, with regard to the step size selection assume we have made a ‘trial’ step with a given step size and computed the corresponding local error estimate. Then we accept the step if for all significant states the local error estimate is smaller than the prescribed local error tolerance. More precisely,


which implies


If the step is not accepted, then we have to decrease the step size. Otherwise we can proceed to the next step where it is likely that we can use an increased step size, because the current one might be smaller than necessary. In both cases, acceptance or rejection, we have to specify by how much the step size is decreased or increased, respectively, and in both cases we do this based on the local error estimate.

Define , denote by the corresponding vector containing the components , and define


It can be easily seen that the largest step size that yields a local error estimate satisfying (18) can be approximated by . Note that if (18) is satisfied and otherwise. This means we can use as a factor in both possible cases, that is, for a too large step size that has been rejected and must be decreased for a retrial, and for an accepted step size to set an increased step size for the next step. In practice we also have to account for the fact that the local error is only estimated. Rejecting a step and re-computing it with a smaller step size should be avoided as much as possible. Therefore, rather than we consider with a safety factor . Besides, step sizes must not be too large and also too large changes of the step size must be avoided since otherwise the above approximation of the largest possible step size is not valid [52]. If the step with step size is accepted, then we set


as the initial trial step size for the next step, where . Otherwise, we decrease the step size for the current step according to


where and is the absolute distance between and the next floating-point number of the same precision as . So, is such that and are different in working precision.

For the computation of the initial trial step size, we first compute (see (5)). This can be done using one stage of Algorithm 2. Then we compute


In the -th iteration, we stretch the time step , if it lies within 10% of . Thus, we set the time step to if . Note that this also covers the case when the time step is too large, and using it would lead to jumping over the final time point .

4 Numerical examples

In this section, we present numerical examples in order to demonstrate the suitability of our approach, its accuracy, run time and the number of significant states to be processed corresponding to the number of differential equations to be integrated. As our first example we consider a birth-death process for which analytical solutions are available such that we can indeed compare our numerical results with exact values. Then we consider a more complex yeast cell polarization model.

We compare the accuracies, run times and numbers of significant states of explicit Euler (referred to as ‘euler’ in the following figures and tables), implicit (backward) Euler (‘beuler’), and an embedded Runge-Kutta (‘rk45’) with weights and coefficients chosen according to [13] together with local error control and adaptive step size selection as described in the previous section. Note that rk45 is similar to the ode45 method of MATLAB, but while with MATLAB’s ode45 only systems of ODEs of moderate size can be solved, here, of course, we consider it in conjunction with our dynamical state space truncation procedure.

For our numerical experiments we fix the relative tolerance . For the dynamical state space truncation, we use , which agrees with the error control property of the ODE solution that the components smaller than ATOL are unimportant. As a safety factor in the time step selection procedure we use . In the solution of the linear system required for the implicit Euler method we set and .

4.1 Birth-Death Process

Our first example is the birth-death process given as

with as the only species and propensity functions . It is clear that the state space is the infinite set of all nonnegative integers so that the corresponding system of differential equations is infinite, too. We chose the rate constants , , the initial state and final time horizon . We analyze the model with different values of . Since reporting the probabilities of single states over time is not of any practical interest we focus on representative properties and on informative measures of the accuracy and the efficiency of our method.

In Figure 1 we plot the average number of species over time as obtained with our approximate numerical integration schemes, where the run times were less than one second. We also plot the exact solution obtained according to [26]. The plots show for all considered values of ATOL that there is no visible difference between the exact values and our approximations, which suggests that our approximations are indeed extremely accurate.

Figure 1: Average species count for birth-death process ( (a) , (b) , (c) ).
Figure 2: L2 error for birth-death model ( (a) , (b) , (c) ).

Figure 2 depicts plots of the L2 error

with . The norm is computed over all states with positive probabilities in the exact solution. Note that if there is no corresponding state in , its probability is taken as 0.

It can be seen that in all cases the L2 error is less than , which confirms the high accuracy of the approximations also formally. It can be further seen that for this example euler is slightly more accurate than beuler and that rk45 is even by orders of magnitude more accurate than the Euler schemes. This is well in accordance with the higher order of rk45.

Figure 3 shows the average step sizes taken by the different integration schemes, where rk45 takes much larger steps than the Euler schemes.

Figure 3: Average step sizes for birth-death process ( (a) , (b) , (c) ).

In Figure 4 we plot the numbers of significant states used during the computation. It can be seen that for all considered values of ATOL all integration schemes only require to handle (integrate) a moderate number of states (differential equations) where the Euler schemes require roughly the same number of states and rk45 requires only slightly more, in any case for any time less than states. Of course, the smaller ATOL (and thus our truncation probability threshold ) the larger the number of significant states but with only a slight increase. This shows that the dynamical truncation procedure indeed substantially reduces the size of the state space and thus renders possible to integrate numerically with – as demonstrated by the previous figures – maintaining a high accuracy of the approximations.

Figure 4: Numbers of significant states for birth-death process ( (a) , (b) , (c) ).

4.2 Yeast cell polarization

As another reference example we consider a stochastic model of the pheromone-induced G-protein cycle in the yeast Saccharomyces cerevisiae [11, 41, 45]

with initial state , where the state vector is given as and the state space is the infinite -dimensional set . Note that in order to keep the meaning of the species here we do not number the species but take the notation from [45].

In Figure 5 we plot the average species counts computed using rk45 with and the numbers of significant states for rk45 over time for different values of .

Figure 5: Yeast cell polarization: (a,b) average species counts, (c) number of significant states for different values of ATOL.

Note that as in the previous example the numbers of significant states as well as the average numbers of species for euler and beuler only slightly differ from those for rk45, so that we omit to include them in the plots. It is clear that in the much more complex yeast cell polarization model there are many more significant states than in the birth-death process, but the numbers of significant states are still in a range that allows accurate approximations in reasonable time.

In Table 1 we list the run times for different values of ATOL. We can see that in this example beuler heavily outperforms the explicit methods euler and rk45. This confirms the advantages of implicit methods over explicit methods for stiff systems. In fact, while the birth-death example is not or only moderately stiff, the yeast cell polarization model constitutes a very stiff system of differential equations.

euler 761s 1481s 5839s
beuler 35s 76s 139s
rk45 5806s 18603s 26126s
Table 1: Run times for yeast cell polarization model.

In Figure 6 we plot the average step size over time and the L2 error for . Since there is no analytical solution, we compute the L2 with respect to the distribution obtained using rk45 with a lower absolute tolerance , applying the rationale that with an even lower error tolerance the method gives nearly exact values. The respective L2 errors for and are similar.

Figure 6: (a) Average step sizes and (b) L2 error for yeast polarization model ().
Figure 7: Relative errors.

5 Conclusion

We have shown that our numerical integration approach with adaptive step size selection based on local error estimates performed in combination with dynamical state space truncation provides a versatile means of approximating the solution of the chemical master equation for complex stochastic reaction networks efficiently and accurately. The state space explosion problem is circumvented by considering in each time step only such states of the overall state space of the reaction network that have at that time step a significant (sufficiently large according to a flexibly adjustable bound) probability, that is, in the course of the integration scheme we keep the number of differential equations to be integrated per step manageable. By a framework that includes explicit as well as implicit integration schemes we offer the flexibility to choose an appropriate integration scheme that is well-suited with regard to the specific dynamics of a given reaction network.

In order to provide meaningful, detailed comparisons of different methods with different parameter choices and to study their impact on accuracy, run times and numbers of significant states to be processed we have considered the explicit Euler method, an explicit Runge-Kutta as an extension of explicit Euler, and the implicit (backward) Euler method, all equipped with a well-suited adaptive step size selection strategy and performed on the dynamically truncated state space. The results show that the proposed approximate numerical integration of the chemical master equation indeed yields satisfactorily accurate results in reasonable time. Future research will be concerned with further advanced integration schemes, to equip them similarly with adaptive step size selection strategies and to study the accuracy and the run times. Of course, also in-depth theoretical investigations of the performance, efficiency and accuracy of the approximate numerical integration approach are highly desirable.


  • [1] D. F. Anderson. Incorporating postleap checks in tau-leaping. Journal of Chemical Physics, 128(5):054103, 2008.
  • [2] A. Andreychenko, W. Sandmann, and V. Wolf. Approximate adative uniformization of continuous-time Markov chains. Applied Mathematical Modelling, 61:561–576, 2018.
  • [3] A. B. Bortz, M. H. Kalos, and J. L. Lebowitz. A new algorithm for Monte Carlo simulation of Ising spin systems. Journal of Computational Physics, 17:10–18, 1975.
  • [4] K. Burrage, M. Hegland, S. MacNamara, and R. Sidje. A Krylov-based finite state projection algorithm for solving the chemical master equation arising in the discrete modelling of biological systems. In A. N. Langville and W. J. Stewart, editors, Proceedings of the Markov Anniversary Meeting: An International Conference to Celebrate the 150th Anniversary of the Birth of A. A. Markov, pages 21–38, 2006.
  • [5] K. Burrage, T. Tian, and P. Burrage. A multi-scaled approach for simulating chemical reaction systems. Progress in Biophysics and Molecular Biology, 85:217–234, 2004.
  • [6] H. Busch, W. Sandmann, and V. Wolf. A numerical aggregation algorithm for the enzyme-catalyzed substrate conversion. In Proceedings of the 2006 International Conference on Computational Methods in Systems Biology, volume 4210 of Lecture Notes in Computer Science, pages 298–311. Springer, 2006.
  • [7] J. C. Butcher. Numerical Methods for Ordinary Differential Equations. John Wiley & Sons, 2nd edition, 2008.
  • [8] Y. Cao, D. T. Gillespie, and L. R. Petzold. Efficient stepsize selection for the tau-leaping simulation. Journal of Chemical Physics, 124:044109–144109–11, 2006.
  • [9] Y. Cao, D. T. Gillespie, and L. R. Petzold. The adaptive explicit-implicit tau-leaping method with automatic tau selection. Journal of Chemical Physics, 126(22):224101, 2007.
  • [10] Y. Cao, H. Li, and L. R. Petzold. Efficient formulation of the stochastic simulation algorithm for chemically reacting systems. Journal of Chemical Physics, 121(9):4059–4067, 2004.
  • [11] C.-S. Chou, Q. Nie, and T.-M. Yi. Modeling robustness tradeoffs in yeast cell polarization induced by spatial gradients. PLoS ONE, 3(9):e3103, 2008.
  • [12] E. de Souza e Silva and P. M. Ochoa.

    State space exploration in Markov models.

    In Proceedings of the ACM SIGMETRICS Joint International Conference on Measurement and Modeling of Computer Systems, pages 152–166, 1992.
  • [13] J. R. Dormand and P. J. Prince. A family of embedded Runge-Kutta formulae. Journal of Computational and Applied Mathematics, 6(1):19–26, 1980.
  • [14] W. E., D. Liu, and E. Vanden-Eijnden. Nested stochastic simulation algorithms for chemical kinetic systems with multiple time scales. Journal of Computational Physics, 221:158–180, 2007.
  • [15] M. A. Gibson and J. Bruck. Efficient exact stochastic simulation of chemical systems with many species and many channels. Journal of Physical Chemistry, A 104:1876–1889, 2000.
  • [16] D. T. Gillespie. A general method for numerically simulating the time evolution of coupled chemical reactions. Journal of Computational Physics, 22:403–434, 1976.
  • [17] D. T. Gillespie. Exact stochastic simulation of coupled chemical reactions. Journal of Physical Chemistry, 71(25):2340–2361, 1977.
  • [18] D. T. Gillespie. A rigorous derivation of the chemical master equation. Physica A, 188:404–425, 1992.
  • [19] D. T. Gillespie. Approximate accelerated stochastic simulation of chemically reacting systems. Journal of Chemical Physics, 115:1716–1732, 2001.
  • [20] E. Hairer, S. P. Nørsett, and G. Wanner. Solving Ordinary Differential Equations I: Nonstiff Problems. Springer-Verlag, Berlin Heidelberg, 2nd edition, 1993.
  • [21] E. Hairer and G. Wanner. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Springer-Verlag, Berlin Heidelberg, 2nd edition, 1996.
  • [22] L. A. Harris and P. Clancy. A partitioned leaping approach for multiscale modeling of chemical reaction dynamics. Journal of Chemical Physics, 125(14):144107, 2006.
  • [23] E. L. Haseltine and J. B. Rawlings. Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics. Journal of Chemical Physics, 117(15):6959–6969, 2002.
  • [24] A. Hellander and P. Lötstedt. Hybrid method for the chemical master equation. Journal of Computational Physics, 227:100–122, 2007.
  • [25] T. Jahnke. An adaptive wavelet method for the chemical master equation. SIAM Journal on Scientific Computing, 31(6):4373 –4394, 2010.
  • [26] T. Jahnke and W. Huisinga. Solving the chemical master equation for monomolecular reactions analytically. Journal of Mathematical Biology, 54(1):1–26, 2007.
  • [27] T. Jahnke and T. Udrescu. Solving chemical master equations by adaptive wavelet compression. Journal of Computational Physics, 229:5724–5741, 2010.
  • [28] T. G. Kurtz. The relationship between stochastic and deterministic models for chemical reactions. Journal of Chemical Physics, 57(7):2976–2978, 1972.
  • [29] I. J. Laurenzi. An analytical solution of the stochastic master equation for reversible bimolecular reaction kinetics. Journal of Chemical Physics, 113(8):3315–3322, 2000.
  • [30] H. Li and L. R. Petzold. Logarithmic direct method for discrete stochastic simulation of chemically reacting systems. Technical Report, Department of Computer Science, University of California, Santa Barbara, 2006.
  • [31] S. MacNamara, A. M. Bersani, K. Burrage, and R. B. Sidje. Stochastic chemical kinetics and the total quasi-steady-state assumption: application to the stochastic simulation algorithm and chemical master equation. Journal of Chemical Physics, 129:095105, 2008.
  • [32] S. MacNamara, K. Burrage, and R. B. Sidje. Multiscale modeling of chemical kinetics via the master equation. Multiscale Modeling & Simulation, 6(4):1146–1168, 2008.
  • [33] M. Mateescu, V. Wolf, F. Didier, and T.A. Henzinger. Fast adaptive uniformisation of the chemical master equation. IET Systems Biology, 4(6):441–452, 2010.
  • [34] J. M. McCollum, G. D. Peterson, C. D. Cox, M. L. Simpson, and N. F. Samatova. The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior. Computational Biology and Chemistry, 30(1):39–49, 2006.
  • [35] L. Mikeev, W. Sandmann, and V. Wolf. Efficient calculation of rare event probabilities in Markovian queueing networks. In Proceedings of the 5th International Conference on Performance Evaluation Methodologies and Tools, VALUETOOLS, pages 186–196, 2011.
  • [36] L. Mikeev, W. Sandmann, and V. Wolf. Numerical approximation of rare event probabilities in biochemically reacting systems. In A. Gupta and T. A. Henzinger, editors, Proceedings of the 11th International Conference on Computational Methods in Systems Biology, CMSB, volume 8130 of Lecture Notes in Bioinformatics, pages 5–18. Springer, 2013.
  • [37] C. Moler and C. Van Loan. Nineteen dubious ways to compute the exponential of a matrix. SIAM Review, 20:801–836, 1978.
  • [38] C. Moler and C. Van Loan. Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM Review, 45:3–49, 2003.
  • [39] B. Munsky and M. Khammash. The finite state projection algorithm for the solution of the chemical master equation. Journal of Chemical Physics, 124:044104, 2006.
  • [40] B. Munsky and M. Khammash. A multiple time interval finite state projection algorithm for the solution to the chemical master equation. Journal of Computational Physics, 226:818–835, 2007.
  • [41] D. Pruyne and A. Bretcher. Polarization of cell growth in yeast I. Establishment and maintenance of polarity states. Journal of Cell Science, 113:365–375, 2000.
  • [42] C. V. Rao and A. Arkin. Stochastic chemical kinetics and the quasi-steady-state-assumption: Application to the Gillespie algorithm. Journal of Chemical Physics, 118:4999–5010, 2003.
  • [43] M. Rathinam and H. El Samad. Reversible-equivalent-monomolecular tau: A leaping method for small number and stiff stochastic chemical systems. Journal of Computational Physics, 224:897–923, 2007.
  • [44] M. Rathinam, L. R. Petzold, Y. Cao, and D. T. Gillespie. Stiffness in stochastic chemically reacting systems: The implicit tau-leaping method. Journal of Chemical Physics, 119:12784–12794, 2003.
  • [45] M. K. Roh, B. J. Daigle Jr., D. T. Gillespie, and L. R. Petzold. State-dependent doubly weighted stochastic simulation algorithm for automatic characterization of stochastic biochemical rare events. Journal of Chemical Physics, 135:234108, 2011.
  • [46] C. Safta, K. Sargsyan, B. Debusschere, and H. N. Najm. Hybrid discrete/continuum algorithms for stochastic reaction networks. Journal of Computational Physics, 281:177–198, 2015.
  • [47] H. Salis and Y. Kaznessis. Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions. Journal of Chemical Physics, 122:054103, 2005.
  • [48] W. Sandmann. Discrete-time stochastic modeling and simulation of biochemical networks. Computational Biology and Chemistry, 32(4):292–297, 2008.
  • [49] W. Sandmann. Sequential estimation for prescribed statistical accuracy in stochastic simulation of biological systems. Mathematical Biosciences, 221(1):43–53, 2009.
  • [50] W. Sandmann. Streamlined formulation of adaptive explicit-implicit tau-leaping with automatic tau selection. In Proceedings of the 2009 Winter Simulation Conference, pages 1104–1112, 2009.
  • [51] W. Sandmann and V. Wolf. Computational probability for systems biology. In Proceedings of the Workshop on Formal Methods in Systems Biology, volume 5054 of Lecture Notes in Computer Science, pages 33–47. Springer, 2008.
  • [52] L. F. Shampine, I. Gladwell, and S. Thompson. Solving ODEs with MATLAB. Cambridge University Press, 2003.
  • [53] T. Tian and K. Burrage. Binomial leap methods for simulating stochastic chemical kinetics. Journal of Chemical Physics, 121:10356–10364, 2004.
  • [54] N. van Kampen. Stochastic Processes in Physics and Chemistry. Elsevier, North-Holland, 1992.
  • [55] K. Vasudeva and U. S. Bhalla. Adaptive stochastic-deterministic chemical kinetics simulation. Bioinformatics, 20(1):78–84, 2004.
  • [56] Z. Xu and W. Cai. Unbiased tau-leaping methods for stochastic simulation of chemically reacting systems. Journal of Chemical Physics, 128(15):154112, 2008.