Approximate Exponential Algorithms to Solve the Chemical Master Equation

This paper discusses new simulation algorithms for stochastic chemical kinetics that exploit the linearity of the chemical master equation and its matrix exponential exact solution. These algorithms make use of various approximations of the matrix exponential to evolve probability densities in time. A sampling of the approximate solutions of the chemical master equation is used to derive accelerated stochastic simulation algorithms. Numerical experiments compare the new methods with the established stochastic simulation algorithm and the tau-leaping method.


Introduction
In many biological systems the small number of participating molecules make the chemical reactions inherently stochastic. The system state is described by probability densities of the numbers of molecules of different species. The evolution of probabilities in time is described by the chemical master equation (CME) [2]. Gillespie proposed the Stochastic Simulation Algorithm (SSA), a Monte Carlo approach that samples from CME [2]. SSA became the standard method for solving well-stirred chemically reacting systems. However, SSA simulates one reaction and is inefficient for most realistic problems. This motivated the quest for approximate sampling techniques to enhance the efficiency.
The first approximate acceleration technique is the tau-leaping method [3] which is able to simulate multiple chemical reactions appearing in a pre-selected time step of length τ . The tau-leap method is accurate if τ is small enough to satisfy the leap condition, meaning that propensity functions remain nearly constant in a time step. The number of firing reactions in a time step is approximated by a Poisson random variable [5]. Explicit tau-leaping method is numerically unstable for stiff systems [15]. Stiffness systems have well-separated "fast" and "slow" time scales present, and the "fast modes" are stable. The implicit tau-leap method [6] overcomes the stability issue but it has a damping effect on the computed variances. More accurate variations of the implicit tauleap method have been proposed to alleviate the damping [3,4,7,11,13,14]. Simulation efficiency has been increased via parallelization [12].
Direct solutions of the CME are computationally important specially in order to estimate moments of the distributions of the chemical species [8]. Various approaches to solve the CME are discussed in [1].
Sandu has explained the explicit tau-leap method as an exact sampling procedure from an approximate solution of the CME [9]. This paper extends that study and proposes new approximations to the CME solution based on various approximations of matrix exponentials. Accelerated stochastic simulation algorithms are the built by performing exact sampling of these approximate probability densities.

2
The paper is organized as follows. Section 2 reviews the stochastic simulation of chemical kinetics. Section 3 developed the new approximation methods. Numerical experiments to illustrate the proposed schemes are carried out in Section 4. Conclusions are drawn in Section 5.
2 Simulation of stochastic chemical kinetics Consider a chemical system in a constant volume container. The system is well-stirred and in thermal equilibrium at some constant temperature. There are N different chemical species S 1 , . . . , S N . Let X i (t) denote the number of molecules of species S i at time t. The state vector x(t) = [X 1 (t), . . . , X N (t)] defines the numbers of molecules of each species present at time t. The chemical network consists of M reaction channels R 1 , . . . , R M . Each individual reaction destroys a number of molecules of reactant species, and produces a number of molecules of the products. Let ν i j be the change in the number of S i molecules caused by a single reaction R j . The state change vector ν j = [ν 1 j , . . . , ν N j ] describes the change in the entire state following R j .
A propensity function a j (x) is associated with each reaction channel R j . The probability that one R j reaction will occur in the next infinitesimal time interval [t, t+dt) is a j (x(t))·dt. The purpose of a stochastic chemical simulation is to trace the time evolution of the system state x(t) given that at the initial timet the system is in the initial state x (t).

Chemical Master Equation
The Chemical Master Equation (CME) [2] has complete information about time evolution of probability of system's state Let Q i be the total possible number of molecules of species S i . The total number of all possible states of the system is: We denote by I(x) the state-space index of state x = [X 1 , . . . , X N ] One firing of reaction R r changes the state from x tox = x − v r . The corresponding change in state space index is:

Approximation to Chemical Master Equation
The discrete solutions of the CME (1) are vectors in the discrete state space, P (t) ∈ R Q . Consider the diagonal matrix A 0 ∈ R Q×Q and the Toeplitz matrices as well as their sum A ∈ R Q×Q with entries where x j denotes the unique state with state space index j = I(x j ). In fact matrix A is a square (Q × Q) matrix which contains all the propensity values for each possible value of all species or let's say all possible states of reaction system. All possible states for a reaction system consists of N species where each specie has at most Q i i = 1, 2, ..., N value.
The CME (1) is a linear ODE on the discrete state space where the system is initially in the known state x(0) =x and therefore the initial probability distribution vector P(0) ∈ R Q is equal to one at I(x) and is zero everywhere else. The exact solution of the linear ODE (3) is follows:

Approximation to Chemical Master Equation
Although the CME (1) fully describes the evolution of probabilities it is difficult to solve in practice due to large state space. Sandu [9] considers the following approximation of the CME: where the arguments of all propensity functions have been changed from x or x − v j tox. In order to obtain an exponential solution to (5) in probability space we consider the diagonal matrixĀ 0 ∈ R Q×Q and the Toeplitz matrices A 1 , ...,Ā M ∈ R Q×Q [9].Ā r matrices are square (Q × Q) matrices are built upon the current state of system in reaction system which is against A r matrices that contain all possible states of reaction system.

