A reduction method for multiple time scale stochastic reaction networks with non-unique equilibrium probability

We develop a reduction method for general closed multiple time scale stochastic reaction networks for which the fast subsystem may have non-unique equilibrium probability. We obtain a reduced ODE system with nonhomogeneous terms whose solutions can approximate the solutions of the full system accurately. We then apply this reduction method to general linear network and nonlinear networks for which the state diagram can be constructed. We illustrate the accuracy of the reduction method by comparing computational results of the full systems with the reduced ODE systems for several examples. Finally, we show how the reduction method may be extended to three or more time scale reaction networks.     

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.     

Reduction and solution of the chemical master equation using time scale separation and finite state projection

The dynamics of chemical reaction networks often takes place on widely differing time scales--from the order of nanoseconds to the order of several days. This is particularly true for gene regulatory networks, which are modeled by chemical kinetics. Multiple time scales in mathematical models often lead to serious computational difficulties, such as numerical stiffness in the case of differential equations or excessively redundant Monte Carlo simulations in the case of stochastic processes. We present a model reduction method for study of stochastic chemical kinetic systems that takes advantage of multiple time scales. The method applies to finite projections of the chemical master equation and allows for effective time scale separation of the system dynamics. We implement this method in a novel numerical algorithm that exploits the time scale separation to achieve model order reductions while enabling error checking and control. We illustrate the efficiency of our method in several examples motivated by recent developments in gene regulatory networks.     

A probability generating function method for stochastic reaction networks

In this paper we present a probability generating function (PGF) approach for analyzing stochastic reaction networks. The master equation of the network can be converted to a partial differential equation for PGF. Using power series expansion of PGF and Padé approximation, we develop numerical schemes for finding probability distributions as well as first and second moments. We show numerical accuracy of the method by simulating chemical reaction examples such as a binding-unbinding reaction, an enzyme-substrate model, Goldbeter-Koshland ultrasensitive switch model, and G(2)/M transition model.     

Two classes of quasi-steady-state model reductions for stochastic kinetics

The quasi-steady-state approximation (QSSA) is a model reduction technique used to remove highly reactive species from deterministic models of reaction mechanisms. In many reaction networks the highly reactive intermediates (QSSA species) have populations small enough to require a stochastic representation. In this work we apply singular perturbation analysis to remove the QSSA species from the chemical master equation for two classes of problems. The first class occurs in reaction networks where all the species have small populations and the QSSA species sample zero the majority of the time. The perturbation analysis provides a reduced master equation in which the highly reactive species can sample only zero, and are effectively removed from the model. The reduced master equation can be sampled with the Gillespie algorithm. This first stochastic QSSA reduction is applied to several example reaction mechanisms (including Michaelis-Menten kinetics) [Biochem. Z. 49, 333 (1913)]. A general framework for applying the first QSSA reduction technique to new reaction mechanisms is derived. The second class of QSSA model reductions is derived for reaction networks where non-QSSA species have large populations and QSSA species numbers are small and stochastic. We derive this second QSSA reduction from a combination of singular perturbation analysis and the Omega expansion. In some cases the reduced mechanisms and reaction rates from these two stochastic QSSA models and the classical deterministic QSSA reduction are equivalent; however, this is not usually the case.     

Efficient stochastic simulations of complex reaction networks on surfaces

Surfaces serve as highly efficient catalysts for a vast variety of chemical reactions. Typically, such surface reactions involve billions of molecules which diffuse and react over macroscopic areas. Therefore, stochastic fluctuations are negligible and the reaction rates can be evaluated using rate equations, which are based on the mean-field approximation. However, in case that the surface is partitioned into a large number of disconnected microscopic domains, the number of reactants in each domain becomes small and it strongly fluctuates. This is, in fact, the situation in the interstellar medium, where some crucial reactions take place on the surfaces of microscopic dust grains. In this case rate equations fail and the simulation of surface reactions requires stochastic methods such as the master equation. However, in the case of complex reaction networks, the master equation becomes infeasible because the number of equations proliferates exponentially. To solve this problem, we introduce a stochastic method based on moment equations. In this method the number of equations is dramatically reduced to just one equation for each reactive species and one equation for each reaction. Moreover, the equations can be easily constructed using a diagrammatic approach. We demonstrate the method for a set of astrophysically relevant networks of increasing complexity. It is expected to be applicable in many other contexts in which problems that exhibit analogous structure appear, such as surface catalysis in nanoscale systems, aerosol chemistry in stratospheric clouds, and genetic networks in cells.     

