Search

Generic open quantum systems are notoriously difficult to simulate unless one looks at specific regimes. In contrast, classical dissipative systems can often be effectively described by stochastic processes, which are generally less computationally expensive. Here, we use the paradigmatic case of a dissipative quantum oscillator to give a pedagogic introduction to the modeling of open quantum systems using quasiclassical methods, i.e., classical stochastic methods that use a "quantum" noise spectrum to capture the influence of the environment on the system. Such quasiclassical methods have the potential to offer insights into the impact of the quantum nature of the environment on the dynamics of the system of interest while still being computationally tractable. [ABSTRACT FROM AUTHOR]

A hybrid stochastic simulation method is developed to study H2–O2 auto-ignition at the microscale. Simulation results show that the discrete and stochastic characteristics of reaction collisions have notable impacts on the ignition process, particularly in the early stages when only a few radicals exist. The statistical properties of ignition delay time, which reflect the accumulated stochasticity during ignition, are obtained and analyzed for different initial temperatures and total molecular numbers. It is found that the average and standard deviation of ignition delay time increase as the total molecular number decreases, with this phenomenon being particularly pronounced near the crossover temperature. When the total molecular number is sufficiently small, the chain initiation reaction becomes crucial to the stochastic properties, as its average firing time exhibits an inverse proportionality to the total molecular number. As the total molecular number increases, the influence of other chain reactions intensifies, causing the power law relation between standard deviation and total molecular number to shift from −1 power to −0.5 power. Owing to different chain reaction paths for high- and low-temperature auto-ignition, the strongest relative fluctuation occurs near the crossover temperature. A theoretical equation for the standard deviation of ignition delay time is obtained based on dimensional analysis, giving excellent agreement with the simulation results in both high- and low-temperature modes. [ABSTRACT FROM AUTHOR]

Tau-leaping methods have recently been proposed for the acceleration of discrete stochastic simulation of chemically reacting systems. This paper considers the numerical stability of these methods. The concept of stochastic absolute stability is defined, discussed, and applied to the following leaping methods: the explicit tau, implicit tau, and trapezoidal tau.(C) 2004 American Institute of Physics.

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. (C) 2004 American Institute of Physics.

Computational models of reaction–diffusion systems involving low copy numbers or strongly heterogeneous molecular spatial distributions, such as those frequently found in cellular signaling pathways, require approaches that account for the stochastic dynamics of individual particles, as opposed to approaches representing them through their average concentrations. Efforts to remedy the high computational cost associated with particle-based stochastic approaches by taking advantage of Green's functions are hampered by the need to draw random numbers from complicated, and therefore costly, non-standard probability distributions to update particle positions. Here, we introduce an approach that permits the reconstruction of entire molecular trajectories, including bimolecular encounters, in retrospect, after a simulated time step, while avoiding inefficient draws from non-standard distributions. This means that highly accurate stochastic simulations can be performed for system sizes that would be prohibitively costly to simulate with conventional Green's function based methods. The algorithm applies equally well to one, two, and three dimensional systems and can be readily extended to include deterministic forces specified by an interaction potential, such as the Coulomb potential. [ABSTRACT FROM AUTHOR]

Understanding phononic heat transport processes in molecular junctions is a central issue in the developing field of nanoscale heat conduction. Here, we present a Langevin dynamics simulation framework to investigate heat transport processes in molecular junctions at and beyond the linear response regime and apply it to saturated and unsaturated linear hydrocarbon chains connecting two gold substrates. Thermal boundary conditions represented by Markovian noise and damping are filtered through several (up to four) gold layers to provide a realistic and controllable bath spectral density. Classical simulations using the full universal force field are compared with quantum calculations that use only the harmonic part of this field. The close agreement found at about room temperature between these very different calculations suggests that heat transport at such temperatures is dominated by lower frequency vibrations whose dynamics is described well by classical mechanics. The results obtained for alkanedithiol molecules connecting gold substrates agree with previous quantum calculations based on the Landauer formula and match recent experimental measurements [e.g., thermal conductance around 20 pW/K for alkanedithiols in single-molecule junctions (SMJs)]. Heat conductance simulations on polyynes of different lengths illuminate the effects of molecular conjugation on thermal transport. The difference between alkanes and polyynes is not large but correlates with the larger rigidity and stronger mode localization that characterize the polyyne structure. This computational approach has been recently used [R. Chen, I. Sharony, and A. Nitzan, J. Phys. Chem. Lett. 11, 4261–4268 (2020)] to unveil local atomic heat currents and phononic interference effect in aromatic-ring based SMJs. [ABSTRACT FROM AUTHOR]

