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.     

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.     

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.     

A hybrid fluctuating hydrodynamics and kinetic Monte Carlo method for modeling chemically-powered nanoscale motion

We present a stochastic multiscale method for modeling heterogeneous catalysis at the nanoscale. The system is decomposed into the fluid domain and the catalyst-fluid interface. We implemented the fluctuating hydrodynamics framework to model the diffusion of the chemical species in the fluid domain, and the chemical master equation to describe the catalytic activity at the interface. The coupling between the domains occurs simultaneously. Using a simple one-dimensional (1D) linear model, we showed that the predictions of our scheme are in excellent agreement with deterministic simulations. The method was specifically developed to model the spatially asymmetric catalysis on the surface of self-propelled nanoswimmers. Numerical simulations showed that our approach can estimate the uncertainty in the swimming velocity resulting from inherent stochastic nature of the chemical reactions at the catalytic interface. Although the method has been applied to simple 1D and 2D models, it can be generalized to handle different geometries and more sophisticated chemical reactions. Therefore, it can serve as a practical mathematical tool to study how the efficiency of chemically powered nanomachines is affected by the interplay between structural complexity, nonlinear reactivity, and nonequilibrium fluctuations.     

Hybrid stochastic simulations of intracellular reaction–diffusion systems

With the observation that stochasticity is important in biological systems, chemical kinetics have begun to receive wider interest. While the use of Monte Carlo discrete event simulations most accurately capture the variability of molecular species, they become computationally costly for complex reaction–diffusion systems with large populations of molecules. On the other hand, continuous time models are computationally efficient but they fail to capture any variability in the molecular species. In this study a hybrid stochastic approach is introduced for simulating reaction–diffusion systems. We developed an adaptive partitioning strategy in which processes with high frequency are simulated with deterministic rate-based equations, and those with low frequency using the exact stochastic algorithm of Gillespie. Therefore the stochastic behavior of cellular pathways is preserved while being able to apply it to large populations of molecules. We describe our method and demonstrate its accuracy and efficiency compared with the Gillespie algorithm for two different systems. First, a model of intracellular viral kinetics with two steady states and second, a compartmental model of the postsynaptic spine head for studying the dynamics of Ca+2 and NMDA receptors.     

Cellular automata models of kinetically and thermodynamically controlled reactions

Cellular automata simulations of the competition between kinetically controlled and thermodynamically controlled products of a reaction are described. The simulations are based on a stochastic first‐order cellular automata model described previously [ 20 ] and demonstrate an alternative to the traditional approach to such problems that relies on solution of a set of coupled differential rate equations. Unlike the traditional approach, the cellular automata models are applicable to finite numbers of elements and yield statistical information on the fluctuations to be expected in such finite cases. The usual deterministic solutions appear as limiting cases involving either very large numbers of reacting ingredients or a large number of trials for smaller sets of ingredients. © 2000 John Wiley & Sons, Inc. Int J Chem Kinet 32: 529–534, 2000     

A kinetic Monte Carlo simulation study of inositol 1,4,5-trisphosphate receptor (IP3R) calcium release channel

Most of the previously theoretical studies about the stochastic nature of the IP3R calcium release channel gating use the chemical master equation (CME) approach. Because of the limitations of this approach we have used a stochastic simulation algorithm (SSA) presented by Gillespie. A single subunit of De Young-Keizer (DYK) model was simulated using Gillespie algorithm. The model has been considered in its complete form with eight states. We investigate the conditions which affect the open state of the model. Calcium concentrations were the subject of fluctuation in the previous works while in this study the population of the states is the subject of stochastic fluctuations. We found out that decreasing open probability is a function of Ca(2+) concentration in fast time domain, while in slow time domain it is a function of IP3 concentration. Studying the population of each state shows a time dependent reaction pattern in fast and medium time domains (10(-4) and 10(-3)s). In this pattern the state of X(010) has a determinative role in selecting the open state path. Also, intensity and frequency of fluctuations and Ca(2+) inhibitions have been studied. The results indicate that Gillespie algorithm can be a better choice for studying such systems, without using any approximation or elimination while having acceptable accuracy. In comparison with the chemical master equation, Gillespie algorithm is also provides a wide area for studying biological systems from other points of view.     

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.     

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.     

Efficient stochastic simulation of simultaneous reaction and diffusion in a gas-liquid interface

A stochastic simulation of simultaneous reaction and diffusion is proposed for the gas-liquid interface formed in the surface of a gas bubble within a liquid. The interface between a carbon dioxide bubble and an aqueous solution of calcium hydroxide was simulated as an application example, taken from the integrated production of calcium carbonate. First Gillespie’s stochastic simulation algorithm was applied in separate reaction and diffusion simulations. The results from these simulations were consistent with deterministic solutions based on differential equations. However it was observed that stochastic diffusion simulations are extremely slow. The sampling of diffusion events was accelerated applying a group molecule transfer scheme based on the binomial distribution function. Simulations of the reaction-diffusion in the gas-liquid interface based on the standard Gillespie’s stochastic algorithm were also slow. However the application of the binomial distribution function scheme allowed to compute the concentration profiles in the gas-liquid interface in a fraction of the time required with the standard Gillespie’s stochastic algorithm.     

A contribution to non-equilibrium chemical kinetics. II. Strongly non-equilibrium “hot” reactions in liquids and solids

