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

Authors: Soyeong Jeong; Pilwon Kim; Chang Hyeong Lee
Journal: Journal of Applied Mathematics
 (2015)

---

## 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.


---

## Body

## 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(1)
∂
∂
t
p
n
,
t
=
∑
k
‍
a
k
n
-
V
k
·
p
n
-
V
k
,
t
-
∑
k
‍
a
k
n
·
p
n
,
t
,where n
(
t
)
=
(
n
1
(
t
)
,
n
2
(
t
)
,
…
,
n
s
(
t
)
), each n
i
(
t
) denotes the molecular number of ith species at time t, a
k is the propensity function or probability intensity for the kth reaction, and V
k denotes the kth column of the stoichiometric matrix V of which (
i
,
j
)th entry is the change in the number of molecules of the ith species by the occurrence of the jth 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 Section2, 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(2)
G
(
z
,
t
)
=
∑
n
=
0
∞
‍
z
n
p
(
n
,
t
)
,where z
=
(
z
1
,
…
,
z
s
), z
i
∈
[
-
1,1
], n
=
(
n
1
,
n
2
,
…
,
n
s
), and (3)
z
n
≡
z
1
n
1
z
2
n
2
⋯
z
s
n
s
.After differentiation of (2) with respect to t and application of (1), one can derive a PDE(4)
∂
G
∂
t
=
∑
|
k
|
≤
m
‍
H
k
(
z
,
t
)
D
k
G
.Here k
=
(
k
1
,
…
,
k
s
), |
k
|
=
k
1
+
⋯
+
k
s, m is the highest order of the reactions in the chemical system, and D
k
G denotes any kth order partial derivatives of G as (5)
D
k
G
=
∂
k
1
∂
z
1
k
1
∂
k
2
∂
z
2
k
2
⋯
∂
k
s
∂
z
s
k
s
G
.One can see that the PDE (4) is linear and H
k is a polynomial function of degree at most order m which satisfies H
k
(
z
=
1
,
t
)
=
0 for |
k
|
=
0,1
,
…
,
m.The conditions of the PDE (4) can be found from the conditions on n [14]. The initial condition is (6)
G
z
,
t
=
0
=
z
n
0
,where n
0 is the initial condition of n, and also it is true that(7)
G
(
z
=
0
,
t
)
=
p
(
n
=
0
,
t
)
,
G
(
z
=
1
,
t
)
=
∑
n
‍
p
(
n
,
t
)
=
1
.Using the solutionG 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 N
i
(
t
) denotes the number of molecules of ith species at time t, the probability that N
i
=
l at time t is(8)
p
(
N
i
=
l
,
t
)
=
1
l
!
∂
l
G
(
z
,
t
)
∂
z
i
l
z
i
=
0
,
z
j
=
1
,
j
≠
i
,and the mean and second moment are obtained by the first and second derivatives of G
(9)
∂
∂
z
i
G
z
=
1
,
t
=
E
N
i
t
,
∂
2
∂
z
i
∂
z
j
G
(
z
=
1
,
t
)
=
E
N
i
t
N
j
t
,
if
i
≠
j
,
E
N
i
2
t
-
E
N
i
t
,
if
i
=
j
,where E
[
X
(
t
)
] denotes the expectation of X
(
t
).


## 3. Perturbation of Probability Generating Function Method

### 3.1. Perturbation of PGF