This work considers the method of uniformization for continuous-time Markov chains in the context of chemical reaction networks. Previous work in the literature has shown that uniformization can be beneficial in the context of time-inhomogeneous models, such as chemical reaction networks incorporating extrinsic noise. This paper lays focus on the understanding of uniformization from the viewpoint of sample paths of chemical reaction networks. In particular, an efficient pathwise stochastic simulation algorithm for time-homogeneous models is presented which is complexity-wise equal to Gillespie's direct method. This new approach therefore enlarges the class of problems for which the uniformization approach forms a computationally attractive choice. Furthermore, as a new application of the uniformization method, we provide a novel variance reduction method for (raw) moment estimators of chemical reaction networks based upon the combination of stratification and uniformization. [ABSTRACT FROM AUTHOR]

The workhorse simulating the dissipative dynamics is mainly based on the harmonic bath potentials together with linear system-bath couplings, but a realistic bath always assumes anharmonicity. In this work, we extend the linear dissipation model to include quadratic couplings and suggest a stochastic simulation scheme for the anharmonic dissipation. We show that the non-Gaussian noises induced by the anharmonic bath can be rigorously constructed, and the resulting stochastic Liouville equation has the same form as that for the linear dissipation model. As a preliminary application, we use this stochastic method to investigate the vibration-induced symmetry breaking in two-level electronic systems and find that the characteristic function of the non-Gaussian noises determines the absorption and fluorescence spectra. [ABSTRACT FROM AUTHOR]

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. (C) 2005 American Institute of Physics.

Rare event probabilities play an important role in the understanding of the behavior of biochemical systems. Due to the intractability of the most natural Markov jump process representation of a system of interest, rare event probabilities are typically estimated using importance sampling. While the resulting algorithm is reasonably well developed, the problem of choosing a suitable importance density is far from straightforward. We therefore leverage recent developments on simulation of conditioned jump processes to propose an importance density that is simple to implement and requires no tuning. Our results demonstrate superior performance over some existing approaches. [ABSTRACT FROM AUTHOR]

The stochastic simulation algorithm has been used to generate exact trajectories of biochemical reaction networks. For each simulation step, the simulation selects a reaction and its firing time according to a probability that is proportional to the reaction propensity. We investigate in this paper new efficient formulations of the stochastic simulation algorithm to improve its computational efficiency.We examine the selection of the next reaction firing and reduce its computational cost by reusing the computation in the previous step. For biochemical reactions with delays, we present a new method for computing the firing time of the next reaction. The principle for computing the firing time of our approach is based on recycling of random numbers. Our newapproach for generating the firing time of the next reaction is not only computationally efficient but also easy to implement.We further analyze and reduce the number of propensity updates when a delayed reaction occurred. We demonstrate the applicability of our improvements by experimenting with concrete biological models. [ABSTRACT FROM AUTHOR]

