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.     

Efficient Stochastic Simulation Algorithm for Chemically Reacting Systems Based on Support Vector Regression

The stochastic simulation algorithm （SSA） accurately depicts spatially homogeneous wellstirred chemically reacting systems 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, a twin support vector regression based stochastic simulations algorithm （TS^3A） is proposed by combining the twin support vector regression and SSA, the former is a well-known robust regression method in machine learning. Numerical results indicate that this proposed algorithm can be applied to a wide range of chemically reacting systems and obtain significant improvements on efficiency and accuracy with fewer simulating runs over the existing methods.     

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.     

A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks

The time evolution of species concentrations in biochemical reaction networks is often modeled using the stochastic simulation algorithm (SSA) [Gillespie, J. Phys. Chem. 81, 2340 (1977)]. The computational cost of the original SSA scaled linearly with the number of reactions in the network. Gibson and Bruck developed a logarithmic scaling version of the SSA which uses a priority queue or binary tree for more efficient reaction selection [Gibson and Bruck, J. Phys. Chem. A 104, 1876 (2000)]. More generally, this problem is one of dynamic discrete random variate generation which finds many uses in kinetic Monte Carlo and discrete event simulation. We present here a constant-time algorithm, whose cost is independent of the number of reactions, enabled by a slightly more complex underlying data structure. While applicable to kinetic Monte Carlo simulations in general, we describe the algorithm in the context of biochemical simulations and demonstrate its competitive performance on small- and medium-size networks, as well as its superior constant-time performance on very large networks, which are becoming necessary to represent the increasing complexity of biochemical data for pathways that mediate cell function.     

A data-integrated method for analyzing stochastic biochemical networks

Variability and fluctuations among genetically identical cells under uniform experimental conditions stem from the stochastic nature of biochemical reactions. Understanding network function for endogenous biological systems or designing robust synthetic genetic circuits requires accounting for and analyzing this variability. Stochasticity in biological networks is usually represented using a continuous-time discrete-state Markov formalism, where the chemical master equation (CME) and its kinetic Monte Carlo equivalent, the stochastic simulation algorithm (SSA), are used. These two representations are computationally intractable for many realistic biological problems. Fitting parameters in the context of these stochastic models is particularly challenging and has not been accomplished for any but very simple systems. In this work, we propose that moment equations derived from the CME, when treated appropriately in terms of higher order moment contributions, represent a computationally efficient framework for estimating the kinetic rate constants of stochastic network models and subsequent analysis of their dynamics. To do so, we present a practical data-derived moment closure method for these equations. In contrast to previous work, this method does not rely on any assumptions about the shape of the stochastic distributions or a functional relationship among their moments. We use this method to analyze a stochastic model of a biological oscillator and demonstrate its accuracy through excellent agreement with CME/SSA calculations. By coupling this moment-closure method with a parameter search procedure, we further demonstrate how a model's kinetic parameters can be iteratively determined in order to fit measured distribution data.     

Michaelis-Menten speeds up tau-leaping under a wide range of conditions

This paper examines the benefits of Michaelis-Menten model reduction techniques in stochastic tau-leaping simulations. Results show that although the conditions for the validity of the reductions for tau-leaping remain the same as those for the stochastic simulation algorithm (SSA), the reductions result in a substantial speed-up for tau-leaping under a different range of conditions than they do for SSA. The reason of this discrepancy is that the time steps for SSA and for tau-leaping are determined by different properties of system dynamics.     

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.     

State-dependent doubly weighted stochastic simulation algorithm for automatic characterization of stochastic biochemical rare events

In recent years there has been substantial growth in the development of algorithms for characterizing rare events in stochastic biochemical systems. Two such algorithms, the state-dependent weighted stochastic simulation algorithm (swSSA) and the doubly weighted SSA (dwSSA) are extensions of the weighted SSA (wSSA) by H. Kuwahara and I. Mura [J. Chem. Phys. 129, 165101 (2008)]. The swSSA substantially reduces estimator variance by implementing system state-dependent importance sampling (IS) parameters, but lacks an automatic parameter identification strategy. In contrast, the dwSSA provides for the automatic determination of state-independent IS parameters, thus it is inefficient for systems whose states vary widely in time. We present a novel modification of the dwSSA--the state-dependent doubly weighted SSA (sdwSSA)--that combines the strengths of the swSSA and the dwSSA without inheriting their weaknesses. The sdwSSA automatically computes state-dependent IS parameters via the multilevel cross-entropy method. We apply the method to three examples: a reversible isomerization process, a yeast polarization model, and a lac operon model. Our results demonstrate that the sdwSSA offers substantial improvements over previous methods in terms of both accuracy and efficiency.     

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.     

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.     

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.     

Efficient formulation of the stochastic simulation algorithm for chemically reacting systems

