Gillespie's Stochastic Simulation Research Articles

Binning for Efficient Stochastic Multiscale Particle Simulations

Gillespie's Stochastic Simulation Algorithm (SSA) is an exact procedure for simulating the evolution of a collection of discrete, interacting entities, such as coalescing aerosol particles or reacting chemical species. The high computational cost of SSA has motivated the development of more efficient variants, such as Tau-Leaping, which sacrifices the exactness of SSA. For models whose interacting entities can be characterized by a continuous parameter, such as a measure of size for aerosol particles, we analyze strategies for accelerating these algorithms by aggregating particles of similar size into bins. We show that for such models an appropriate binning strategy can dramatically enhance efficiency, and in particular can make SSA computationally competitive without sacrificing exactness. These strategies are especially effective for highly multiscale problems. We formulate binned versions of both the SSA and Tau-Leaping algorithms and analyze and demonstrate their performance.

Accelerating the Gillespie Exact Stochastic Simulation Algorithm Using Hybrid Parallel Execution on Graphics Processing Units

The Gillespie Stochastic Simulation Algorithm (GSSA) and its variants are cornerstone techniques to simulate reaction kinetics in situations where the concentration of the reactant is too low to allow deterministic techniques such as differential equations. The inherent limitations of the GSSA include the time required for executing a single run and the need for multiple runs for parameter sweep exercises due to the stochastic nature of the simulation. Even very efficient variants of GSSA are prohibitively expensive to compute and perform parameter sweeps. Here we present a novel variant of the exact GSSA that is amenable to acceleration by using graphics processing units (GPUs). We parallelize the execution of a single realization across threads in a warp (fine-grained parallelism). A warp is a collection of threads that are executed synchronously on a single multi-processor. Warps executing in parallel on different multi-processors (coarse-grained parallelism) simultaneously generate multiple trajectories. Novel data-structures and algorithms reduce memory traffic, which is the bottleneck in computing the GSSA. Our benchmarks show an 8×−120× performance gain over various state-of-the-art serial algorithms when simulating different types of models.

Efficient calculation of steady state probability distribution for stochastic biochemical reaction network

The Steady State (SS) probability distribution is an important quantity needed to characterize the steady state behavior of many stochastic biochemical networks. In this paper, we propose an efficient and accurate approach to calculating an approximate SS probability distribution from solution of the Chemical Master Equation (CME) under the assumption of the existence of a unique deterministic SS of the system. To find the approximate solution to the CME, a truncated state-space representation is used to reduce the state-space of the system and translate it to a finite dimension. The subsequent ill-posed eigenvalue problem of a linear system for the finite state-space can be converted to a well-posed system of linear equations and solved. The proposed strategy yields efficient and accurate estimation of noise in stochastic biochemical systems. To demonstrate the approach, we applied the method to characterize the noise behavior of a set of biochemical networks of ligand-receptor interactions for Bone Morphogenetic Protein (BMP) signaling. We found that recruitment of type II receptors during the receptor oligomerization by itself doesn't not tend to lower noise in receptor signaling, but regulation by a secreted co-factor may provide a substantial improvement in signaling relative to noise. The steady state probability approximation method shortened the time necessary to calculate the probability distributions compared to earlier approaches, such as Gillespie's Stochastic Simulation Algorithm (SSA) while maintaining high accuracy.

A Stochastic Model of DNA Fragments Rejoining

When cells are exposed to ionizing radiation, DNA damages in the form of single strand breaks (SSBs), double strand breaks (DSBs), base damage or their combinations are frequent events. It is known that the complexity and severity of DNA damage depends on the quality of radiation, and the microscopic dose deposited in small segments of DNA, which is often related to the linear transfer energy (LET) of the radiation. Experimental studies have suggested that under the same dose, high LET radiation induces more small DNA fragments than low-LET radiation, which affects Ku efficiently binding with DNA end and might be a main reason for high-LET radiation induced RBE [1] since DNA DSB is a major cause for radiation-induced cell death. In this work, we proposed a mathematical model of DNA fragments rejoining according to non-homologous end joining (NHEJ) mechanism. By conducting Gillespie's stochastic simulation, we found several factors that impact the efficiency of DNA fragments rejoining. Our results demonstrated that aberrant DNA damage repair can result predominantly from the occurrence of a spatial distribution of DSBs leading to short DNA fragments. Because of the low efficiency that short DNA fragments recruit repair protein and release the protein residue after fragments rejoining, Ku-dependent NHEJ is significantly interfered with short fragments. Overall, our work suggests that inhibiting the Ku-dependent NHEJ may significantly contribute to the increased efficiency for cell death and mutation observed for high LET radiation.