The reaction-diffusion master equation (RDME) is a model that allows for efficient on-lattice simulation of spatially resolved stochastic chemical kinetics. Compared to off-lattice hard-sphere simulations with Brownian dynamics or Green's function reaction dynamics, the RDME can be orders of magnitude faster if the lattice spacing can be chosen coarse enough. However, strongly diffusion-controlled reactions mandate a very fine mesh resolution for acceptable accuracy. It is common that reactions in the same model differ in their degree of diffusion control and therefore require different degrees of mesh resolution. This renders mesoscopic simulation inefficient for systems with multiscale properties. Mesoscopic-microscopic hybrid methods address this problem by resolving the most challenging reactions with a microscale, off-lattice simulation. However, all methods to date require manual partitioning of a system, effectively limiting their usefulness as "black-box" simulation codes. In this paper, we propose a hybrid simulation algorithm with automatic system partitioning based on indirect a priori error estimates.We demonstrate the accuracy and efficiency of the method on models of diffusion-controlled networks in 3D. [ABSTRACT FROM AUTHOR]

In this paper, we present a novel hybrid method to simulate discrete stochastic reaction-diffusion models arising in biochemical signaling pathways. We study moderately stiff systems, for which we can partition each reaction or diffusion channel into either a slow or fast subset, based on its propensity. Numerical approaches missing this distinction are often limited with respect to computational run time or approximation quality. We design an approximate scheme that remedies these pitfalls by using a new blending strategy of the well-established inhomogeneous stochastic simulation algorithm and the tau-leaping simulation method. The advantages of our hybrid simulation algorithm are demonstrated on three benchmarking systems, with special focus on approximation accuracy and efficiency. [ABSTRACT FROM AUTHOR]

At the molecular level, biochemical processes are governed by random interactions between reactant molecules, and the dynamics of such systems are inherently stochastic. When the copy numbers of reactants are large, a deterministic description is adequate, but when they are small, such systems are often modeled as continuous-time Markov jump processes that can be described by the chemical master equation. Gillespie's Stochastic Simulation Algorithm (SSA) generates exact trajectories of these systems, but the amount of computational work required for each step of the original SSA is proportional to the number of reaction channels, leading to computational complexity that scales linearly with the problem size. The original SSA is therefore inefficient for large problems, which has prompted the development of several alternative formulations with improved scaling properties. We describe an exact SSA that uses a table data structure with event time binning to achieve constant computational complexity with respect to the number of reaction channels for weakly coupled reaction networks. We present a novel adaptive binning strategy and discuss optimal algorithm parameters. We compare the computational efficiency of the algorithm to existing methods and demonstrate excellent scaling for large problems. This method is well suited for generating exact trajectories of large weakly coupled models, including those that can be described by the reaction-diffusion master equation that arises from spatially discretized reaction-diffusion processes. [ABSTRACT FROM AUTHOR]

Over decades, the theoretical study of the quantum dissipative dynamics was mainly based on the linear dissipation model. The study of the nonlinear dissipative dynamics in condensed phases, where there exist an infinite number of bath modes, is extremely difficult even if not impossible. This work put forward a stochastic scheme for the simulation of the nonlinear dissipative dynamics. In the linear response regime, the second-order cumulant expansion becomes exact to reproduce the effect of the bath on the evolution of the reduced system. Consequently, a Hermitian stochastic Liouville equation is derived without explicit treatment of the bath. Stochastic simulations for an anharmonic model illustrate that the dynamics dissipated by anharmonic bath exhibits substantial difference on temperature dependence compared to that with the Caldeira-Leggett model. [ABSTRACT FROM AUTHOR]

We have developed a method for modeling spatial stochastic biochemical reactions in complex, threedimensional, and time-dependent domains using the reaction-diffusion master equation formalism. In particular, we look to address the fully coupled problems that arise in systems biology where the shape and mechanical properties of a cell are determined by the state of the biochemistry and vice versa. To validate our method and characterize the error involved, we compare our results for a carefully constructed test problem to those of a microscale implementation. We demonstrate the effectiveness of our method by simulating a model of polarization and shmoo formation during the mating of yeast. The method is generally applicable to problems in systems biology where biochemistry and mechanics are coupled, and spatial stochastic effects are critical. [ABSTRACT FROM AUTHOR]

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. [ABSTRACT FROM AUTHOR]