For the PGF method recently proposed in [14], G
(
z
,
t
) is represented by a power expansion with respect to t as(10)
G
(
z
,
t
)
=
∑
n
=
0
∞
‍
f
n
z
t
n
,where f
n, n
=
0,1
,
…, are coefficient functions of z. Note that the initial condition determines the first coefficient as f
0
(
z
)
=
z
n
0.Suppose that the reaction rates in the network do not change with time and therefore the coefficients of the PDE do not involvet explicitly. This implies H
k
(
z
,
t
)
=
H
k
(
z
) in (4). Plugging (10) into the both sides of (4) and comparing the power series lead to a recursive relation between f
n
+
1
(
z
) and f
n
(
z
), n
=
0,1
,
…. One can derive coefficient functions f
1
(
z
)
,
f
2
(
z
)
,
…
,
f
J
(
z
) from f
1
(
z
)
=
z
n
0 and then construct an approximation ∑
n
=
0
J
‍
f
n
(
z
)
t
n for sufficiently large J. Since the initial function f
0
(
z
) and the coefficient of the PDE (2) are polynomials, all f
n
(
z
) 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 thatc
≪
1 is a small reaction rate in the system; then one can expand the PGF as(11)
G
z
,
t
=
g
0
z
,
t
+
c
g
1
z
,
t
+
c
2
g
2
z
,
t
+
⋯
,
=
∑
j
=
0
∞
‍
c
j
g
j
(
z
,
t
)
,where g
j
(
z
,
t
), j
=
0,1
,
2
,
…, are functions to be determined. As in (10), one can set for each g
j
(12)
g
j
(
z
,
t
)
=
∑
k
=
0
∞
‍
f
j
,
k
(
z
)
t
k
.Now, instead of (10), one can seek G
(
z
,
t
) in the following form:(13)
G
z
,
t
=
∑
j
=
0
∞
‍
c
j
∑
k
=
0
∞
‍
f
j
,
k
z
t
n
=
∑
j
=
0
M
‍
c
j
∑
k
=
0
N
‍
f
j
,
k
(
z
)
t
k
+
O
c
M
+
1
,
t
N
+
1for large M and N. By plugging (13) into the both sides of (4) and comparing the power series, one can find recurrence relations for f
j
,
k
,
f
j
+
1
,
k and f
j
+
1
,
k
+
1, where j
,
k
=
0,1
,
2
,
…. This decomposition enables us to handle stiffness caused by c properly, by adjusting the perturbation order M 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(14)
E
+
S
⇄
c
-
1
c
1
E
S
→
c
2
E
+
P
,where E, S, E
S, and P denote enzyme, substrate, enzyme-substrate complex, and product, respectively, and c
1, c
-
1, and c
2 are probability constants for reactions. If we denote the molecular numbers of E, S, E
S, and P by n
1, n
2, n
3, and n
4, respectively, the governing equation of the stochastic enzyme-substrate system is obtained [16] as(15)
∂
p
n
1
,
n
2
,
n
3
,
n
4
,
t
∂
t
=
c
1
n
1
+
1
n
2
+
1
·
p
(
n
1
+
1
,
n
2
+
1
,
n
3
-
1
,
n
4
,
t
)
+
c
-
1
(
n
3
+
1
)
p
(
n
1
-
1
,
n
2
-
1
,
n
3
+
1
,
n
4
,
t
)
+
c
2
(
n
3
+
1
)
p
(
n
1
-
1
,
n
2
,
n
3
+
1
,
n
4
-
1
,
t
)
-
(
c
1
n
1
n
2
+
c
-
1
n
3
+
c
2
n
3
)
p
(
n
1
,
n
2
,
n
3
,
n
4
,
t
)
.One can see that there are two conservation quantities in the system and thus we can remove two of variables, sayn
3, n
4. After the removal, we can derive a PDE of G
(
z
1
,
z
2
,
t
) from the master equation [14] (16)
G
t
=
c
1
1
-
z
1
z
2
G
12
+
c
-
1
+
c
2
z
1
-
c
-
1
z
1
2
z
2
-
c
2
z
1
2
G
1
+
A
c
-
1
z
1
z
2
+
c
2
z
1
-
c
-
1
+
c
2
G
.Here A is the conserved quantity n
1
+
n
3 and we use the notations (17)
G
i
=
∂
G
∂
z
i
,
G
i
j
=
∂
2
G
∂
z
i
∂
z
j
,
i
,
j
=
1,2
,
…
.In the example, we use parametersc
1
=
0.01
s
-
1, c
-
1
=
1
s
-
1, and c
3
=
1
s
-
1 and let the perturbation be expanded in c
1. 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 thatn
1
(
0
)
=
5, n
2
(
0
)
=
5, and A
=
5. (a) shows the comparison of exact values and approximate values obtained from the PPGF method for mean and mean ± standard deviations (±
σ) of n
1 and (b) shows the comparison of the exact probability and the approximate probability obtained from the PPGF method. In the figure, P
(
s
=
i
) denotes the probability that the number n
2 of the substrate is i. The exact probability, mean, and variance are obtained by rewriting (15) into the linear ODE system d
p
/
d
t
=
K
p under the given initial condition and then solving it by MATLAB.
(a)
(b)


