Perturbation and Truncation of Probability Generating Function Methods for Stiff Chemical Reactions

Abstract

One can reformulate chemical master equations of the stochastic reaction network into a partial differential equation (PDE) of a probability generating function (PGF). In this paper, we present two improvements in such PGF-PDE approach, based on perturbation and double-truncation, respectively. The stiff system that involves fast and slow reactions together often requires high computational cost. By applying the perturbation method to PGF-PDEs, we expand the equation in terms of a small reaction rate which is often responsible for such stiffness of the system. Also by doubly truncating, we dump relatively small terms and reduce the computational load significantly at each time step. The terms corresponding to rare events are sieved out through truncations of Taylor expansion. It is shown through numerical examples of enzyme kinetics, transition model, and Brusselator model that the suggested method is accurate and efficient for approximation of the state probabilities.

Academic Editor:Lotfollah Najjar

Department of Mathematical Sciences, School of Natural Science, Ulsan National Institute of Science and Technology (UNIST), Ulsan Metropolitan City 689-798, Republic of Korea

Received 15 September 2014; Revised 3 February 2015; Accepted 6 February 2015; 28 May 2015

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Stochastic chemical kinetics arises mostly in the modeling of chemical reactions in relatively small biochemicalbhlt systems such as gene regulatory networks, cell systems, and protein folding models [1]. The governing equation, chemical master equation (CME), describes the probability distribution of the system under the Markov assumption as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , each [figure omitted; refer to PDF] denotes the molecular number of [figure omitted; refer to PDF] th species at time [figure omitted; refer to PDF] , [figure omitted; refer to PDF] is the propensity function or probability intensity for the [figure omitted; refer to PDF] th reaction, and [figure omitted; refer to PDF] denotes the [figure omitted; refer to PDF] th column of the stoichiometric matrix [figure omitted; refer to PDF] of which [figure omitted; refer to PDF] th entry is the change in the number of molecules of the [figure omitted; refer to PDF] th species by the occurrence of the [figure omitted; refer to PDF] th reaction [2].

Since (1) is a high dimensional system in most real applications, it is very challenging to handle both analytically and computationally. One of the key methods to help simplify the computation of the stochastic dynamics is done by using the stochastic simulation algorithm (SSA) [3]. One realization of the SSA shows the trajectory of the time-evolution of the states. Since it is basically a Monte-Carlo type algorithm, one has to perform many realizations to obtain important statistical quantities such as mean and variance. Moreover, it often requires a large number of realizations, especially in cases of the systems consisting of large numbers of molecules or very different time scale reactions. There have been many approaches to improve the SSA including Tau-leaping methods [4-6], reduction methods on the slow time scale [7-10], and moment closure methods [11-13].

In the recent papers [14, 15], the authors presented numerical schemes based on probability generating functions (PGF) for finding the probability distribution and moments. One can reformulate chemical master equations of the stochastic reaction network into a partial differential equation (PDE). Although such PDEs are mostly hard to deal with due to variable coefficients and lack of proper boundary conditions, several improvements in approximation of the moments based on Pade approximation and numerical superimposition have been made.

In this paper, we present two improvements in such PGF-PDE approach. The first method is an approximation method called "perturbed probability generating function (PPGF) method." By applying the perturbation method to PDEs, we expand the equation in terms of a small reaction rate which is often responsible for the stiffness of the system. This enables us to handle the system more efficiently with less iteration and therefore capture same accuracy faster.

We also develop the double truncation method (DTM) that efficiently approximates the state probabilities. As suggested from the name, the method combines two truncation procedures for the series solution to PGF-PDEs.

Through consecutive truncations of the Taylor expansion of such PDEs in time, one can dump relatively rare events which correspond to excessive occurrence of reactions during a short time interval. This elimination is justified from the fact that the coefficient of each term in the PGF is the probability of a specific state. This enables us to ignore a great deal of small terms which do not affect the system significantly.

An outline of the paper is as follows. In Section 2, we review the properties of PDEs that PGF should satisfy. The perturbation of PGF is discussed in Section 3 with its numerical accuracy of through some examples. Section 4 presents the double truncation method for the PGF-PDEs. In Sections 3 and 4, we also illustrate numerical accuracy of the method in several examples.

2. PGF-PDE Approach