Stochastic simulation of reaction-diffusion systems enables the investigation of stochastic events arising from the small numbers and heterogeneous distribution of molecular species in biological cells. Stochastic variations in intracellular microdomains and in diffusional gradients play a significant part in the spatiotemporal activity and behavior of cells. Although an exact stochastic simulation that simulates every individual reaction and diffusion event gives a most accurate trajectory of the system's state over time, it can be too slow for many practical applications. We present an accelerated algorithm for discrete stochastic simulation of reaction-diffusion systems designed to improve the speed of simulation by reducing the number of time-steps required to complete a simulation run. This method is unique in that it employs two strategies that have not been incorporated in existing spatial stochastic simulation algorithms. First, diffusive transfers between neighboring subvolumes are based on concentration gradients. This treatment necessitates sampling of only the net or observed diffusion events from higher to lower concentration gradients rather than sampling all diffusion events regardless of local concentration gradients. Second, we extend the non-negative Poisson tau-leaping method that was originally developed for speeding up nonspatial or homogeneous stochastic simulation algorithms. This method calculates each leap time in a unified step for both reaction and diffusion processes while satisfying the leap condition that the propensities do not change appreciably during the leap and ensuring that leaping does not cause molecular populations to become negative. Numerical results are presented that illustrate the improvement in simulation speed achieved by incorporating these two new strategies. [ABSTRACT FROM AUTHOR]

We present the partial-propensity stochastic simulation algorithm with composition-rejection sampling (PSSA-CR). It is an exact formulation of the stochastic simulation algorithm (SSA) for well-stirred systems of coupled chemical reactions. The new formulation is a partial-propensity variant [R. Ramaswamy, N. González-Segredo, and I. F. Sbalzarini, J. Chem. Phys. 130, 244104 (2009)] of the composition- rejection SSA [A. Slepoy, A. P. Thompson, and S. J. Plimpton, J. Chem. Phys. 128, 205101 (2008)]. The computational cost of this new formulation is bounded by a constant for weakly coupled reaction networks, and it increases at most linearly with the number of chemical species for strongly coupled reaction networks. PSSA-CR thus combines the advantages of partial-propensity methods and the composition-rejection SSA, providing favorable scaling of the computational cost for all classes of reaction networks. [ABSTRACT FROM AUTHOR]

Gillespie’s direct method (DM) [D. Gillespie, J. Chem. Phys. 81, 2340 (1977)] for exact stochastic simulation of chemical reaction systems has been widely adopted. It is easy to implement but requires large computation for relatively large systems. Recently, two more efficient methods, next reaction method (NRM) [M. A. Gibson and J. Bruck, J. Phys. Chem. A 105, 1876 (2000)] and optimized DM (ODM) [Y. Cao et al., J. Chem. Phys. 121, 4059 (2004)], have been developed to improve simulation speed. It has been demonstrated that the ODM is the state-of-the-art most efficient method for exact stochastic simulation of most practical reaction systems. In this paper, we first develop an exact stochastic simulation algorithm named ODMK that is more efficient than the ODM. We then develop an approximate method named K-skip method to further accelerate simulation. Using two chemical reaction systems, we demonstrate that our ODMK and K-skip method can save 20%–30% and 70%–80% simulation time, respectively, comparing to the ODM. We also show that our ODMK and K-skip method provide almost the same simulation accuracy as the ODM. [ABSTRACT FROM AUTHOR]

The weighted stochastic simulation algorithm (wSSA) recently introduced by Kuwahara and Mura [J. Chem. Phys. 129, 165101 (2008)] is an innovative variation on the stochastic simulation algorithm (SSA). It enables one to estimate, with much less computational effort than was previously thought possible using a Monte Carlo simulation procedure, the probability that a specified event will occur in a chemically reacting system within a specified time when that probability is very small. This paper presents some procedural extensions to the wSSA that enhance its effectiveness in practical applications. The paper also attempts to clarify some theoretical issues connected with the wSSA, including its connection to first passage time theory and its relation to the SSA. [ABSTRACT FROM AUTHOR]