### 3.3. Example 2: Brusselator Model

We consider the Brusselator model:(18)
A
→
c
1
X
,
2
X
+
Y
→
c
2
3
X
,
B
+
X
→
c
3
Y
+
D
,
X
→
c
4
E
.Here X, Y are molecular species and A, B, D, and E 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 X, Y by n
1, n
2, respectively, we obtain the PDE for the probability generating function as follows:(19)
G
t
=
c
1
A
z
1
-
1
G
+
c
2
2
z
1
2
z
1
-
z
2
G
112
+
c
3
B
(
z
2
-
z
1
)
G
1
+
c
4
(
1
-
z
1
)
G
1
.Here A
=
B
=
1 and the reaction rates c
1
=
1
s
-
1, c
2
=
0.0001
s
-
1, and c
3
=
c
4
=
0.1 are used. We use PPGF with respect to the smallest parameter c
2. 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 conditionn
1
(
0
)
=
0 and n
2
(
0
)
=
0 are assumed. (a) shows a typical run of Y in 500 seconds from a Monte Carlo simulation of the Brusselator model. In (b), P
(
y
=
i
) denotes the probability that the number n
2 of the substrate is i. (c) and (d) compare mean and mean ± standard deviations (±
σ) of n
1, n
2, and n
3 obtained from the SSA and the PPGF. The results from the SSA are based on 50,000 realizations.
(a)
(b)
(c)
(d)


## 3.1. Perturbation of PGF