Efficient procedures for the numerical simulation of mid-size RNA kinetics

MotivationMethods for simulating the kinetic folding of RNAs by numerically solving the chemical master equation have been developed since the late 90's, notably the programs Kinfold and Treekin with Barriers that are available in the Vienna RNA package. Our goal is to formulate extensions to the algorithms used, starting from the Gillespie algorithm, that will allow numerical simulations of mid-size (~ 60–150 nt) RNA kinetics in some practical cases where numerous distributions of folding times are desired. These extensions can contribute to analyses and predictions of RNA folding in biologically significant problems.ResultsBy describing in a particular way the reduction of numerical simulations of RNA folding kinetics into the Gillespie stochastic simulation algorithm for chemical reactions, it is possible to formulate extensions to the basic algorithm that will exploit memoization and parallelism for efficient computations. These can be used to advance forward from the small examples demonstrated to larger examples of biological interest.SoftwareThe implementation that is described and used for the Gillespie algorithm is freely available by contacting the authors, noting that the efficient procedures suggested may also be applicable along with Vienna's Kinfold.

Reaction Event Counting Statistics of Biopolymer Reaction Systems with Dynamic Heterogeneity.

We investigate the reaction event counting statistics (RECS) of an elementary biopolymer reaction in which the rate coefficient is dependent on states of the biopolymer and the surrounding environment and discover a universal kinetic phase transition in the RECS of the reaction system with dynamic heterogeneity. From an exact analysis for a general model of elementary biopolymer reactions, we find that the variance in the number of reaction events is dependent on the square of the mean number of the reaction events when the size of measurement time is small on the relaxation time scale of rate coefficient fluctuations, which does not conform to renewal statistics. On the other hand, when the size of the measurement time interval is much greater than the relaxation time of rate coefficient fluctuations, the variance becomes linearly proportional to the mean reaction number in accordance with renewal statistics. Gillespie's stochastic simulation method is generalized for the reaction system with a rate coefficient fluctuation. The simulation results confirm the correctness of the analytic results for the time dependent mean and variance of the reaction event number distribution. On the basis of the obtained results, we propose a method of quantitative analysis for the reaction event counting statistics of reaction systems with rate coefficient fluctuations, which enables one to extract information about the magnitude and the relaxation times of the fluctuating reaction rate coefficient, without a bias that can be introduced by assuming a particular kinetic model of conformational dynamics and the conformation dependent reactivity. An exact relationship is established between a higher moment of the reaction event number distribution and the multitime correlation of the reaction rate for the reaction system with a nonequilibrium initial state distribution as well as for the system with the equilibrium initial state distribution.

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.

Error Bound for Piecewise Deterministic Processes Modeling Stochastic Reaction Systems

Biological processes involving the random interaction of $d$ species with integer particle numbers are often modeled by a Markov jump process on $\mathbb{N}_{0}^{d}$. Realizations of this process can, in principle, be generated with Gillespie's classical stochastic simulation algorithm, but for very reactive systems this method is usually inefficient. Hybrid models based on piecewise deterministic processes offer an attractive alternative which can decrease the simulation time considerably in applications where species with rather low particle numbers interact with very abundant species. We investigate the convergence of the hybrid model to the original one for a class of reaction systems with two distinct scales. Our main result is an error bound which states that, under suitable assumptions, the hybrid model approximates the marginal distribution of the discrete species and the conditional moments of the continuous species up to an error of $\mathcal{O}\!\left(M^{-1}\right)$, where $M$ is the scaling pa...

StochKit2: software for discrete stochastic simulation of biochemical systems with events