An exact method for stochastic simulation of chemical reaction networks, which accelerates the stochastic simulation algorithm (SSA), is proposed. The present “ER-leap” algorithm is derived from analytic upper and lower bounds on the multireaction probabilities sampled by SSA, together with rejection sampling and an adaptive multiplicity for reactions. The algorithm is tested on a number of well-quantified reaction networks and is found experimentally to be very accurate on test problems including a chaotic reaction network. At the same time ER-leap offers a substantial speedup over SSA with a simulation time proportional to the 2/3 power of the number of reaction events in a Galton–Watson process. [ABSTRACT FROM AUTHOR]

The Inhomogeneous Stochastic Simulation Algorithm (ISSA) is a variant of the stochastic simulation algorithm in which the spatially inhomogeneous volume of the system is divided into homogeneous subvolumes, and the chemical reactions in those subvolumes are augmented by diffusive transfers of molecules between adjacent subvolumes. The ISSA can be prohibitively slow when the system is such that diffusive transfers occur much more frequently than chemical reactions. In this paper we present the Multinomial Simulation Algorithm (MSA), which is designed to, on the one hand, outperform the ISSA when diffusive transfer events outnumber reaction events, and on the other, to handle small reactant populations with greater accuracy than deterministic-stochastic hybrid algorithms. The MSA treats reactions in the usual ISSA fashion, but uses appropriately conditioned binomial random variables for representing the net numbers of molecules diffusing from any given subvolume to a neighbor within a prescribed distance. Simulation results illustrate the benefits of the algorithm. [ABSTRACT FROM AUTHOR]

The slow-scale stochastic simulation algorithm (ssSSA) proposed by Cao et al. [J. Chem. Phys. 122, 014116 (2005)] and the nested stochastic simulation algorithm (nSSA) proposed by E et al. [J. Chem. Phys. 123, 194107 (2005)] are closely related approximate simulation procedures aimed at speeding up the stochastic simulation of stiff chemical systems, i.e., systems that evolve through fast and slow dynamical modes with the fast mode being stable. This Comment aims to clarify some misconceptions that have arisen over the relationship between the ssSSA and the nSSA as regards both their theoretical foundations and their practical implementations. [ABSTRACT FROM AUTHOR]

Leap methods are very promising for accelerating stochastic simulation of a well stirred chemically reacting system, while providing acceptable simulation accuracy. In Gillespie’s τ-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. [ABSTRACT FROM AUTHOR]

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. [ABSTRACT FROM AUTHOR]

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. [ABSTRACT FROM AUTHOR]

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. © 2005 American Institute of Physics. [ABSTRACT FROM AUTHOR]

A new algorithm, 'HiER-leap' (hierarchical exact reaction-leaping), is derived which improves on the computational properties of the ER-leap algorithm for exact accelerated simulation of stochastic chemical kinetics. Unlike ER-leap, HiER-leap utilizes a hierarchical or divide-and-conquer organization of reaction channels into tightly coupled 'blocks' and is thereby able to speed up systems with many reaction channels. Like ER-leap, HiER-leap is based on the use of upper and lower bounds on the reaction propensities to define a rejection sampling algorithm with inexpensive early rejection and acceptance steps. But in HiER-leap, large portions of intra-block sampling may be done in parallel. An accept/reject step is used to synchronize across blocks. This method scales well when many reaction channels are present and has desirable asymptotic properties. The algorithm is exact, parallelizable and achieves a significant speedup over the stochastic simulation algorithm and ER-leap on certain problems. This algorithm offers a potentially important step towards efficient in silico modeling of entire organisms. [ABSTRACT FROM AUTHOR]

Numerical simulation of the time evolution of a spatially homogeneous chemical system is always of great interest. Gillespie first developed the exact stochastic simulation algorithm (SSA), which is accurate but time-consuming. Recently, many approximate schemes of the SSA are proposed to speed up simulation. Presented here is the N-leap method, which guarantees the validity of the leap condition and at the same time keeps the efficiency. In many cases, N-leap has better performance than the widely-used τ-leap method. The details of the N-leap method are described and several examples are presented to show its validity. [ABSTRACT FROM AUTHOR]