The probability generating function is defined as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] After differentiation of (2) with respect to [figure omitted; refer to PDF] and application of (1), one can derive a PDE [figure omitted; refer to PDF] Here [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] is the highest order of the reactions in the chemical system, and [figure omitted; refer to PDF] denotes any [figure omitted; refer to PDF] th order partial derivatives of [figure omitted; refer to PDF] as [figure omitted; refer to PDF] One can see that the PDE (4) is linear and [figure omitted; refer to PDF] is a polynomial function of degree at most order [figure omitted; refer to PDF] which satisfies [figure omitted; refer to PDF] for [figure omitted; refer to PDF] .

The conditions of the PDE (4) can be found from the conditions on [figure omitted; refer to PDF] [14]. The initial condition is [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the initial condition of [figure omitted; refer to PDF] , and also it is true that [figure omitted; refer to PDF]

Using the solution [figure omitted; refer to PDF] of the PDE (4), one can find the important information for stochastic reaction network such as marginal probability, mean, and the second moment [14]. If [figure omitted; refer to PDF] denotes the number of molecules of [figure omitted; refer to PDF] th species at time [figure omitted; refer to PDF] , the probability that [figure omitted; refer to PDF] at time [figure omitted; refer to PDF] is [figure omitted; refer to PDF] and the mean and second moment are obtained by the first and second derivatives of [figure omitted; refer to PDF] [figure omitted; refer to PDF] where [figure omitted; refer to PDF] denotes the expectation of [figure omitted; refer to PDF] .

3. Perturbation of Probability Generating Function Method

3.1. Perturbation of PGF

For the PGF method recently proposed in [14], [figure omitted; refer to PDF] is represented by a power expansion with respect to [figure omitted; refer to PDF] as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , are coefficient functions of [figure omitted; refer to PDF] . Note that the initial condition determines the first coefficient as [figure omitted; refer to PDF] .

Suppose that the reaction rates in the network do not change with time and therefore the coefficients of the PDE do not involve [figure omitted; refer to PDF] explicitly. This implies [figure omitted; refer to PDF] in (4). Plugging (10) into the both sides of (4) and comparing the power series lead to a recursive relation between [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , [figure omitted; refer to PDF] . One can derive coefficient functions [figure omitted; refer to PDF] from [figure omitted; refer to PDF] and then construct an approximation [figure omitted; refer to PDF] for sufficiently large [figure omitted; refer to PDF] . Since the initial function [figure omitted; refer to PDF] and the coefficient of the PDE (2) are polynomials, all [figure omitted; refer to PDF] generated are polynomials as well.

Now, let us apply the perturbation method to PGF. Since many practical chemical networks that interest us involve reactions in a broad range of time scales, from seconds to microseconds, it is important to efficiently handle such stiffness in simulation of the systems. Suppose that [figure omitted; refer to PDF] is a small reaction rate in the system; then one can expand the PGF as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , are functions to be determined. As in (10), one can set for each [figure omitted; refer to PDF] [figure omitted; refer to PDF] Now, instead of (10), one can seek [figure omitted; refer to PDF] in the following form: [figure omitted; refer to PDF] for large [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . By plugging (13) into the both sides of (4) and comparing the power series, one can find recurrence relations for [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , where [figure omitted; refer to PDF] . This decomposition enables us to handle stiffness caused by [figure omitted; refer to PDF] properly, by adjusting the perturbation order [figure omitted; refer to PDF] as will be shown in the numerical examples in the next sections.

3.2. Example 1: Enzyme Kinetics

As the first application, we consider a fundamental chemical kinetic model, the enzyme-substrate reaction system [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] denote enzyme, substrate, enzyme-substrate complex, and product, respectively, and [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] are probability constants for reactions. If we denote the molecular numbers of [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] by [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , respectively, the governing equation of the stochastic enzyme-substrate system is obtained [16] as [figure omitted; refer to PDF]

One can see that there are two conservation quantities in the system and thus we can remove two of variables, say [figure omitted; refer to PDF] , [figure omitted; refer to PDF] . After the removal, we can derive a PDE of [figure omitted; refer to PDF] from the master equation [14] [figure omitted; refer to PDF] Here [figure omitted; refer to PDF] is the conserved quantity [figure omitted; refer to PDF] and we use the notations [figure omitted; refer to PDF]

In the example, we use parameters [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] and let the perturbation be expanded in [figure omitted; refer to PDF] . Figure 1 shows comparison of the exact solution of the master equation (15) and the approximate solution obtained from the PPGF method.

Figure 1: Comparison of mean, standard deviation and marginal probability for the enzyme-substrate model. In the three figures, we assume that [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . (a) shows the comparison of exact values and approximate values obtained from the PPGF method for mean and mean [figure omitted; refer to PDF] standard deviations ( [figure omitted; refer to PDF] ) of [figure omitted; refer to PDF] and (b) shows the comparison of the exact probability and the approximate probability obtained from the PPGF method. In the figure, [figure omitted; refer to PDF] denotes the probability that the number [figure omitted; refer to PDF] of the substrate is [figure omitted; refer to PDF] . The exact probability, mean, and variance are obtained by rewriting (15) into the linear ODE system [figure omitted; refer to PDF] under the given initial condition and then solving it by MATLAB.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

3.3. Example 2: Brusselator Model

We consider the Brusselator model: [figure omitted; refer to PDF] Here [figure omitted; refer to PDF] , [figure omitted; refer to PDF] are molecular species and [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] are constant species. The Brusselator is a generic model for a type of autocatalytic reaction. It is well known that the model reveals simple chemical oscillations. Refer to Figure 2(a) for a typical run of the model. If we denote the number of [figure omitted; refer to PDF] , [figure omitted; refer to PDF] by [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , respectively, we obtain the PDE for the probability generating function as follows: [figure omitted; refer to PDF] Here [figure omitted; refer to PDF] and the reaction rates [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] are used. We use PPGF with respect to the smallest parameter [figure omitted; refer to PDF] . Figure 2 compares the corresponding numerical result to those from SSA.

Figure 2: Comparison of mean and standard deviation for the Brusselator model. The initial condition [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are assumed. (a) shows a typical run of [figure omitted; refer to PDF] in 500 seconds from a Monte Carlo simulation of the Brusselator model. In (b), [figure omitted; refer to PDF] denotes the probability that the number [figure omitted; refer to PDF] of the substrate is [figure omitted; refer to PDF] . (c) and (d) compare mean and mean ± standard deviations ( [figure omitted; refer to PDF] ) of [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] obtained from the SSA and the PPGF. The results from the SSA are based on 50,000 realizations.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

4. Double Truncation Method

4.1. Truncations in [figure omitted; refer to PDF] and [figure omitted; refer to PDF]

The truncation method was firstly introduced in [15] to deal with PGF-PDEs. Although the method provides a fast way to evaluate the mean and the variance, most of high order terms need to be thrown away for efficient computations at each time step. As a result, the statistical information that is carried in high order terms, such as the higher order moments and the probabilities of the states, is generally not available from the method.

In this section, we develop an extension of the truncation method which enables us to carry high order terms in the computation. This is done by another truncation in [figure omitted; refer to PDF] at each step. Note that this is not just additional manipulation to simply make further elimination after the first truncation. By introducing the second truncation, we can alleviate the restriction on the original truncation in [figure omitted; refer to PDF] and allow higher order terms to remain, which is essential for computation of high order moments and probabilities.

In the second truncation, the terms with relatively small coefficients are dropped. It is justified from the fact that the coefficient of each term in the PGF is the probability of a specific state. We can therefore carry over high order terms, while keeping a reasonable number of terms only. This is important since evaluation of the high order moments and the probabilities is generally expensive, especially when the system involves stiff chemical reactions. Since DTM can preserve the high order terms without increasing the total number of the terms, it enables us to easily evaluate the probabilities for the states with many molecules.

The First Truncation. To overcome gradually increasing computational load in the series expansion (10), we restrict the series solution onto a short time interval. That is, on each interval [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , we construct the separate Taylor expansion [figure omitted; refer to PDF] at the center [figure omitted; refer to PDF] of finite order [figure omitted; refer to PDF] . We expect that side-by-side PGFs on consecutive intervals, [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , are equal on the boundary [figure omitted; refer to PDF] as [figure omitted; refer to PDF] More generally, [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is a "carry-over" functional. Our motivation of this successive superimpositions is to keep approximation order [figure omitted; refer to PDF] low in favor of computational efficiency while carrying over essential information like first several moments to a next function over the following time interval. Note that truncation of high order temporal terms can be justified in that, in PGF computation, [figure omitted; refer to PDF] represents how many reaction would happen for that interval [figure omitted; refer to PDF] , [figure omitted; refer to PDF] . Elimination of such terms has an effect of excluding excessive occurrence of reactions during such short interval.

One of the issues in the first truncation is how to determine carry-over functional [figure omitted; refer to PDF] . Although the best approximation in this framework would be obtained with the identity functional, it is not a practical choice since the number of terms involved in the computations keeps increasing rapidly. In [15], [figure omitted; refer to PDF] was suggested as a projection onto [figure omitted; refer to PDF] , a set of polynomials of degree three or less, which preserves the mean and the variance. However, such functional does not preserve the information related to high order derivatives of PGF and therefore approximation of quantities such as the state probabilities is not possible. Here we propose a different type of carry-over functional that preserves high order information. This functional indeed works as another truncation on what remains after the first truncation.

The Second Truncation. At time [figure omitted; refer to PDF] , the separate PGF [figure omitted; refer to PDF] becomes [figure omitted; refer to PDF] where [figure omitted; refer to PDF] means the probability for the state being [figure omitted; refer to PDF] at time [figure omitted; refer to PDF] . Now we define the carry-over functional [figure omitted; refer to PDF] such that it maps (22) into [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is a function defined as [figure omitted; refer to PDF] Application of such [figure omitted; refer to PDF] implies that we sieve the terms out of which coefficients are less than [figure omitted; refer to PDF] . One can see from (2) that those terms correspond to rare events of which coefficients probability is less than [figure omitted; refer to PDF] . Considering that the iterative process in the PGF approach is a linear operation with respect to [figure omitted; refer to PDF] , manipulation on a term only affects entailed terms, each of which represents possible consequence of the corresponding state. Since the probabilities of all the entailed events and their sum as well cannot exceed the probability of the eliminated event, the effect of elimination does not propagate or grow in probability.

Note that this truncation process is clearly consistent in that the scheme reverts to the first truncation as [figure omitted; refer to PDF] . However, one should be careful when applying the second truncation to approximation of the moments, since the moments do not solely depend on the state probabilities. For example, suppose that a rare state eventually gives a birth to large number of a species in an open system. If the mean is what is sought, one cannot simply omit such state for the reason that the corresponding probability is small, since the elimination may substantially change the resulted mean or variance in approximation. Hence, for the moment approximation of the open systems, it is better to use the above functional [figure omitted; refer to PDF] together with other carry-over functional as suggested in [15].

4.2. Example 1: Enzyme Kinetics

Once again, we use the enzyme kinetic model to test the performance of DTM. Figure 3 shows comparison of approximate solutions obtained from DTM and those from SSA. Note that, in the perturbation method in Section 3.2, evaluation of the probabilities was available only for the states that correspond to small number of that species, [figure omitted; refer to PDF] . This is because the probability computation (8) requires as many times of differentiation as the number of the species at the state. However, one can see that DTM successfully resolves this problem and produces the probabilities of the states corresponding as many as [figure omitted; refer to PDF] .

Figure 3: The initial condition [figure omitted; refer to PDF] and parameters [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] are assumed. In (c), each curve [figure omitted; refer to PDF] denotes the time-dependent probability solution that [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , respectively. The terms which have a coefficient less than [figure omitted; refer to PDF] are dropped out.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

4.3. Example 2: [figure omitted; refer to PDF] Transition Model

We consider the [figure omitted; refer to PDF] transition network in the eukaryotic cell cycle as follows [14, 17]: [figure omitted; refer to PDF]

We denote the number of molecules of [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] by [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , respectively.

We obtain a PDE of the moment generating function [figure omitted; refer to PDF] [14] [figure omitted; refer to PDF]

In Figure 4, we compare the simulation results obtained from the SSA and the DTM.

Figure 4: The initial condition [figure omitted; refer to PDF] and parameters [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] are assumed. In (c), each curve [figure omitted; refer to PDF] denotes the time-dependent probability solution that [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , respectively. The terms which have a coefficient less than [figure omitted; refer to PDF] are dropped out.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

5. Conlusion

In this work, we presented the double truncation method (DTM) that improves the probability generating function (PGF) method for computation of the solutions of stochastic reaction networks. Advantage of the PGF method is that, rather than struggling with the chemical master equation of extremely high order, one can handle a partial differential equation (PDE) of low order. However, it still requires intense computation for high order approximation, especially for the systems with stiff parameters. The method proposed in this paper applies the perturbation method to the corresponding PGF for reaction networks. That is, we expand the PGF in the parameter that is responsible for stiffness of the system. Then the PGF-PDE leads to a recurrence relation of the coefficient functions with respect to both the time and the corresponding parameter, which we can easily solve symbolically. It turns out that the suggested method requires a lower approximation order than the original PGF method, while it produces as good approximation as such.

While mean and variance of species are easily traceable in single truncation suggested in [15], efficient evaluation of the state probability was still a challenging task. In Section 4, we showed that the probability evaluation though PGF-PDEs can be achieved by doubly truncating the Taylor series at each time step. Through consecutive truncation processes, we dump relatively rare events, most of which are generated from excessive occurrence of reactions during a short time interval. This also implies elimination of a great deal of terms and significant reduction in computation.

It turns out in several numerical examples that the suggested method produces good approximation for the state probabilities. We expect that the DTM will be used for a fast and accurate approximation for computation of real and complex stochastic chemical reaction networks.

Acknowledgments

Pilwon Kim acknowledges that this work was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2011-0023486). Chang Hyeong Lee was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2014R1A1A2054976).

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] C. V. Rao, D. M. Wolf, A. P. Arkin, "Control, exploitation and tolerance of intracellular noise," Nature , vol. 420, no. 6912, pp. 231-237, 2002.

[2] D. T. Gillespie, "A rigorous derivation of the chemical master equation," Physica A: Statistical Mechanics and its Applications , vol. 188, no. 1-3, pp. 404-425, 1992.

[3] D. T. Gillespie, "Exact stochastic simulation of coupled chemical reactions," Journal of Physical Chemistry , vol. 81, no. 25, pp. 2340-2361, 1977.

[4] D. T. Gillespie, "Approximate accelerated stochastic simulation of chemically reacting systems," Journal of Chemical Physics , vol. 115, no. 4, pp. 1716-1733, 2001.

[5] Y. Cao, D. T. Gillespie, L. R. Petzold, "Efficient step size selection for the tauleaping simulation method," The Journal of Chemical Physics , vol. 124, no. 4, 2006.

[6] D. F. Anderson, "Incorporating postleap checks in tau-leaping," Journal of Chemical Physics , vol. 128, no. 5, 2008.

[7] Y. Cao, D. T. Gillespie, L. R. Petzold, "The slow-scale stochastic simulation algorithm," Journal of Chemical Physics , vol. 122, no. 1, 2005.

[8] B. Munsky, M. Khammash, "The finite state projection algorithm for the solution of the chemical master equation," Journal of Chemical Physics , vol. 124, no. 4, 2006.

[9] C. H. Lee, R. Lui, "A reduction method for multiple time scale stochastic reaction networks," Journal of Mathematical Chemistry , vol. 46, no. 4, pp. 1292-1321, 2009.

[10] C. H. Lee, R. Lui, "A reduction method for multiple time scale stochastic reaction networks with non-unique equilibrium probability," Journal of Mathematical Chemistry , vol. 47, no. 2, pp. 750-770, 2010.

[11] C. H. Lee, K.-H. Kim, P. Kim, "A moment closure method for stochastic reaction networks," Journal of Chemical Physics , vol. 130, no. 13, 2009.

[12] P. Milner, C. S. Gillespie, D. J. Wilkinson, "Moment closure approximations for stochastic kinetic models with rational rate laws," Mathematical Biosciences , vol. 231, no. 2, pp. 99-104, 2011.

[13] C. H. Lee, "A moment closure method for stochastic chemical reaction networks with general kinetics," MATCH. Communications in Mathematical and in Computer Chemistry , vol. 70, no. 3, pp. 785-800, 2013.

[14] P. Kim, C. H. Lee, "A probability generating function method for stochastic reaction networks," Journal of Chemical Physics , vol. 136, no. 23, 2012.

[15] P. Kim, C. H. Lee, "Fast probability generating function method for stochastic chemical reaction networks," MATCH: Communications in Mathematical and in Computer Chemistry , vol. 71, no. 1, pp. 57-69, 2014.

[16] D. A. McQuarrie, "Stochastic approach to chemical kinetics," Journal of Applied Probability , vol. 4, pp. 413-478, 1967.

[17] A. Kumar, K. Josic, "Reduced models of networks of coupled enzymatic reactions," Journal of Theoretical Biology , vol. 278, pp. 87-106, 2011.