The method proposed in part I for non-equilibrium chemical kinetics is applied to processes provoked by non-equilibrium assemblies of energetic particles in liquids and solids. The movement of such an energetic particle belonging to a certain energy group is considered as a stochastic process when the direction of the velocity is changed stochastically at each step. On the ground of this consideration a simplified model of such a process is introduced: the stochastic movement of a particle is replaced by the deterministic movement of the corresponding quasi-particle having parameters determined through corresponding averages of the stochastic process. By use of this model, group constants of kinetic equations of our abovementioned work were expressed through parameters of microscopic processes in solids and liquids, and systems of non-equilibrium chemical kinetics' equations were written for different case. The proposed approach also permits us to consider the non-equilibrium of the crystalline lattice created by energetic particles. “Hot spot” reactions were considered as an example and a method to distinguish between direct and “hot spot” reactions was indicated. The proposed approach and obtained kinetic equations can be applied to recoil atoms (ions), fission products, hot particles produced in radiation chemistry, photochemistry, by laser beams, flash-photolysis etc. The destruction of the crystalline lattice by laser beams can also be considered by use of these equations.     

Developing Itô stochastic differential equation models for neuronal signal transduction pathways

Mathematical modeling and simulation of dynamic biochemical systems are receiving considerable attention due to the increasing availability of experimental knowledge of complex intracellular functions. In addition to deterministic approaches, several stochastic approaches have been developed for simulating the time-series behavior of biochemical systems. The problem with stochastic approaches, however, is the larger computational time compared to deterministic approaches. It is therefore necessary to study alternative ways to incorporate stochasticity and to seek approaches that reduce the computational time needed for simulations, yet preserve the characteristic behavior of the system in question. In this work, we develop a computational framework based on the It? stochastic differential equations for neuronal signal transduction networks. There are several different ways to incorporate stochasticity into deterministic differential equation models and to obtain It? stochastic differential equations. Two of the developed models are found most suitable for stochastic modeling of neuronal signal transduction. The best models give stable responses which means that the variances of the responses with time are not increasing and negative concentrations are avoided. We also make a comparative analysis of different kinds of stochastic approaches, that is the It? stochastic differential equations, the chemical Langevin equation, and the Gillespie stochastic simulation algorithm. Different kinds of stochastic approaches can be used to produce similar responses for the neuronal protein kinase C signal transduction pathway. The fine details of the responses vary slightly, depending on the approach and the parameter values. However, when simulating great numbers of chemical species, the Gillespie algorithm is computationally several orders of magnitude slower than the It? stochastic differential equations and the chemical Langevin equation. Furthermore, the chemical Langevin equation produces negative concentrations. The It? stochastic differential equations developed in this work are shown to overcome the problem of obtaining negative values.     

First-principles molecular dynamics simulations of condensed-phase V-type nerve agent reaction pathways and energy barriers

Computational studies of condensed-phase chemical reactions are challenging in part because of complexities in understanding the effects of the solvent environment on the reacting chemical species. Such studies are further complicated due to the demanding computational resources required to implement high-level ab initio quantum chemical methods when considering the solvent explicitly. Here, we use first-principles molecular dynamics simulations to examine condensed-phase decontamination reactions of V-type nerve agents in an explicit aqueous solvent. Our results include a detailed study of hydrolysis, base-hydrolysis, and nucleophilic oxidation of both VX and R-VX, as well as their protonated counterparts (i.e., VXH(+) and R-VXH(+)). The decontamination mechanisms and chemical reaction energy barriers, as determined from our simulations, are found to be in good agreement with experiment. The results demonstrate the applicability of using such simulations to assist in understanding new decontamination technologies or other applications that require computational screening of condensed-phase chemical reaction mechanisms.     

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.     

COAST: Controllable approximative stochastic reaction algorithm

We present an approximative algorithm for stochastic simulations of chemical reaction systems, called COAST, based on three different modeling levels: for small numbers of particles an exact stochastic model; for intermediate numbers an approximative, but computationally more efficient stochastic model based on discrete Gaussian distributions; and for large numbers the deterministic reaction kinetics. In every simulation time step, the subdivision of the reaction channels into the three different modeling levels is done automatically, where all approximations applied can be controlled by a single error parameter for which an appropriate value can easily be found. Test simulations show that the results of COAST simulations agree well with the outcomes of exact algorithms; however, the asymptotic run times of COAST are asymptotically proportional to smaller powers of the particle numbers than exact algorithms.     

一混沌模型系综模拟研究

用随机模拟方法研究了化学混沌模型的介观动力学。对该混沌模型的系综模拟发现,在这种不稳定运动中存在强烈的内部涨落,然而由于混沌运动整体上的稳定性,使得系综中的代表点被限制在混沌吸引子上,并且单个代表点形成的随机轨道很好地保持了确定性混沌吸引子的基本特征。     

A constrained approach to multiscale stochastic simulation of chemically reacting systems

Stochastic simulation of coupled chemical reactions is often computationally intensive, especially if a chemical system contains reactions occurring on different time scales. In this paper, we introduce a multiscale methodology suitable to address this problem, assuming that the evolution of the slow species in the system is well approximated by a Langevin process. It is based on the conditional stochastic simulation algorithm (CSSA) which samples from the conditional distribution of the suitably defined fast variables, given values for the slow variables. In the constrained multiscale algorithm (CMA) a single realization of the CSSA is then used for each value of the slow variable to approximate the effective drift and diffusion terms, in a similar manner to the constrained mean-force computations in other applications such as molecular dynamics. We then show how using the ensuing Fokker-Planck equation approximation, we can in turn approximate average switching times in stochastic chemical systems.