An analytical approach to solutions of master equations for stochastic nonlinear reactions

In this paper we consider a class of nonlinear reactions which are important in stochastic reaction networks. We find the exact solution of the chemical master equation for a class of irreversible and reversible nonlinear reactions. We also present the explicit form of the equilibrium probability solution of the reactions. The results can be used for analyzing stochastic dynamics of important reactions such as binding/unbinding reaction and protein dimerization.     

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.     

Multiple‐Well,multiple‐path unimolecular reaction systems. I. MultiWell computer program suite

Unimolecular reaction systems in which multiple isomers undergo simultaneous reactions via multiple decomposition reactions and multiple isomerization reactions are of fundamental interest in chemical kinetics. The computer program suite described here can be used to treat such coupled systems, including the effects of collisional energy transfer (weak collisions). The program suite consists of MultiWell, which solves the internal energy master equation for complex unimolecular reactions systems; DenSum, which calculates sums and densities of states by an exact‐count method; MomInert, which calculates external principal moments of inertia and internal rotation reduced moments of inertia; and Thermo, which calculates equilibrium constants and other thermodynamics quantities. MultiWell utilizes a hybrid master equation approach, which performs like an energy‐grained master equation at low energies and a continuum master equation in the vibrational quasicontinuum. An adaptation of Gillespie's exact stochastic method is used for the solution. The codes are designed for ease of use. Details are presented of various methods for treating weak collisions with virtually any desired collision step‐size distribution and for utilizing RRKM theory for specific unimolecular rate constants. © 2001 John Wiley & Sons, Inc. Int J Chem Kinet 33: 232–245, 2001     

A numerical characteristic method for probability generating functions on stochastic first-order reaction networks

We propose an efficient and accurate numerical scheme for solving probability generating functions arising in stochastic models of general first-order reaction networks by using the characteristic curves. A partial differential equation derived by a probability generating function is the transport equation with variable coefficients. We apply the idea of characteristics for the estimation of statistical measures, consisting of the mean, variance, and marginal probability. Estimation accuracy is obtained by the Newton formulas for the finite difference and time accuracy is obtained by applying the fourth order Runge–Kutta scheme for the characteristic curve and the Simpson method for the integration on the curve. We apply our proposed method to motivating biological examples and show the accuracy by comparing simulation results from the characteristic method with those from the stochastic simulation algorithm.     

Path ensembles and path sampling in nonequilibrium stochastic systems

Markovian models based on the stochastic master equation are often encountered in single molecule dynamics, reaction networks, and nonequilibrium problems in chemistry, physics, and biology. An efficient and convenient method to simulate these systems is the kinetic Monte Carlo algorithm which generates continuous-time stochastic trajectories. We discuss an alternative simulation method based on sampling of stochastic paths. Utilizing known probabilities of stochastic paths, it is possible to apply Metropolis Monte Carlo in path space to generate a desired ensemble of stochastic paths. The method is a generalization of the path sampling idea to stochastic dynamics, and is especially suited for the analysis of rare paths which are not often produced in the standard kinetic Monte Carlo procedure. Two generic examples are presented to illustrate the methodology.     

Lattice microbes: High‐performance stochastic simulation method for the reaction‐diffusion master equation

Spatial stochastic simulation is a valuable technique for studying reactions in biological systems. With the availability of high‐performance computing (HPC), the method is poised to allow integration of data from structural, single‐molecule and biochemical studies into coherent computational models of cells. Here, we introduce the Lattice Microbes software package for simulating such cell models on HPC systems. The software performs either well‐stirred or spatially resolved stochastic simulations with approximated cytoplasmic crowding in a fast and efficient manner. Our new algorithm efficiently samples the reaction‐diffusion master equation using NVIDIA graphics processing units and is shown to be two orders of magnitude faster than exact sampling for large systems while maintaining an accuracy of ～ 0.1%. Display of cell models and animation of reaction trajectories involving millions of molecules is facilitated using a plug‐in to the popular VMD visualization platform. The Lattice Microbes software is open source and available for download at http://www.scs.illinois.edu/schulten/lm © 2012 Wiley Periodicals, Inc.     

The stochastic quasi-steady-state assumption: reducing the model but not the noise

