Journal of Applied Mathematics

Volume 2015, Article ID 864238, 8 pages

http://dx.doi.org/10.1155/2015/864238

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

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

Academic Editor: Lotfollah Najjar

Copyright © 2015 Soyeong Jeong et al. 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.

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

#### 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 aswhere , each denotes the molecular number of th species at time , is the propensity function or probability intensity for the th reaction, and denotes the th column of the stoichiometric matrix of which th entry is the change in the number of molecules of the th species by the occurrence of the 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 aswhere , , , and After differentiation of (2) with respect to and application of (1), one can derive a PDEHere , , is the highest order of the reactions in the chemical system, and denotes any th order partial derivatives of as One can see that the PDE (4) is linear and is a polynomial function of degree at most order which satisfies for .

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

Using the solution 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 denotes the number of molecules of th species at time , the probability that at time isand the mean and second moment are obtained by the first and second derivatives of where denotes the expectation of .

#### 3. Perturbation of Probability Generating Function Method

##### 3.1. Perturbation of PGF

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

Suppose that the reaction rates in the network do not change with time and therefore the coefficients of the PDE do not involve explicitly. This implies in (4). Plugging (10) into the both sides of (4) and comparing the power series lead to a recursive relation between and , . One can derive coefficient functions from and then construct an approximation for sufficiently large . Since the initial function and the coefficient of the PDE (2) are polynomials, all 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 is a small reaction rate in the system; then one can expand the PGF aswhere , , are functions to be determined. As in (10), one can set for each Now, instead of (10), one can seek in the following form:for large and . By plugging (13) into the both sides of (4) and comparing the power series, one can find recurrence relations for and , where . This decomposition enables us to handle stiffness caused by properly, by adjusting the perturbation order 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 systemwhere , , , and denote enzyme, substrate, enzyme-substrate complex, and product, respectively, and , , and are probability constants for reactions. If we denote the molecular numbers of , , , and by , , , and , respectively, the governing equation of the stochastic enzyme-substrate system is obtained [16] as

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

In the example, we use parameters , , and and let the perturbation be expanded in . Figure 1 shows comparison of the exact solution of the master equation (15) and the approximate solution obtained from the PPGF method.