StochKit2 is the first major upgrade of the popular StochKit stochastic simulation software package. StochKit2 provides highly efficient implementations of several variants of Gillespie's stochastic simulation algorithm (SSA), and tau-leaping with automatic step size selection. StochKit2 features include automatic selection of the optimal SSA method based on model properties, event handling, and automatic parallelism on multicore architectures. The underlying structure of the code has been completely updated to provide a flexible framework for extending its functionality. StochKit2 runs on Linux/Unix, Mac OS X and Windows. It is freely available under GPL version 3 and can be downloaded from http://sourceforge.net/projects/stochkit/. petzold@engineering.ucsb.edu.

Computational study of noise in a large signal transduction network

BackgroundBiochemical systems are inherently noisy due to the discrete reaction events that occur in a random manner. Although noise is often perceived as a disturbing factor, the system might actually benefit from it. In order to understand the role of noise better, its quality must be studied in a quantitative manner. Computational analysis and modeling play an essential role in this demanding endeavor.ResultsWe implemented a large nonlinear signal transduction network combining protein kinase C, mitogen-activated protein kinase, phospholipase A2, and β isoform of phospholipase C networks. We simulated the network in 300 different cellular volumes using the exact Gillespie stochastic simulation algorithm and analyzed the results in both the time and frequency domain. In order to perform simulations in a reasonable time, we used modern parallel computing techniques. The analysis revealed that time and frequency domain characteristics depend on the system volume. The simulation results also indicated that there are several kinds of noise processes in the network, all of them representing different kinds of low-frequency fluctuations. In the simulations, the power of noise decreased on all frequencies when the system volume was increased.ConclusionsWe concluded that basic frequency domain techniques can be applied to the analysis of simulation results produced by the Gillespie stochastic simulation algorithm. This approach is suited not only to the study of fluctuations but also to the study of pure noise processes. Noise seems to have an important role in biochemical systems and its properties can be numerically studied by simulating the reacting system in different cellular volumes. Parallel computing techniques make it possible to run massive simulations in hundreds of volumes and, as a result, accurate statistics can be obtained from computational studies.

Equivalence of on-lattice stochastic chemical kinetics with the well-mixed chemical master equation in the limit of fast diffusion

Well-mixed and lattice-based descriptions of stochastic chemical kinetics have been extensively used in the literature. Realizations of the corresponding stochastic processes are obtained by the Gillespie stochastic simulation algorithm and lattice kinetic Monte Carlo algorithms, respectively. However, the two frameworks have remained disconnected. We show the equivalence of these frameworks whereby the stochastic lattice kinetics reduces to effective well-mixed kinetics in the limit of fast diffusion. In the latter, the lattice structure appears implicitly, as the lumped rate of bimolecular reactions depends on the number of neighbors of a site on the lattice. Moreover, we propose a mapping between the stochastic propensities and the deterministic rates of the well-mixed vessel and lattice dynamics that illustrates the hierarchy of models and the key parameters that enable model reduction.

STOCHSIMGPU: parallel stochastic simulation for the Systems Biology Toolbox 2 for MATLAB

The importance of stochasticity in biological systems is becoming increasingly recognized and the computational cost of biologically realistic stochastic simulations urgently requires development of efficient software. We present a new software tool STOCHSIMGPU that exploits graphics processing units (GPUs) for parallel stochastic simulations of biological/chemical reaction systems and show that significant gains in efficiency can be made. It is integrated into MATLAB and works with the Systems Biology Toolbox 2 (SBTOOLBOX2) for MATLAB. The GPU-based parallel implementation of the Gillespie stochastic simulation algorithm (SSA), the logarithmic direct method (LDM) and the next reaction method (NRM) is approximately 85 times faster than the sequential implementation of the NRM on a central processing unit (CPU). Using our software does not require any changes to the user's models, since it acts as a direct replacement of the stochastic simulation software of the SBTOOLBOX2. The software is open source under the GPL v3 and available at http://www.maths.ox.ac.uk/cmb/STOCHSIMGPU. The web site also contains supplementary information. klingbeil@maths.ox.ac.uk Supplementary data are available at Bioinformatics online.

CellMC—a multiplatform model compiler for the Cell Broadband Engine and ×86