In recent years, computer simulations have become increasingly useful when trying to understand the complex dynamics of biochemical networks, particularly in stochastic systems. In such situations stochastic simulation is vital in gaining an understanding of the inherent stochasticity present, as these models are rarely analytically tractable. However, a stochastic approach can be computationally prohibitive for many models. A number of approximations have been proposed that aim to speed up stochastic simulations. However, the majority of these approaches are fundamentally serial in terms of central processing unit (CPU) usage. In this paper, we propose a novel simulation algorithm that utilises the potential of multi-core machines. This algorithm partitions the model into smaller sub-models. These sub-models are then simulated, in parallel, on separate CPUs. We demonstrate that this method is accurate and can speed-up the simulation by a factor proportional to the number of processors available. [ABSTRACT FROM AUTHOR]

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. [ABSTRACT FROM AUTHOR]

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. [ABSTRACT FROM AUTHOR]

Several real-world systems, such as gene expression networks in biological cells, contain coupled chemical reactions with a time delay between reaction initiation and completion. The non-Markovian kinetics of such reaction networks can be exactly simulated using the delay stochastic simulation algorithm (dSSA). The computational cost of dSSA scales with the total number of reactions in the network. We reduce this cost to scale at most with the smaller number of species by using the concept of partial reaction propensities. The resulting delay partial-propensity direct method (dPDM) is an exact dSSA formulation for well-stirred systems of coupled chemical reactions with delays. We detail dPDM and present a theoretical analysis of its computational cost. Furthermore, we demonstrate the implications of the theoretical cost analysis in two prototypical benchmark applications. The dPDM formulation is shown to be particularly efficient for strongly coupled reaction networks, where the number of reactions is much larger than the number of species. [ABSTRACT FROM AUTHOR]

The weighted stochastic simulation algorithm (wSSA) was developed by Kuwahara and Mura [J. Chem. Phys. 129, 165101 (2008)] to efficiently estimate the probabilities of rare events in discrete stochastic systems. The wSSA uses importance sampling to enhance the statistical accuracy in the estimation of the probability of the rare event. The original algorithm biases the reaction selection step with a fixed importance sampling parameter. In this paper, we introduce a novel method where the biasing parameter is state-dependent. The new method features improved accuracy, efficiency, and robustness. [ABSTRACT FROM AUTHOR]

The Chemical Langevin Equation (CLE), which is a stochastic differential equation driven by a multidimensional Wiener process, acts as a bridge between the discrete stochastic simulation algorithm and the deterministic reaction rate equation when simulating (bio)chemical kinetics. The CLE model is valid in the regime where molecular populations are abundant enough to assume their concentrations change continuously, but stochastic fluctuations still play a major role. The contribution of this work is that we observe and explore that the CLE is not a single equation, but a parametric family of equations, all of which give the same finite-dimensional distribution of the variables. On the theoretical side, we prove that as many Wiener processes are sufficient to formulate the CLE as there are independent variables in the equation, which is just the rank of the stoichiometric matrix. On the practical side, we show that in the case where there are m1 pairs of reversible reactions and m2 irreversible reactions there is another, simple formulation of the CLE with only m1+m2 Wiener processes, whereas the standard approach uses 2m1+m2. We demonstrate that there are considerable computational savings when using this latter formulation. Such transformations of the CLE do not cause a loss of accuracy and are therefore distinct from model reduction techniques. We illustrate our findings by considering alternative formulations of the CLE for a human ether a-go-go related gene ion channel model and the Goldbeter–Koshland switch. [ABSTRACT FROM AUTHOR]

We introduce an alternative formulation of the exact stochastic simulation algorithm (SSA) for sampling trajectories of the chemical master equation for a well-stirred system of coupled chemical reactions. Our formulation is based on factored-out, partial reaction propensities. This novel exact SSA, called the partial-propensity direct method (PDM), is highly efficient and has a computational cost that scales at most linearly with the number of chemical species, irrespective of the degree of coupling of the reaction network. In addition, we propose a sorting variant, SPDM, which is especially efficient for multiscale reaction networks. [ABSTRACT FROM AUTHOR]