For the PGF method recently proposed in [14], G
(
z
,
t
) is represented by a power expansion with respect to t as(10)
G
(
z
,
t
)
=
∑
n
=
0
∞
‍
f
n
z
t
n
,where f
n, n
=
0,1
,
…, are coefficient functions of z. Note that the initial condition determines the first coefficient as f
0
(
z
)
=
z
n
0.Suppose that the reaction rates in the network do not change with time and therefore the coefficients of the PDE do not involvet explicitly. This implies H
k
(
z
,
t
)
=
H
k
(
z
) in (4). Plugging (10) into the both sides of (4) and comparing the power series lead to a recursive relation between f
n
+
1
(
z
) and f
n
(
z
), n
=
0,1
,
…. One can derive coefficient functions f
1
(
z
)
,
f
2
(
z
)
,
…
,
f
J
(
z
) from f
1
(
z
)
=
z
n
0 and then construct an approximation ∑
n
=
0
J
‍
f
n
(
z
)
t
n for sufficiently large J. Since the initial function f
0
(
z
) and the coefficient of the PDE (2) are polynomials, all f
n
(
z
) 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 thatc
≪
1 is a small reaction rate in the system; then one can expand the PGF as(11)
G
z
,
t
=
g
0
z
,
t
+
c
g
1
z
,
t
+
c
2
g
2
z
,
t
+
⋯
,
=
∑
j
=
0
∞
‍
c
j
g
j
(
z
,
t
)
,where g
j
(
z
,
t
), j
=
0,1
,
2
,
…, are functions to be determined. As in (10), one can set for each g
j
(12)
g
j
(
z
,
t
)
=
∑
k
=
0
∞
‍
f
j
,
k
(
z
)
t
k
.Now, instead of (10), one can seek G
(
z
,
t
) in the following form:(13)
G
z
,
t
=
∑
j
=
0
∞
‍
c
j
∑
k
=
0
∞
‍
f
j
,
k
z
t
n
=
∑
j
=
0
M
‍
c
j
∑
k
=
0
N
‍
f
j
,
k
(
z
)
t
k
+
O
c
M
+
1
,
t
N
+
1for large M and N. By plugging (13) into the both sides of (4) and comparing the power series, one can find recurrence relations for f
j
,
k
,
f
j
+
1
,
k and f
j
+
1
,
k
+
1, where j
,
k
=
0,1
,
2
,
…. This decomposition enables us to handle stiffness caused by c properly, by adjusting the perturbation order M 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(14)
E
+
S
⇄
c
-
1
c
1
E
S
→
c
2
E
+
P
,where E, S, E
S, and P denote enzyme, substrate, enzyme-substrate complex, and product, respectively, and c
1, c
-
1, and c
2 are probability constants for reactions. If we denote the molecular numbers of E, S, E
S, and P by n
1, n
2, n
3, and n
4, respectively, the governing equation of the stochastic enzyme-substrate system is obtained [16] as(15)
∂
p
n
1
,
n
2
,
n
3
,
n
4
,
t
∂
t
=
c
1
n
1
+
1
n
2
+
1
·
p
(
n
1
+
1
,
n
2
+
1
,
n
3
-
1
,
n
4
,
t
)
+
c
-
1
(
n
3
+
1
)
p
(
n
1
-
1
,
n
2
-
1
,
n
3
+
1
,
n
4
,
t
)
+
c
2
(
n
3
+
1
)
p
(
n
1
-
1
,
n
2
,
n
3
+
1
,
n
4
-
1
,
t
)
-
(
c
1
n
1
n
2
+
c
-
1
n
3
+
c
2
n
3
)
p
(
n
1
,
n
2
,
n
3
,
n
4
,
t
)
.One can see that there are two conservation quantities in the system and thus we can remove two of variables, sayn
3, n
4. After the removal, we can derive a PDE of G
(
z
1
,
z
2
,
t
) from the master equation [14] (16)
G
t
=
c
1
1
-
z
1
z
2
G
12
+
c
-
1
+
c
2
z
1
-
c
-
1
z
1
2
z
2
-
c
2
z
1
2
G
1
+
A
c
-
1
z
1
z
2
+
c
2
z
1
-
c
-
1
+
c
2
G
.Here A is the conserved quantity n
1
+
n
3 and we use the notations (17)
G
i
=
∂
G
∂
z
i
,
G
i
j
=
∂
2
G
∂
z
i
∂
z
j
,
i
,
j
=
1,2
,
…
.In the example, we use parametersc
1
=
0.01
s
-
1, c
-
1
=
1
s
-
1, and c
3
=
1
s
-
1 and let the perturbation be expanded in c
1. 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 thatn
1
(
0
)
=
5, n
2
(
0
)
=
5, and A
=
5. (a) shows the comparison of exact values and approximate values obtained from the PPGF method for mean and mean ± standard deviations (±
σ) of n
1 and (b) shows the comparison of the exact probability and the approximate probability obtained from the PPGF method. In the figure, P
(
s
=
i
) denotes the probability that the number n
2 of the substrate is i. The exact probability, mean, and variance are obtained by rewriting (15) into the linear ODE system d
p
/
d
t
=
K
p under the given initial condition and then solving it by MATLAB.
(a)
(b)


## 3.3. Example 2: Brusselator Model

We consider the Brusselator model:(18)
A
→
c
1
X
,
2
X
+
Y
→
c
2
3
X
,
B
+
X
→
c
3
Y
+
D
,
X
→
c
4
E
.Here X, Y are molecular species and A, B, D, and E 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 X, Y by n
1, n
2, respectively, we obtain the PDE for the probability generating function as follows:(19)
G
t
=
c
1
A
z
1
-
1
G
+
c
2
2
z
1
2
z
1
-
z
2
G
112
+
c
3
B
(
z
2
-
z
1
)
G
1
+
c
4
(
1
-
z
1
)
G
1
.Here A
=
B
=
1 and the reaction rates c
1
=
1
s
-
1, c
2
=
0.0001
s
-
1, and c
3
=
c
4
=
0.1 are used. We use PPGF with respect to the smallest parameter c
2. 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 conditionn
1
(
0
)
=
0 and n
2
(
0
)
=
0 are assumed. (a) shows a typical run of Y in 500 seconds from a Monte Carlo simulation of the Brusselator model. In (b), P
(
y
=
i
) denotes the probability that the number n
2 of the substrate is i. (c) and (d) compare mean and mean ± standard deviations (±
σ) of n
1, n
2, and n
3 obtained from the SSA and the PPGF. The results from the SSA are based on 50,000 realizations.
(a)
(b)
(c)
(d)


