Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions

The dynamical solution of a well-mixed, nonlinear stochastic chemical kinetic system, described by the Master equation, may be exactly computed using the stochastic simulation algorithm. However, because the computational cost scales with the number of reaction occurrences, systems with one or more "fast" reactions become costly to simulate. This paper describes a hybrid stochastic method that partitions the system into subsets of fast and slow reactions, approximates the fast reactions as a continuous Markov process, using a chemical Langevin equation, and accurately describes the slow dynamics using the integral form of the "Next Reaction" variant of the stochastic simulation algorithm. The key innovation of this method is its mechanism of efficiently monitoring the occurrences of slow, discrete events while simultaneously simulating the dynamics of a continuous, stochastic or deterministic process. In addition, by introducing an approximation in which multiple slow reactions may occur within a time step of the numerical integration of the chemical Langevin equation, the hybrid stochastic method performs much faster with only a marginal decrease in accuracy. Multiple examples, including a biological pulse generator and a large-scale system benchmark, are simulated using the exact and proposed hybrid methods as well as, for comparison, a previous hybrid stochastic method. Probability distributions of the solutions are compared and the weak errors of the first two moments are computed. In general, these hybrid methods may be applied to the simulation of the dynamics of a system described by stochastic differential, ordinary differential, and Master equations.     

Overcoming stiffness in stochastic simulation stemming from partial equilibrium: a multiscale Monte Carlo algorithm

In this paper the problem of stiffness in stochastic simulation of singularly perturbed systems is discussed. Such stiffness arises often from partial equilibrium or quasi-steady-state type of conditions. A multiscale Monte Carlo method is discussed that first assesses whether partial equilibrium is established using a simple criterion. The exact stochastic simulation algorithm (SSA) is next employed to sample among fast reactions over short time intervals (microscopic time steps) in order to compute numerically the proper probability distribution function for sampling the slow reactions. Subsequently, the SSA is used to sample among slow reactions and advance the time by large (macroscopic) time steps. Numerical examples indicate that not only long times can be simulated but also fluctuations are properly captured and substantial computational savings result.     

On the origins of approximations for stochastic chemical kinetics

This paper considers the derivation of approximations for stochastic chemical kinetics governed by the discrete master equation. Here, the concepts of (1) partitioning on the basis of fast and slow reactions as opposed to fast and slow species and (2) conditional probability densities are used to derive approximate, partitioned master equations, which are Markovian in nature, from the original master equation. Under different conditions dictated by relaxation time arguments, such approximations give rise to both the equilibrium and hybrid (deterministic or Langevin equations coupled with discrete stochastic simulation) approximations previously reported. In addition, the derivation points out several weaknesses in previous justifications of both the hybrid and equilibrium systems and demonstrates the connection between the original and approximate master equations. Two simple examples illustrate situations in which these two approximate methods are applicable and demonstrate the two methods' efficiencies.     

Stiffness detection and reduction in discrete stochastic simulation of biochemical systems

Typical multiscale biochemical models contain fast-scale and slow-scale reactions, where "fast" reactions fire much more frequently than "slow" ones. This feature often causes stiffness in discrete stochastic simulation methods such as Gillespie's algorithm and the Tau-Leaping method leading to inefficient simulation. This paper proposes a new strategy to automatically detect stiffness and identify species that cause stiffness for the Tau-Leaping method, as well as two stiffness reduction methods. Numerical results on a stiff decaying dimerization model and a heat shock protein regulation model demonstrate the efficiency and accuracy of the proposed methods for multiscale biochemical systems.     

Nested stochastic simulation algorithm for chemical kinetic systems with disparate rates

An efficient simulation algorithm for chemical kinetic systems with disparate rates is proposed. This new algorithm is quite general, and it amounts to a simple and seamless modification of the classical stochastic simulation algorithm (SSA), also known as the Gillespie [J. Comput. Phys. 22, 403 (1976); J. Phys. Chem. 81, 2340 (1977)] algorithm. The basic idea is to use an outer SSA to simulate the slow processes with rates computed from an inner SSA which simulates the fast reactions. Averaging theorems for Markov processes can be used to identify the fast and slow variables in the system as well as the effective dynamics over the slow time scale, even though the algorithm itself does not rely on such information. This nested SSA can be easily generalized to systems with more than two separated time scales. Convergence and efficiency of the algorithm are discussed using the established error estimates and illustrated through examples.     