Highly reactive species at small copy numbers play an important role in many biological reaction networks. We have described previously how these species can be removed from reaction networks using stochastic quasi-steady-state singular perturbation analysis (sQSPA). In this paper we apply sQSPA to three published biological models: the pap operon regulation, a biochemical oscillator, and an intracellular viral infection. These examples demonstrate three different potential benefits of sQSPA. First, rare state probabilities can be accurately estimated from simulation. Second, the method typically results in fewer and better scaled parameters that can be more readily estimated from experiments. Finally, the simulation time can be significantly reduced without sacrificing the accuracy of the solution.     

Steady-state fluctuations of a genetic feedback loop: an exact solution

Genetic feedback loops in cells break detailed balance and involve bimolecular reactions; hence, exact solutions revealing the nature of the stochastic fluctuations in these loops are lacking. We here consider the master equation for a gene regulatory feedback loop: a gene produces protein which then binds to the promoter of the same gene and regulates its expression. The protein degrades in its free and bound forms. This network breaks detailed balance and involves a single bimolecular reaction step. We provide an exact solution of the steady-state master equation for arbitrary values of the parameters, and present simplified solutions for a number of special cases. The full parametric dependence of the analytical non-equilibrium steady-state probability distribution is verified by direct numerical solution of the master equations. For the case where the degradation rate of bound and free protein is the same, our solution is at variance with a previous claim of an exact solution [J. E. M. Hornos, D. Schultz, G. C. P. Innocentini, J. Wang, A. M. Walczak, J. N. Onuchic, and P. G. Wolynes, Phys. Rev. E 72, 051907 (2005), and subsequent studies]. We show explicitly that this is due to an unphysical formulation of the underlying master equation in those studies.     

Master equation approach to vibrational energy transfer between two diatomic molecules

In this paper a semiclassical non-markovian master equation is derived. We begin by using the well-known tetradic form of the Liouville equation for a reduced density operator. By projecting the diagonal matrix elements of the operator, we obtain an infinite-order master equation. This equation is then applied in the lowest-order approximation to collinear collisions between the diatomic molecules: H₂H₂, N₂N₂ and Cl₂Cl₂. With an assumed form of the interaction potential for such a problem we have also derived an analytical expression for the V—V transition probabilities. They are then calculated over a wide range of velocities of the colliding molecules and compared with exact semiclassical ones. An excellent agreement of the results is found for small velocities (i.e. υ ≈ 10⁴ cm/s). For larger values of υ (≈ 10⁵ cm/s) the results obtained from the master equation approach agree with the exact ones only in the low-velocity range for light molecules and low oscillatory states.     

受高斯白噪声影响的有限化学反应体系的随机热力学——外噪声对化学反应体系影响的有效热力学量度

通过分析噪声对跃迁概率的不同影响,借助Novikov定理及MSR理论,建立了受多源白噪声影响的有限化学反应体系的有效主方程及有效熵平衡方程,导出非平衡定态时这类噪声体系熵产生的一般表达式,揭示涨落熵产生的统计内涵及噪声贡献,并针对非宏观量级的外噪声,借助扰动按分布参数分离法及有效主方程的Kramers-Moyal展开,进一步对简单加合性噪声建立了非平衡定态宏观稳定性判据的随机模拟,论证了噪声对化学反应体系定态稳定性的弱化作用.     

A variational approach to the stochastic aspects of cellular signal transduction

Cellular signaling networks have evolved to cope with intrinsic fluctuations, coming from the small numbers of constituents, and the environmental noise. Stochastic chemical kinetics equations govern the way biochemical networks process noisy signals. The essential difficulty associated with the master equation approach to solving the stochastic chemical kinetics problem is the enormous number of ordinary differential equations involved. In this work, we show how to achieve tremendous reduction in the dimensionality of specific reaction cascade dynamics by solving variationally an equivalent quantum field theoretic formulation of stochastic chemical kinetics. The present formulation avoids cumbersome commutator computations in the derivation of evolution equations, making the physical significance of the variational method more transparent. We propose novel time-dependent basis functions which work well over a wide range of rate parameters. We apply the new basis functions to describe stochastic signaling in several enzymatic cascades and compare the results so obtained with those from alternative solution techniques. The variational Ansatz gives probability distributions that agree well with the exact ones, even when fluctuations are large and discreteness and nonlinearity are important. A numerical implementation of our technique is many orders of magnitude more efficient computationally compared with the traditional Monte Carlo simulation algorithms or the Langevin simulations.     

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.     

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.