Gillespie's stochastic simulation algorithm (SSA) is often the most tractable method to study stochastic models of biochemical systems. The algorithm itself is very simple and a natural target for implementation on specialized architectures such as the Cell Broadband Engine (Cell/BE). We have developed CellMC, a multiplatform SBML model compiler implementing a vectorized version of SSA for use on Cell/BE or x86 PCs. The code is freely available from http://www.cellmc.org. It will run on a wide variety of x86 computers running Linux/MacOSX (Darwin) and on Cell/BE computers such as the Sony PlayStation3 (PS3) and the IBM BladeCenter QS22. CellMC requires gcc, libxml2 and libxslt, all of which are installed by default on most of the supported platforms. Supplementary data are available at Bioinformatics online.

Chapter 5 Discrete Stochastic Simulation Methods for Chemically Reacting Systems

Discrete stochastic chemical kinetics describe the time evolution of a chemically reacting system by taking into account the fact that, in reality, chemical species are present with integer populations and exhibit some degree of randomness in their dynamical behavior. In recent years, with the development of new techniques to study biochemistry dynamics in a single cell, there are increasing studies using this approach to chemical kinetics in cellular systems, where the small copy number of some reactant species in the cell may lead to deviations from the predictions of the deterministic differential equations of classical chemical kinetics. This chapter reviews the fundamental theory related to stochastic chemical kinetics and several simulation methods based on that theory. We focus on nonstiff biochemical systems and the two most important discrete stochastic simulation methods: Gillespie's stochastic simulation algorithm (SSA) and the tau-leaping method. Different implementation strategies of these two methods are discussed. Then we recommend a relatively simple and efficient strategy that combines the strengths of the two methods: the hybrid SSA/tau-leaping method. The implementation details of the hybrid strategy are given here and a related software package is introduced. Finally, the hybrid method is applied to simple biochemical systems as a demonstration of its application.

Detailed comparison between StochSim and SSA

Morton-Firth and Bray's stochastic simulator (StochSim) and Gillespie's stochastic simulation algorithm (SSA) are two important methods for stochastic modelling and simulation of biochemical systems. They have been widely applied to different biological problems. A key question is discussed here: Are these two methods equivalent? These two methods are compared using fundamental probability analysis. The analysis clearly shows that, when the time step in the StochSim is chosen very small, the StochSim can be viewed as a first-order approximation to the SSA. The analysis also explains why the SSA is usually much more efficient than the StochSim for biochemical systems. However, when multistate species present in a system, the StochSim clearly shows its advantage. The Complexity analysis is used to explain this advantage. The hybrid SSA (HSSA) is proposed to combine the advantages of both the StochSim and SSA. When the populations of the multistate species are small, the HSSA is very efficient. Numerical experiments are presented to verify the analysis.

A Dominated Coupling From The Past algorithm for the stochastic simulation of networks of biochemical reactions

BackgroundIn recent years, stochastic descriptions of biochemical reactions based on the Master Equation (ME) have become widespread. These are especially relevant for models involving gene regulation. Gillespie’s Stochastic Simulation Algorithm (SSA) is the most widely used method for the numerical evaluation of these models. The SSA produces exact samples from the distribution of the ME for finite times. However, if the stationary distribution is of interest, the SSA provides no information about convergence or how long the algorithm needs to be run to sample from the stationary distribution with given accuracy. ResultsWe present a proof and numerical characterization of a Perfect Sampling algorithm for the ME of networks of biochemical reactions prevalent in gene regulation and enzymatic catalysis. Our algorithm combines the SSA with Dominated Coupling From The Past (DCFTP) techniques to provide guaranteed sampling from the stationary distribution. The resulting DCFTP-SSA is applicable to networks of reactions with uni-molecular stoichiometries and sub-linear, (anti-) monotone propensity functions. We showcase its applicability studying steady-state properties of stochastic regulatory networks of relevance in synthetic and systems biology.ConclusionThe DCFTP-SSA provides an extension to Gillespie’s SSA with guaranteed sampling from the stationary solution of the ME for a broad class of stochastic biochemical networks.

Uncertainty Quantification for Atomistic Reaction Models: An Equation-Free Stochastic Simulation Algorithm Example