The slow-scale stochastic simulation algorithm

Reactions in real chemical systems often take place on vastly different time scales, with "fast" reaction channels firing very much more frequently than "slow" ones. These firings will be interdependent if, as is usually the case, the fast and slow reactions involve some of the same species. An exact stochastic simulation of such a system will necessarily spend most of its time simulating the more numerous fast reaction events. This is a frustratingly inefficient allocation of computational effort when dynamical stiffness is present, since in that case a fast reaction event will be of much less importance to the system's evolution than will a slow reaction event. For such situations, this paper develops a systematic approximate theory that allows one to stochastically advance the system in time by simulating the firings of only the slow reaction events. Developing an effective strategy to implement this theory poses some challenges, but as is illustrated here for two simple systems, when those challenges can be overcome, very substantial increases in simulation speed can be realized.     

Look before you leap: a confidence-based method for selecting species criticality while avoiding negative populations in τ-leaping

The stochastic simulation algorithm was introduced by Gillespie and in a different form by Kurtz. There have been many attempts at accelerating the algorithm without deviating from the behavior of the simulated system. The crux of the explicit τ-leaping procedure is the use of Poisson random variables to approximate the number of occurrences of each type of reaction event during a carefully selected time period, τ. This method is acceptable providing the leap condition, that no propensity function changes "significantly" during any time-step, is met. Using this method there is a possibility that species numbers can, artificially, become negative. Several recent papers have demonstrated methods that avoid this situation. One such method classifies, as critical, those reactions in danger of sending species populations negative. At most, one of these critical reactions is allowed to occur in the next time-step. We argue that the criticality of a reactant species and its dependent reaction channels should be related to the probability of the species number becoming negative. This way only reactions that, if fired, produce a high probability of driving a reactant population negative are labeled critical. The number of firings of more reaction channels can be approximated using Poisson random variables thus speeding up the simulation while maintaining the accuracy. In implementing this revised method of criticality selection we make use of the probability distribution from which the random variable describing the change in species number is drawn. We give several numerical examples to demonstrate the effectiveness of our new method.     

Binomial tau-leap spatial stochastic simulation algorithm for applications in chemical kinetics

In cell biology, cell signaling pathway problems are often tackled with deterministic temporal models, well mixed stochastic simulators, and/or hybrid methods. But, in fact, three dimensional stochastic spatial modeling of reactions happening inside the cell is needed in order to fully understand these cell signaling pathways. This is because noise effects, low molecular concentrations, and spatial heterogeneity can all affect the cellular dynamics. However, there are ways in which important effects can be accounted without going to the extent of using highly resolved spatial simulators (such as single-particle software), hence reducing the overall computation time significantly. We present a new coarse grained modified version of the next subvolume method that allows the user to consider both diffusion and reaction events in relatively long simulation time spans as compared with the original method and other commonly used fully stochastic computational methods. Benchmarking of the simulation algorithm was performed through comparison with the next subvolume method and well mixed models (MATLAB), as well as stochastic particle reaction and transport simulations (CHEMCELL, Sandia National Laboratories). Additionally, we construct a model based on a set of chemical reactions in the epidermal growth factor receptor pathway. For this particular application and a bistable chemical system example, we analyze and outline the advantages of our presented binomial tau-leap spatial stochastic simulation algorithm, in terms of efficiency and accuracy, in scenarios of both molecular homogeneity and heterogeneity.     

K-leap method for accelerating stochastic simulation of coupled chemical reactions