In robust biological systems, wide deviations from highly controlled normal behavior may be rare, yet they may result in catastrophic complications. While in silico analysis has gained an appreciation as a tool to offer insights into system-level properties of biological systems, analysis of such rare events provides a particularly challenging computational problem. This paper proposes an efficient stochastic simulation method to analyze rare events in biochemical systems. Our new approach can substantially increase the frequency of the rare events of interest by appropriately manipulating the underlying probability measure of the system, allowing high-precision results to be obtained with substantially fewer simulation runs than the conventional direct Monte Carlo simulation. Here, we show the algorithm of our new approach, and we apply it to the analysis of rare deviant transitions of two systems, resulting in several orders of magnitude speedup in generating high-precision estimates compared with the conventional Monte Carlo simulation. [ABSTRACT FROM AUTHOR]

Recently the application of the quasi-steady-state approximation (QSSA) to the stochastic simulation algorithm (SSA) was suggested for the purpose of speeding up stochastic simulations of chemical systems that involve both relatively fast and slow chemical reactions [Rao and Arkin, J. Chem. Phys. 118, 4999 (2003)] and further work has led to the nested and slow-scale SSA. Improved numerical efficiency is obtained by respecting the vastly different time scales characterizing the system and then by advancing only the slow reactions exactly, based on a suitable approximation to the fast reactions. We considerably extend these works by applying the QSSA to numerical methods for the direct solution of the chemical master equation (CME) and, in particular, to the finite state projection algorithm [Munsky and Khammash, J. Chem. Phys. 124, 044104 (2006)], in conjunction with Krylov methods. In addition, we point out some important connections to the literature on the (deterministic) total QSSA (tQSSA) and place the stochastic analogue of the QSSA within the more general framework of aggregation of Markov processes. We demonstrate the new methods on four examples: Michaelis–Menten enzyme kinetics, double phosphorylation, the Goldbeter–Koshland switch, and the mitogen activated protein kinase cascade. Overall, we report dramatic improvements by applying the tQSSA to the CME solver. [ABSTRACT FROM AUTHOR]

The τ-leap method first developed by Gillespie [D. T. Gillespie, J. Chem. Phys. 115, 1716 (2001)] can significantly speed up stochastic simulation of certain chemically reacting systems with acceptable losses in accuracy. Recently, several improved τ-leap methods, including the binomial, multinomial, and modified τ-leap methods, have been developed. However, in all these τ-leap methods, the mean of the number of times, Km, that the mth reaction channel fires during a leap is not equal to the true mean. Therefore, all existing τ-leap methods produce biased simulation results, which limit the simulation accuracy and speed. In this paper, we analyze the mean of Km based on the chemical master equation. Using this analytical result, we develop unbiased Poisson and binomial τ-leap methods. Moreover, we analyze the variance of Km, and then develop an unbiased Poisson/Gaussian/binomial τ-leap method to correct the errors in both the mean and variance of Km. Simulation results demonstrate that our unbiased τ-leap method can significantly improve simulation accuracy without sacrificing speed. [ABSTRACT FROM AUTHOR]

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 τ-leap spatial stochastic simulation algorithm, in terms of efficiency and accuracy, in scenarios of both molecular homogeneity and heterogeneity. [ABSTRACT FROM AUTHOR]