We describe a computational framework linking uncertainty quantification (UQ) methods for continuum problems depending on random parameters with equation-free (EF) methods for performing continuum deterministic numerics by acting directly on atomistic/stochastic simulators. Our illustrative example is a heterogeneous catalytic reaction mechanism with an uncertain atomistic kinetic parameter; the “inner” dynamic simulator of choice is a Gillespie stochastic simulation algorithm (SSA). We demonstrate UQ computations at the coarse-grained level, in a nonintrusive way, through the design of brief, appropriately initialized computational experiments with the SSA code. The system is thus observed at three levels: (a) a fine scale for each stochastic simulation at each value of the uncertain parameter; (b) an intermediate coarse-grained state for the expected behavior of the SSA at each value of the uncertain parameter; and (c) the desired fully coarse-grained level: distributions of the coarse-grained behavior ...

GillespieSSA: Implementing the Stochastic Simulation Algorithm inR

The deterministic dynamics of populations in continuous time are traditionally described using coupled, first-order ordinary differential equations. While this approach is accurate for large systems, it is often inadequate for small systems where key species may be present in small numbers or where key reactions occur at a low rate. The Gillespie stochastic simulation algorithm (SSA) is a procedure for generating time-evolution trajectories of finite populations in continuous time and has become the standard algorithm for these types of stochastic models. This article presents a simple-to-use and flexible framework for implementing the SSA using the high-level statistical computing language R and the package GillespieSSA. Using three ecological models as examples (logistic growth, Rosenzweig-MacArthur predator-prey model, and Kermack-McKendrick SIRS metapopulation model), this paper shows how a deterministic model can be formulated as a finite-population stochastic model within the framework of SSA theory and how it can be implemented in R. Simulations of the stochastic models are performed using four different SSA Monte Carlo methods: one exact method (Gillespie's direct method); and three approximate methods (explicit, binomial, and optimized tau-leap methods). Comparison of simulation results confirms that while the time-evolution trajectories obtained from the different SSA methods are indistinguishable, the approximate methods are up to four orders of magnitude faster than the exact methods.

Multiscale Modeling of Chemical Kinetics via the Master Equation

We present numerical methods for both the direct solution and simulation of the chemical master equation (CME), and, compared to popular methods in current use, such as the Gillespie stochastic simulation algorithm (SSA) and τ-Leap approximations, this new approach has the advantage of being able to detect when the system has settled down to equilibrium. This improved performance is due to the incorporation of information from the associated CME, a valuable complementary approach to the SSA that has often been felt to be too computationally inefficient. Hybrid methods, that combine these complementary approaches and so are able to detect equilibrium while maintaining the efficiency of the leap methods, are also presented. Amongst CME-solvers the recently suggested finite state projection algorithm is especially well suited to this purpose and has been adapted here for the task, leading to a type of “exact τ-Leap.” It is also observed that a CMEsolver is often more efficient than an SSA or even a τ-Leap approach for computing moments of the solution such as the mean and variance. These techniques are demonstrated on a test suite of five biologically inspired models, namely, stochastic models of the genetic toggle, receptor oligomerization, the Schl¨ogl reactions, Goutsias’ model of regulated gene transcription, and a decaying-dimerizing reaction set. For the gene toggle it is observed that important experimentally measurable traits such as the percentage of cells that undergo so-called switching may also be more efficiently approximated via a CME-based approach.

Antagonism and bistability in protein interaction networks

A protein interaction network (PIN) is a set of proteins that modulate one another's activities by regulated synthesis and degradation, by reversible binding to form complexes, and by catalytic reactions (e.g., phosphorylation and dephosphorylation). Most PINs are so complex that their dynamical characteristics cannot be deduced accurately by intuitive reasoning alone. To predict the properties of such networks, many research groups have turned to mathematical models (differential equations based on standard biochemical rate laws, e.g., mass-action, Michaelis–Menten, Hill). When using Michaelis–Menten rate expressions to model PINs, care must be exercised to avoid making inconsistent assumptions about enzyme–substrate complexes. We show that an appealingly simple model of a PIN that functions as a bistable switch is compromised by neglecting enzyme–substrate intermediates. When the neglected intermediates are put back into the model, bistability of the switch is lost. The theory of chemical reaction networks predicts that bistability can be recovered by adding specific reaction channels to the molecular mechanism. We explore two very different routes to recover bistability. In both cases, we show how to convert the original ‘phenomenological’ model into a consistent set of mass-action rate laws that retains the desired bistability properties. Once an equivalent model is formulated in terms of elementary chemical reactions, it can be simulated accurately either by deterministic differential equations or by Gillespie's stochastic simulation algorithm.