Leap methods are very promising for accelerating stochastic simulation of a well stirred chemically reacting system, while providing acceptable simulation accuracy. In Gillespie's tau-leap method [D. Gillespie, J. Phys. Chem. 115, 1716 (2001)], the number of firings of each reaction channel during a leap is a Poisson random variable, whose sample values are unbounded. This may cause large changes in the populations of certain molecular species during a leap, thereby violating the leap condition. In this paper, we develop an alternative leap method called the K-leap method, in which we constrain the total number of reactions occurring during a leap to be a number K calculated from the leap condition. As the number of firings of each reaction channel during a leap is upper bounded by a properly chosen number, our K-leap method can better satisfy the leap condition, thereby improving simulation accuracy. Since the exact stochastic simulation algorithm (SSA) is a special case of our K-leap method when K=1, our K-leap method can naturally change from the exact SSA to an approximate leap method during simulation, whenever the leap condition allows to do so.     

Hybrid modeling and simulation of stochastic effects on progression through the eukaryotic cell cycle

The eukaryotic cell cycle is regulated by a complicated chemical reaction network. Although many deterministic models have been proposed, stochastic models are desired to capture noise in the cell resulting from low numbers of critical species. However, converting a deterministic model into one that accurately captures stochastic effects can result in a complex model that is hard to build and expensive to simulate. In this paper, we first apply a hybrid (mixed deterministic and stochastic) simulation method to such a stochastic model. With proper partitioning of reactions between deterministic and stochastic simulation methods, the hybrid method generates the same primary characteristics and the same level of noise as Gillespie's stochastic simulation algorithm, but with better efficiency. By studying the results generated by various partitionings of reactions, we developed a new strategy for hybrid stochastic modeling of the cell cycle. The new approach is not limited to using mass-action rate laws. Numerical experiments demonstrate that our approach is consistent with characteristics of noisy cell cycle progression, and yields cell cycle statistics in accord with experimental observations.     

Stochastic hybrid modeling of intracellular calcium dynamics

Deterministic models of biochemical processes at the subcellular level might become inadequate when a cascade of chemical reactions is induced by a few molecules. Inherent randomness of such phenomena calls for the use of stochastic simulations. However, being computationally intensive, such simulations become infeasible for large and complex reaction networks. To improve their computational efficiency in handling these networks, we present a hybrid approach, in which slow reactions and fluxes are handled through exact stochastic simulation and their fast counterparts are treated partially deterministically through chemical Langevin equation. The classification of reactions as fast or slow is accompanied by the assumption that in the time-scale of fast reactions, slow reactions do not occur and hence do not affect the probability of the state. Our new approach also handles reactions with complex rate expressions such as Michaelis-Menten kinetics. Fluxes which cannot be modeled explicitly through reactions, such as flux of Ca(2+) from endoplasmic reticulum to the cytosol through inositol 1,4,5-trisphosphate receptor channels, are handled deterministically. The proposed hybrid algorithm is used to model the regulation of the dynamics of cytosolic calcium ions in mouse macrophage RAW 264.7 cells. At relatively large number of molecules, the response characteristics obtained with the stochastic and deterministic simulations coincide, which validates our approach in the limit of large numbers. At low doses, the response characteristics of some key chemical species, such as levels of cytosolic calcium, predicted with stochastic simulations, differ quantitatively from their deterministic counterparts. These observations are ubiquitous throughout dose response, sensitivity, and gene-knockdown response analyses. While the relative differences between the peak-heights of the cytosolic [Ca(2+)] time-courses obtained from stochastic (mean of 16 realizations) and deterministic simulations are merely 1%-4% for most perturbations, it is specially sensitive to levels of G(βγ) (relative difference as large as 90% at very low G(βγ)).     

Accelerated stochastic simulation of the stiff enzyme-substrate reaction