Gillespie’s exact stochastic simulation algorithm (SSA) [J. Phys. Chem. 81, 2350 (1977)] has been widely used to simulate the stochastic dynamics of chemically reacting systems. In this algorithm, it is assumed that all reactions occur instantly. While this is true in many cases, it is also possible that some chemical reactions, such as gene transcription and translation in living cells, take certain time to finish after they are initiated. Thus, the product of such reactions will emerge after certain delays. Apparently, Gillespie’s SSA is not an exact algorithm for chemical reaction systems with delays. In this paper, the author develops an exact SSA for chemical reaction systems with delays, based upon the same fundamental premise of stochastic kinetics used by Gillespie in the development of his SSA. He then shows that an algorithm modified from Gillespie’s SSA by Barrio et al. [PLOS Comput. Biol. 2, 1017 (2006)] is also an exact SSA for chemical reaction systems with delays, but it needs to generate more random variables than the author’s algorithm. [ABSTRACT FROM AUTHOR]

A novel algorithm is proposed for the acceleration of the exact stochastic simulation algorithm by a predefined number of reaction firings (R-leaping) that may occur across several reaction channels. In the present approach, the numbers of reaction firings are correlated binomial distributions and the sampling procedure is independent of any permutation of the reaction channels. This enables the algorithm to efficiently handle large systems with disparate rates, providing substantial computational savings in certain cases. Several mechanisms for controlling the accuracy and the appearance of negative species are described. The advantages and drawbacks of R-leaping are assessed by simulations on a number of benchmark problems and the results are discussed in comparison with established methods. [ABSTRACT FROM AUTHOR]

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. [ABSTRACT FROM AUTHOR]

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. [ABSTRACT FROM AUTHOR]

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. © 2004 American Institute of Physics. [ABSTRACT FROM AUTHOR]

In numerically simulating the time evolution of a well-stirred chemically reacting system, the recently introduced “tau-leaping” procedure attempts to accelerate the exact stochastic simulation algorithm by using a special Poisson approximation to leap over sequences of noncritical reaction events. Presented here is an improved procedure for determining the maximum leap size for a specified degree of accuracy. © 2003 American Institute of Physics. [ABSTRACT FROM AUTHOR]

We propose a stochastic simulation for solving the Boltzmann equation for the Lorentz gas and develop a formalism to rigorously show the equivalence between the simulation and the Boltzmann equation. We also propose several simple simulations which offer a substantial improvement over the Boltzmann equation and discuss the simulations in the context of repeated ring equations. [ABSTRACT FROM AUTHOR]

The results of stochastic simulations of growth and evaporation of small clusters in vapor are reported. Energy dependent growth rates are determined from the monomer-cluster collision rate and decay rates are found from a detailed balance, with the equilibrium size and energy distribution of clusters calculated using the capillarity approximation and the equilibrium vapor pressure. These rates are used in simulations of two-dimensional random walks in size and energy space to determine the fraction of clusters in supersaturated vapor of size (imin+1) that reach a size imax. By assuming that clusters of size imin are in equilibrium, this fraction can be related to the nonisothermal nucleation rate. The simulated rates show good agreement with the previously published analytical results. In the absence of an inert carrier gas, the nonisothermal nucleation rates are typically between 1% and 5% of the isothermal rates. [ABSTRACT FROM AUTHOR]

The master equation of a lattice gas reaction tracks the probability of visiting all spatial configurations. The large number of unique spatial configurations on a lattice renders master equation simulations infeasible for even small lattices. In this work, a reduced master equation is derived for the probability distribution of the coverages in the infinite diffusion limit. This derivation justifies the widely used assumption that the adlayer is in equilibrium for the current coverages and temperature when all reactants are highly mobile. Given the reduced master equation, two novel and efficient simulation methods of lattice gas reactions in the infinite diffusion limit are derived. The first method involves solving the reduced master equation directly for small lattices, which is intractable in configuration space. The second method involves reducing the master equation further in the large lattice limit to a set of differential equations that tracks only the species coverages. Solution of the reduced master equation and differential equations requires information that can be obtained through short, diffusion-only kinetic Monte Carlo simulation runs at each coverage. These simulations need to be run only once because the data can be stored and used for simulations with any set of kinetic parameters, gas-phase concentrations, and initial conditions. An idealized CO oxidation reaction mechanism with strong lateral interactions is used as an example system for demonstrating the reduced master equation and deterministic simulation techniques. [ABSTRACT FROM AUTHOR]