Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels

Abstract

There are two fundamental ways to view coupled systems of chemical equations: as continuous, represented by differential equations whose variables are concentrations, or as discrete, represented by stochastic processes whose variables are numbers of molecules. Although the former is by far more common, systems with very small numbers of molecules are important in some applications (e.g., in small biological cells or in surface processes). In both views, most complicated systems with multiple reaction channels and multiple chemical species cannot be solved analytically. There are exact numerical simulation methods to simulate trajectories of discrete, stochastic systems, (methods that are rigorously equivalent to the Master Equation approach) but these do not scale well to systems with many reaction pathways. This paper presents the Next Reaction Method, an exact algorithm to simulate coupled chemical reactions that is also efficient: it (a) uses only a single random number per simulation event, and (b) takes time proportional to the logarithm of the number of reactions, not to the number of reactions itself. The Next Reaction Method is extended to include timedependent rate constants and non-Markov processes and is applied to a sample application in biology (the lysis/lysogeny decision circuit of lambda phage). The performance of the Next Reaction Method on this application is compared with one standard method and an optimized version of that standard method.

Mass fluctuation kinetics: analysis and computation of equilibria and local dynamics

“…The tau-leaping algorithm, which approximates the SSA, alleviates this drawback and is presented in the next section. Yet other approaches to e¢ cient stochastic simulation may be found in the literature, see for example [4].…”

“…The thermodynamic limit of the system is de…ned as its limiting behavior with large molecule numbers, x(t), and system size, , i.e., thermodynamic limit 4 = lim fx(t); g!1 fsystemg (1.13)…”

An integrative simulation model linking major biochemical reactions of actin-polymerization to structural properties of actin filaments

“…In particular, the exact Gillespie's "direct" and "first reaction" and the "next reaction" algorithms [2]. The inexact "τ-leap" algorithm has been tested as well.…”

“…• class CSimpleCell is an abstract parent class, in which the numbers of molecules for reagents are stored and the methods of stochastic chemical reactions simulations are implemented (the "direct" method, the "first reaction" method [4], the "next reaction" method [2]). The functions of the class operate by a list of 21 reagents in the same way for both nSRF and nSRF models.…”

“…Several stochastic simulation strategies have been implemented and compared. More specifically, we developed new simulation schemes based on the updates of the Gillespie's "direct" method and its modifications: the "first reaction" method [34], the "next reaction" method of Gibson-Bruck [35] and the τ-leap algorithm [36]. The best results were obtained for the one generated on the base of the Gillespie's "direct" …”