The enzyme-catalyzed conversion of a substrate into a product is a common reaction motif in cellular chemical systems. In the three reactions that comprise this process, the intermediate enzyme-substrate complex is usually much more likely to decay into its original constituents than to produce a product molecule. This condition makes the reaction set mathematically "stiff." We show here how the simulation of this stiff reaction set can be dramatically speeded up relative to the standard stochastic simulation algorithm (SSA) by using a recently introduced procedure called the slow-scale SSA. The speedup occurs because the slow-scale SSA explicitly simulates only the relatively rare conversion reactions, skipping over occurrences of the other two less interesting but much more frequent reactions. We describe, explain, and illustrate this simulation procedure for the isolated enzyme-substrate reaction set, and then we show how the procedure extends to the more typical case in which the enzyme-substrate reactions occur together with other reactions and species. Finally, we explain the connection between this slow-scale SSA approach and the Michaelis-Menten [Biochem. Z. 49, 333 (1913)] formula, which has long been used in deterministic chemical kinetics to describe the enzyme-substrate reaction.     

An equation-free probabilistic steady-state approximation: dynamic application to the stochastic simulation of biochemical reaction networks

Stochastic chemical kinetics more accurately describes the dynamics of "small" chemical systems, such as biological cells. Many real systems contain dynamical stiffness, which causes the exact stochastic simulation algorithm or other kinetic Monte Carlo methods to spend the majority of their time executing frequently occurring reaction events. Previous methods have successfully applied a type of probabilistic steady-state approximation by deriving an evolution equation, such as the chemical master equation, for the relaxed fast dynamics and using the solution of that equation to determine the slow dynamics. However, because the solution of the chemical master equation is limited to small, carefully selected, or linear reaction networks, an alternate equation-free method would be highly useful. We present a probabilistic steady-state approximation that separates the time scales of an arbitrary reaction network, detects the convergence of a marginal distribution to a quasi-steady-state, directly samples the underlying distribution, and uses those samples to accurately predict the state of the system, including the effects of the slow dynamics, at future times. The numerical method produces an accurate solution of both the fast and slow reaction dynamics while, for stiff systems, reducing the computational time by orders of magnitude. The developed theory makes no approximations on the shape or form of the underlying steady-state distribution and only assumes that it is ergodic. We demonstrate the accuracy and efficiency of the method using multiple interesting examples, including a highly nonlinear protein-protein interaction network. The developed theory may be applied to any type of kinetic Monte Carlo simulation to more efficiently simulate dynamically stiff systems, including existing exact, approximate, or hybrid stochastic simulation techniques.     

The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior

A key to advancing the understanding of molecular biology in the post-genomic age is the development of accurate predictive models for genetic regulation, protein interaction, metabolism, and other biochemical processes. To facilitate model development, simulation algorithms must provide an accurate representation of the system, while performing the simulation in a reasonable amount of time. Gillespie's stochastic simulation algorithm (SSA) accurately depicts spatially homogeneous models with small populations of chemical species and properly represents noise, but it is often abandoned when modeling larger systems because of its computational complexity. In this work, we examine the performance of different versions of the SSA when applied to several biochemical models. Through our analysis, we discover that transient changes in reaction execution frequencies, which are typical of biochemical models with gene induction and repression, can dramatically affect simulator performance. To account for these shifts, we propose a new algorithm called the sorting direct method that maintains a loosely sorted order of the reactions as the simulation executes. Our measurements show that the sorting direct method performs favorably when compared to other well-known exact stochastic simulation algorithms.     

Abstract Next Subvolume Method: a logical process-based approach for spatial stochastic simulation of chemical reactions

The spatial stochastic simulation of biochemical systems requires significant calculation efforts. Parallel discrete-event simulation is a promising approach to accelerate the execution of simulation runs. However, achievable speedup depends on the parallelism inherent in the model. One of our goals is to explore this degree of parallelism in the Next Subvolume Method type simulations. Therefore we introduce the Abstract Next Subvolume Method, in which we decouple the model representation from the sequential simulation algorithms, and prove that state trajectories generated by its executions statistically accord with those generated by the Next Subvolume Method. The experimental performance analysis shows that optimistic synchronization algorithms, together with careful controls over the speculative execution, are necessary to achieve considerable speedup and scalability in parallel spatial stochastic simulation of chemical reactions. Our proposed method facilitates a flexible incorporation of different synchronization algorithms, and can be used to select the proper synchronization algorithm to achieve the efficient parallel simulation of chemical reactions.     