Tau-leaping method
together with their sumĀ =Ā 0 + · · · +Ā M . The approximate CME (5) can be written as the linear ODE and has an exact solution

Tau-leaping method
In tau-leap method the number of times a reaction fires is a random variable from a Poisson distribution with parameter a r (x) τ . Since each reaction fires independently, the probability that each reaction R r fires exactly k r times, r = 1, 2, · · · , M , is the product of M Poisson probabilities.
Then the state vector after these reactions will change as follows: The probability to go from statex att to state x att + τ , P (X (t + τ )) = x, is the sum of all possible firing reactions which is: (a r (xT )) kr K r ! Equation (7) can be approximated by product of each matrix exponential: It has been shown in [9] that the probability given by the tau-leaping method is exactly the probability evolved by the approximate solution (9).
3 Approximations to the exponential solution

Column based splitting
In column based splitting the matrix A (2) is decomposed in a sum of columns Each matrix A j has the same j-th column as the matrix A, and is zero everywhere else. Here c j is the j th column of matrix A and e j is the canonical vector which is zero every where except the j th component. The exponential of τ A j is: Since e T j c j is equal to the j-th diagonal entry of matrix A: Consequently the matrix exponential (11) becomes We have and the approximation to the CME solution reads

Accelerated tau-leaping
In this approximation method we build the matrices where a r (x) are the propensity functions. The matrix A in (2) can be written as

Symmetric accelerated tau-leaping
The solution of the linear CME (4) can be approximated by Note that the evolution of state probability by e τ Bj ·P (t) describes the change in probability when only reaction j fires in the time interval τ . The corresponding evolution of the number of molecules that samples the evolved probability is where K (a j (x (t)) τ ) is a random number drawn from a Poisson distribution with parameter a j (x (t)) τ , and V j is the j-th column of stoichiometry matrix. The approximate solution (12) accounts for the change in probability due to a sequential firing of reactions M , M − 1, down to 1. Sampling from the resulting probability density can be done by changing the system state sequentially consistent with the firing of each reaction. This results in the following accelerated tau-leaping algorithm: Moreover, (12) can also be written as: Then, (13) can be written as:

Symmetric accelerated tau-leaping
A more accurate version of accelerated tau-leaping can be constructed by using symmetric Strang splitting (10) to approximate the matrix exponential in (12). Following the procedure used to derive (13) leads to the following sampling algorithm:X (16)

Numerical experiments
The above approximation techniques are used to solve two test systems, reversible isomer and the Schlogl reactions [15]. The experimental results are presented in following sections.

Isomer reaction
The reversible isomer reaction system is [15] x 1 The stoichiometry matrix and the propensity functions are: The exact exponential solution of CME obtained from (4) is a joint probability distribution vector for the two species at final time. Figure 1(a) shows that the histogram of 10,000 SSA solutions is very close to the exact exponential solution. The approximate solution using the sum of exponentials (7) is illustrated in Figure 1(b). This approximation is not very accurate since it uses only the current state of the system. Other approximation methods based on the product of exponentials (9) and Strang splitting (10) are not very strong approximations as the exact solution hence, the results are not reported.
The results reported in Figure 2 indicate that for small time steps τ the accelerated tau-leap (13) solution is very close to the results provided by traditional explicit tau-leap. Symmetric accelerated tau-leap method (16) yields even better results, as shown in Figure 3. For small time steps the traditional and symmetric accelerated methods give similar results, however, for large time steps, the results of the symmetric accelerated method is considerably more stable.

Schlogl reaction
We next consider the Schlogl reaction system [15] whose solution has a bi-stable distribution. Let N 1 , N 2 be the numbers of molecules of species B 1 and B 2 , respectively. The reaction stoichiometry matrix and the propensity functions are:   accuracy can be improved to some extent using the strategies described in (14) and (15).

Conclusions
This study proposes new numerical solvers for stochastic simulations of chemical kinetics. The proposed approach exploits the linearity of the CME and the exponential form of its exact solution. The matrix exponential appearing in the CME solution is approximated as a product of simpler matrix exponentials. This leads to an approximate ("numerical") solution of the probability density evolved to a future time. The solution algorithms sample exactly this approximate probability density and provide extensions of the traditional tauleap approach.
Different approximations of the matrix exponential lead to different numerical algorithms: Strang splitting, column splitting, accelerated tau-leap, and symmetric accelerated tau-leap. Current work by the authors focuses on improving the accuracy of these novel approximation techniques for stochastic chemical kinetics.
The matrix A ∈ R Q·Q×Q·Q As an example for a maximum number of species Q 1 = 2, Q 2 = 2 the matrix A is:
As an example matrix A for maximum number of molecules Q = 5 is the following tridiagonal matrix: