Moment Propagation Methods for Stochastic Simulation of Complex Biochemical Systems

Abstract

We are interested in predicting the time dependent behavior of biochemical networks such as interaction between proteins. These networks are represented by a system of chemical reactions. In the forward problem, we have knowledge of all the reactions and the associated rate constants. We want to determine the joint probability density of the populations of all the molecular species at any time instant. The given biochemical system is a discrete state, continuous time Markov process, and the time evolution of its probability density function is described by the chemical master equation (CME). We want to obtain the solution of the master equation in complex biochemical systems with a large number of species and reactions. Analytical solutions of the chemical master equation for first and second order reactions have only been obtained for select cases with a few species and reactions. Another approach to solving the problem is to approximate the CME with the Fokker-Planck equations. But this would require solving partial differential equations with a large number of variables. The present state of the art approaches for stochastic simulation of such systems are based on Monte Carlo methods. One such popular method is the Stochastic Simulation Algorithm (SSA) derived by Gillespie in 1976. Several authors have developed accelerated versions of the SSA such as the Next Reaction Method and time leaping methods, in order to reduce the computation time of SSA. The Monte Carlo methods provide approximations of the complete distribution, but they require simulations of many realizations of the Markov process and many time steps. Hence their computation times are prohibitively long for very complex systems. Methods for modeling the biochemical networks based on moment propagation is a relatively unexplored area. We propose a new method for propagating the first two moments of the joint probability distribution of the number of molecules. In many systems, the distribution can be approximated as Gaussians and therefore computing the first two moments is sufficient. Simulation results show that our method yields accurate results for first order and second order reactions. Compared with the Monte Carlo methods, our method yields significant savings in computation time. Compared with other moment propagation methods, the recursive expressions in our method can be implemented by specifying rate constants and stoichiometries, without having to derive or solve any differential equations. Whereas other moment propagation methods with similar accuracy have been demonstrated for a few species, we demonstrate our method for complex biochemical systems with hundreds of species.