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

TLDR

Chemical systems can be modeled either as continuous differential equations or as discrete stochastic processes, but most complex systems with many species and reaction channels lack analytic solutions and exact simulation methods struggle to scale. This paper introduces the Next Reaction Method, an exact algorithm that uses a single random number per event and scales with the logarithm of the number of reactions. The method is extended to handle time‑dependent rate constants and non‑Markov processes and is applied to the lambda phage lysis/lysogeny decision circuit. Its performance on this biological application was compared with a standard method and an optimized variant, demonstrating competitive efficiency.

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 time-dependent 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.