Solution reaction space Hamiltonian based on an electrostatic potential representation of solvent dynamics

Quantum chemical solvation models usually rely on the equilibrium solvation condition and is thus not immediately applicable to the study of nonequilibrium solvation dynamics, particularly those associated with chemical reactions. Here we address this problem by considering an effective Hamiltonian for solution-phase reactions based on an electrostatic potential (ESP) representation of solvent dynamics. In this approach a general ESP field of solvent is employed as collective solvent coordinate, and an effective Hamiltonian is constructed by treating both solute geometry and solvent ESP as dynamical variables. A harmonic bath is then attached onto the ESP variables in order to account for the stochastic nature of solvent dynamics. As an illustration we apply the above method to the proton transfer of a substituted phenol-amine complex in a polar solvent. The effective Hamiltonian is constructed by means of the reference interaction site model self-consistent field method (i.e., a type of quantum chemical solvation model), and a mixed quantum/classical simulation is performed in the space of solute geometry and solvent ESP. The results suggest that important dynamical features of proton transfer in solution can be captured by the present approach, including spontaneous fluctuations of solvent ESP that drives the proton from reactant to product potential wells.     

Hierarchical Order Parameters for Macromolecular Assembly Simulations I: Construction and Dynamical Properties of Order Parameters

Macromolecular assemblies often display a hierarchical organization of macromolecules or their sub-assemblies. To model this, we have formulated a space warping method that enables capturing overall macromolecular structure and dynamics via a set of coarse-grained order parameters (OPs). This article is the first of two describing the construction and computational implementation of an additional class of OPs that has built into them the hierarchical architecture of macromolecular assemblies. To accomplish this, first, the system is divided into subsystems, each of which is described via a representative set of OPs. Then, a global set of variables is constructed from these subsystem-centered OPs to capture overall system organization. Dynamical properties of the resulting OPs are compared to those of our previous nonhierarchical ones, and implied conceptual and computational advantages are discussed for a 100ns, 2 million atom solvated Human Papillomavirus-like particle simulation. In the second article, the hierarchical OPs are shown to enable a multiscale analysis that starts with the N-atom Liouville equation and yields rigorous Langevin equations of stochastic OP dynamics. The latter is demonstrated via a force-field based simulation algorithm that probes key structural transition pathways, simultaneously accounting for all-atom details and overall structure.     

A comparative study of the master equation approach and the discrete galerkin method for the simulation of free radical polymerization

The present article deals with the mathematical treatment of free radical polymerization reactions. As a typical example the synthesis of poly(methyl methacrylate) under realistic experimental conditions is investigated. Since the mathematical treatment of the kinetic rate equations raises severe numerical problems, alternative approaches are required. In this paper two of these methods, i.e. the discrete Galerkin method and the master equation approach, are compared. The discrete Galerkin method circumvents difficulties encountered by the direct integration of the kinetic rate equations but requires much a priori knowledge of the chemical reaction system. Within the framework of the master equation approach the polymerization reaction is regarded as a stochastic process. For the simulation of this stochastic process a modified algorithm is presented. The example of the polymerization of methyl methacrylate shows that the master equation approach is an efficient tool in the simulation of free radical polymerization reactions.     

Stochastic Simulation of Imperfect Mixing in Free Radical Polymerization

A modified algorithm for the stochastic simulation of chemical reactions subject to mass transfer limitation (imperfect mixing) is presented. This algorithm takes into account the mixing by diffusion of the reacting species between two consecutive reactions. The method is used to simulate the effect of mass transfer limitation in free-radical polymerization. Since this is a stiff reaction network, a hybrid stochastic-deterministic approach is considered. The hybrid stochastic algorithm under imperfect mixing (HSSA-IM) is applied to the bulk polymerization of methyl methacrylate up to high conversions. The accuracy of the algorithm relies on the precise determination of diffusion coefficients during the reaction.