Double Truncation Method for Simulation of Stochastic Chemical Reaction Networks

Abstract

We develop Double Truncation Method(DTM) to understand the stochasticity of chemical reaction network on the mesoscopic scale. The Chemical Master Equation(CME) describes the probability distribution of the system accurately under the markov assumption. But solving CME is computationally heavy. It faces the curse of dimensionality because it considers reactions on every different states. To settle the issue, stochastic approach came out, typically Gillspie's stochastic simulation algorithm(SSA). SSA lifts the curse of dimensionality, but it needs too many realizations, which makes it less practical, in cases of the system consisting of large numbers of molecules or very different time scale reactions. A recently developed Probability Generating Function(PGF) method supplements those weaknesses. It is a deterministic description and sparks the reactions in stead of considering those for all the states. By doing that, though it expresses the system efficiently, but it implements the symbolic computation and converges still slowly. So here we suggest DTM to speed up. As suggested from the name, DTM has two truncations for time and for coefficients based on PGF method. We perform the first truncation for a short time and second truncation for small coefficients at each time step. First truncation or superimposition can be performed underlying the power series expansion on time. And next truncation can be conducted with the elimination of relatively small coefficient terms. Since the coefficient of PGF means the probability of a specific state, the sum of coefficient can be understood as an weight for the system. This observation enables us to ignore a great deal of small terms which do not affect the system significantly. The method is procedurally simple and powerful, especially for mesoscopic scale problems. It works well even for open systems, such as brusselator. We apply the method to simulation of binding reactions, enzyme kinetics, transition model and brusselator and compare the results with those of SSA or matrix exponential.