In this paper we examine the different formulations of Gillespie's stochastic simulation algorithm (SSA) [D. Gillespie, J. Phys. Chem. 81, 2340 (1977)] with respect to computational efficiency, and propose an optimization to improve the efficiency of the direct method. Based on careful timing studies and an analysis of the time-consuming operations, we conclude that for most practical problems the optimized direct method is the most efficient formulation of SSA. This is in contrast to the widely held belief that Gibson and Bruck's next reaction method [M. Gibson and J. Bruck, J. Phys. Chem. A 104, 1876 (2000)] is the best way to implement the SSA for large systems. Our analysis explains the source of the discrepancy.     

The finite state projection algorithm for the solution of the chemical master equation

This article introduces the finite state projection (FSP) method for use in the stochastic analysis of chemically reacting systems. One can describe the chemical populations of such systems with probability density vectors that evolve according to a set of linear ordinary differential equations known as the chemical master equation (CME). Unlike Monte Carlo methods such as the stochastic simulation algorithm (SSA) or tau leaping, the FSP directly solves or approximates the solution of the CME. If the CME describes a system that has a finite number of distinct population vectors, the FSP method provides an exact analytical solution. When an infinite or extremely large number of population variations is possible, the state space can be truncated, and the FSP method provides a certificate of accuracy for how closely the truncated space approximation matches the true solution. The proposed FSP algorithm systematically increases the projection space in order to meet prespecified tolerance in the total probability density error. For any system in which a sufficiently accurate FSP exists, the FSP algorithm is shown to converge in a finite number of steps. The FSP is utilized to solve two examples taken from the field of systems biology, and comparisons are made between the FSP, the SSA, and tau leaping algorithms. In both examples, the FSP outperforms the SSA in terms of accuracy as well as computational efficiency. Furthermore, due to very small molecular counts in these particular examples, the FSP also performs far more effectively than tau leaping methods.     

Exact on-lattice stochastic reaction-diffusion simulations using partial-propensity methods

Stochastic reaction-diffusion systems frequently exhibit behavior that is not predicted by deterministic simulation models. Stochastic simulation methods, however, are computationally expensive. We present a more efficient stochastic reaction-diffusion simulation algorithm that samples realizations from the exact solution of the reaction-diffusion master equation. The present algorithm, called partial-propensity stochastic reaction-diffusion (PSRD) method, uses an on-lattice discretization of the reaction-diffusion system and relies on partial-propensity methods for computational efficiency. We describe the algorithm in detail, provide a theoretical analysis of its computational cost, and demonstrate its computational performance in benchmarks. We then illustrate the application of PSRD to two- and three-dimensional pattern-forming Gray-Scott systems, highlighting the role of intrinsic noise in these systems.     

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.     

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.     

Dissipative particle dynamics at isothermal, isobaric, isoenergetic, and isoenthalpic conditions using Shardlow-like splitting algorithms

Numerical integration schemes based upon the Shardlow-splitting algorithm (SSA) are presented for dissipative particle dynamics (DPD) approaches at various fixed conditions, including a constant-enthalpy (DPD-H) method that is developed by combining the equations-of-motion for a barostat with the equations-of-motion for the constant-energy (DPD-E) method. The DPD-H variant is developed for both a deterministic (Hoover) and stochastic (Langevin) barostat, where a barostat temperature is defined to satisfy the fluctuation-dissipation theorem for the Langevin barostat. For each variant, the Shardlow-splitting algorithm is formulated for both a velocity-Verlet scheme and an implicit scheme, where the velocity-Verlet scheme consistently performed better. The application of the Shardlow-splitting algorithm is particularly critical for the DPD-E and DPD-H variants, since it allows more temporally practical simulations to be carried out. The equivalence of the DPD variants is verified using both a standard DPD fluid model and a coarse-grain solid model. For both models, the DPD-E and DPD-H variants are further verified by instantaneously heating a slab of particles in the simulation cell, and subsequent monitoring of the evolution of the corresponding thermodynamic variables as the system approaches an equilibrated state while maintaining their respective constant-energy and constant-enthalpy conditions. The original SSA formulated for systems of equal-mass particles has been extended to systems of unequal-mass particles. The Fokker-Planck equation and derivations of the fluctuation-dissipation theorem for each DPD variant are also included for completeness.     

"All possible steps" approach to the accelerated use of Gillespie's algorithm

Many physical and biological processes are stochastic in nature. Computational models and simulations of such processes are a mathematical and computational challenge. The basic stochastic simulation algorithm was published by Gillespie about three decades ago [J. Phys. Chem. 81, 2340 (1977)]. Since then, intensive work has been done to make the algorithm more efficient in terms of running time. All accelerated versions of the algorithm are aimed at minimizing the running time required to produce a stochastic trajectory in state space. In these simulations, a necessary condition for reliable statistics is averaging over a large number of simulations. In this study the author presents a new accelerating approach which does not alter the stochastic algorithm, but reduces the number of required runs. By analysis of collected data the author demonstrates high precision levels with fewer simulations. Moreover, the suggested approach provides a good estimation of statistical error, which may serve as a tool for determining the number of required runs.