## 4. Double Truncation Method

### 4.1. Truncations int and z

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 inz at each step. Note that this isnot 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 t 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 [
t
r
,
t
r
+
1
], r
=
0,1
,
2
,
…, t
0
=
0
<
t
1
<
t
2
<
⋯, we construct the separate Taylor expansion G
[
r
]
(
z
,
t
) at the center t
=
t
r of finite order N. We expect that side-by-side PGFs on consecutive intervals, G
[
r
]
(
z
,
t
) and G
[
r
+
1
]
(
z
,
t
), are equal on the boundary t
=
t
r
+
1 as (20)
G
[
r
]
(
z
,
t
r
+
1
)
≈
G
[
r
+
1
]
(
z
,
t
r
+
1
)
.More generally, (21)
G
0
z
,
t
=
∑
n
=
0
N
‍
f
n
0
z
t
n
,
G
0
z
,
0
=
z
n
,
G
[
1
]
(
z
,
t
)
=
∑
n
=
0
N
‍
f
n
[
1
]
(
z
)
t
-
t
1
n
,
G
1
z
,
t
1
=
F
G
0
z
,
t
1
,
G
2
(
z
,
t
)
=
∑
n
=
0
N
‍
f
n
2
(
z
)
t
-
t
2
n
,
G
2
z
,
t
2
=
F
G
1
z
,
t
2
,
aaaaaaaaaa
⋮
G
[
r
]
(
z
,
t
)
=
∑
n
=
0
N
‍
f
n
[
r
]
(
z
)
t
-
t
r
n
,
G
r
z
,
t
r
=
F
G
r
-
1
z
,
t
r
,
aaaaaaaaaaa
⋮where F is a “carry-over” functional. Our motivation of this successive superimpositions is to keep approximation order N 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, N represents how many reaction would happen for that interval [
t
r
,
t
r
+
1
], r
=
0,1
,
2
,
…. 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 functionalF. 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], F was suggested as a projection onto P
3
(
z
), 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 t
=
t
r
+
1, the separate PGF G
[
r
] becomes(22)
G
[
r
]
(
z
,
t
r
+
1
)
=
∑
n
=
0
∞
‍
z
n
p
[
r
]
(
n
,
t
r
+
1
)
,where p
[
r
]
(
n
,
t
r
+
1
) means the probability for the state being n at time t
r
+
1. Now we define the carry-over functional F such that it maps (22) into (23)
∑
n
=
0
∞
‍
z
n
I
ϵ
p
r
n
,
t
r
+
1
,where H
ϵ is a function defined as (24)
I
ϵ
(
x
)
=
x
,
if
x
>
ϵ
,
0
,
if
x
≤
ϵ
.Application of such F implies that we sieve the terms out of which coefficients are less than ϵ. One can see from (2) that those terms correspond to rare events of which coefficients probability is less than ϵ. Considering that the iterative process in the PGF approach is a linear operation with respect to z
n, 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ϵ
→
0. 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 F 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. Figure3 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, n
2
=
0,1
,
2. 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 n
2
=
10,15.Figure 3
The initial conditionn
=
[
10,20,20,0
] and parameters c
1
=
0.1, c
-
1
=
1, and c
2
=
0.5 are assumed. In (c), each curve p
(
n
2
=
i
) denotes the time-dependent probability solution that n
2
=
i, i
=
5,10,17, respectively. The terms which have a coefficient less than 1
0
-
7 are dropped out.
(a)
(b)
(c)


### 4.3. Example 2:G
2
/
M Transition Model

We consider theG
2
/
M transition network in the eukaryotic cell cycle as follows [14, 17]: (25)
X
+
Y
p
⇄
c
2
c
1
C
x
→
c
3
X
p
+
Y
p
,
E
1
+
X
p
⇄
c
5
c
4
C
x
e
→
c
6
X
+
E
1
X
p
+
Y
⇄
c
8
c
7
C
y
→
c
9
X
p
+
Y
p
,
E
2
+
Y
p
⇄
c
11
c
10
C
y
e
→
c
12
Y
+
E
2
.We denote the number of molecules ofX
p, Y
p, X, Y, E
1, E
2, C
x, C
x
e, C
y, and C
y
e by n
1, n
2, n
3, n
4, n
5, n
6, n
7, n
8, n
9, and n
10, respectively.We obtain a PDE of the moment generating functionG
(
z
,
t
) [14](26)
G
t
=
c
1
G
23
z
7
-
z
2
z
3
+
c
2
G
7
z
2
z
3
-
z
7
+
c
3
G
7
z
1
z
2
-
z
7
+
c
4
G
15
z
8
-
z
1
z
5
+
c
5
G
8
(
z
1
z
5
-
z
8
)
+
c
6
G
8
(
z
3
z
5
-
z
8
)
+
c
7
G
14
(
z
9
-
z
1
z
4
)
+
c
8
G
9
(
z
1
z
4
-
z
9
)
+
c
9
G
9
z
1
z
2
-
z
9
+
c
10
G
26
z
10
-
z
2
z
6
+
c
11
G
10
(
z
2
z
6
-
z
10
)
+
c
12
G
10
(
z
4
z
6
-
z
10
)
.In Figure4, we compare the simulation results obtained from the SSA and the DTM.Figure 4
The initial conditionn
=
[
5,4
,
2,3
,
2,3
,
2,0
,
0,0
] and parameters c
1
=
c
4
=
c
7
=
c
10
=
0.2, c
2
=
c
5
=
c
8
=
c
11
=
1, and c
3
=
c
6
=
c
9
=
c
12
=
0.1 are assumed. In (c), each curve p
(
n
2
=
i
) denotes the time-dependent probability solution that n
2
=
i, i
=
1,3
,
5, respectively. The terms which have a coefficient less than 1
0
-
8 are dropped out.
(a)
(b)
(c)


## 4.1. Truncations int and z

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 inz at each step. Note that this isnot 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 t 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 [
t
r
,
t
r
+
1
], r
=
0,1
,
2
,
…, t
0
=
0
<
t
1
<
t
2
<
⋯, we construct the separate Taylor expansion G
[
r
]
(
z
,
t
) at the center t
=
t
r of finite order N. We expect that side-by-side PGFs on consecutive intervals, G
[
r
]
(
z
,
t
) and G
[
r
+
1
]
(
z
,
t
), are equal on the boundary t
=
t
r
+
1 as (20)
G
[
r
]
(
z
,
t
r
+
1
)
≈
G
[
r
+
1
]
(
z
,
t
r
+
1
)
.More generally, (21)
G
0
z
,
t
=
∑
n
=
0
N
‍
f
n
0
z
t
n
,
G
0
z
,
0
=
z
n
,
G
[
1
]
(
z
,
t
)
=
∑
n
=
0
N
‍
f
n
[
1
]
(
z
)
t
-
t
1
n
,
G
1
z
,
t
1
=
F
G
0
z
,
t
1
,
G
2
(
z
,
t
)
=
∑
n
=
0
N
‍
f
n
2
(
z
)
t
-
t
2
n
,
G
2
z
,
t
2
=
F
G
1
z
,
t
2
,
aaaaaaaaaa
⋮
G
[
r
]
(
z
,
t
)
=
∑
n
=
0
N
‍
f
n
[
r
]
(
z
)
t
-
t
r
n
,
G
r
z
,
t
r
=
F
G
r
-
1
z
,
t
r
,
aaaaaaaaaaa
⋮where F is a “carry-over” functional. Our motivation of this successive superimpositions is to keep approximation order N 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, N represents how many reaction would happen for that interval [
t
r
,
t
r
+
1
], r
=
0,1
,
2
,
…. 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 functionalF. 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], F was suggested as a projection onto P
3
(
z
), 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 t
=
t
r
+
1, the separate PGF G
[
r
] becomes(22)
G
[
r
]
(
z
,
t
r
+
1
)
=
∑
n
=
0
∞
‍
z
n
p
[
r
]
(
n
,
t
r
+
1
)
,where p
[
r
]
(
n
,
t
r
+
1
) means the probability for the state being n at time t
r
+
1. Now we define the carry-over functional F such that it maps (22) into (23)
∑
n
=
0
∞
‍
z
n
I
ϵ
p
r
n
,
t
r
+
1
,where H
ϵ is a function defined as (24)
I
ϵ
(
x
)
=
x
,
if
x
>
ϵ
,
0
,
if
x
≤
ϵ
.Application of such F implies that we sieve the terms out of which coefficients are less than ϵ. One can see from (2) that those terms correspond to rare events of which coefficients probability is less than ϵ. Considering that the iterative process in the PGF approach is a linear operation with respect to z
n, 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ϵ
→
0. 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 F 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. Figure3 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, n
2
=
0,1
,
2. 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 n
2
=
10,15.Figure 3
The initial conditionn
=
[
10,20,20,0
] and parameters c
1
=
0.1, c
-
1
=
1, and c
2
=
0.5 are assumed. In (c), each curve p
(
n
2
=
i
) denotes the time-dependent probability solution that n
2
=
i, i
=
5,10,17, respectively. The terms which have a coefficient less than 1
0
-
7 are dropped out.
(a)
(b)
(c)


## 4.3. Example 2:G
2
/
M Transition Model

We consider theG
2
/
M transition network in the eukaryotic cell cycle as follows [14, 17]: (25)
X
+
Y
p
⇄
c
2
c
1
C
x
→
c
3
X
p
+
Y
p
,
E
1
+
X
p
⇄
c
5
c
4
C
x
e
→
c
6
X
+
E
1
X
p
+
Y
⇄
c
8
c
7
C
y
→
c
9
X
p
+
Y
p
,
E
2
+
Y
p
⇄
c
11
c
10
C
y
e
→
c
12
Y
+
E
2
.We denote the number of molecules ofX
p, Y
p, X, Y, E
1, E
2, C
x, C
x
e, C
y, and C
y
e by n
1, n
2, n
3, n
4, n
5, n
6, n
7, n
8, n
9, and n
10, respectively.We obtain a PDE of the moment generating functionG
(
z
,
t
) [14](26)
G
t
=
c
1
G
23
z
7
-
z
2
z
3
+
c
2
G
7
z
2
z
3
-
z
7
+
c
3
G
7
z
1
z
2
-
z
7
+
c
4
G
15
z
8
-
z
1
z
5
+
c
5
G
8
(
z
1
z
5
-
z
8
)
+
c
6
G
8
(
z
3
z
5
-
z
8
)
+
c
7
G
14
(
z
9
-
z
1
z
4
)
+
c
8
G
9
(
z
1
z
4
-
z
9
)
+
c
9
G
9
z
1
z
2
-
z
9
+
c
10
G
26
z
10
-
z
2
z
6
+
c
11
G
10
(
z
2
z
6
-
z
10
)
+
c
12
G
10
(
z
4
z
6
-
z
10
)
.In Figure4, we compare the simulation results obtained from the SSA and the DTM.Figure 4
The initial conditionn
=
[
5,4
,
2,3
,
2,3
,
2,0
,
0,0
] and parameters c
1
=
c
4
=
c
7
=
c
10
=
0.2, c
2
=
c
5
=
c
8
=
c
11
=
1, and c
3
=
c
6
=
c
9
=
c
12
=
0.1 are assumed. In (c), each curve p
(
n
2
=
i
) denotes the time-dependent probability solution that n
2
=
i, i
=
1,3
,
5, respectively. The terms which have a coefficient less than 1
0
-
8 are dropped out.
(a)
(b)
(c)


## 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.



---

*Source: 864238-2015-05-28.xml*