Gillespie's Stochastic Simulation Research Articles

A hybrid tau-leap for simulating chemical kinetics with applications to parameter estimation.

We consider the problem of efficiently simulating stochastic models of chemical kinetics. The Gillespie stochastic simulation algorithm (SSA) is often used to simulate these models; however, in many scenarios of interest, the computational cost quickly becomes prohibitive. This is further exacerbated in the Bayesian inference context when estimating parameters of chemical models, as the intractability of the likelihood requires multiple simulations of the underlying system. To deal with issues of computational complexity in this paper, we propose a novel hybrid τ-leap algorithm for simulating well-mixed chemical systems. In particular, the algorithm uses τ-leap when appropriate (high population densities), and SSA when necessary (low population densities, when discrete effects become non-negligible). In the intermediate regime, a combination of the two methods, which uses the properties of the underlying Poisson formulation, is employed. As illustrated through a number of numerical experiments, the hybrid τ offers significant computational savings when compared with SSA without, however, sacrificing the overall accuracy. This feature is particularly welcomed in the Bayesian inference context, as it allows for parameter estimation of stochastic chemical kinetics at reduced computational cost.

Highly accelerated kinetic Monte Carlo models for depolymerisation systems

Kinetic Monte Carlo (kMC) models are a well-established modelling framework for the simulation of complex free-radical kinetic systems. kMC models offer the advantage of discretely monitoring every chain sequence in the system, providing full accounting of the chain molecular weight distribution. These models are marred by the necessity to simulate a minimum number of molecules, which confers significant computational burden. This paper adapts and creates a highly generalizable methodology for scaling dilute radical populations in discrete stochastic models, such as Gillespie's Stochastic Simulation Algorithm (SSA). The methodology is then applied to a kMC simulation of polystyrene (PS) pyrolysis, using a modelling framework adapted from literature. The results show that the required number of simulated molecules can be successfully reduced by up to three orders of magnitude with minimal loss of convergent behaviour, corresponding to a wall-clock simulation speed reduction of between 95.2 to 99.6 % at common pyrolysis temperatures.

A Compact and Power-Efficient Noise Generator for Stochastic Simulations.

This paper describes an adaptive noise generator circuit suitable for on-chip simulations of stochastic chemical kinetics. The circuit uses amplified BJT white noise and adaptive low-pass filtering to emulate the power spectrum and autocorrelation of random telegraph signals (RTS) with Poisson-distributed level transitions. A current-mode implementation in the AMS 0.35 μm BiCMOS process shows excellent agreement with theoretical results from the Gillespie stochastic simulation algorithm over a 60 dB range in mean current levels (modeling molecule count numbers). The circuit has an estimated layout area of 0.032 mm2 and typically consumes 400 μA, which are 73% and 50% less, respectively, than prior implementations. Moreover, it does not require any off-chip capacitors. Experimental results from a discrete board-level implementation of the circuit are in good agreement with theoretical predictions.

TiPS: Rapidly simulating trajectories and phylogenies from compartmental models

Abstract Stochastic population dynamics simulations are essential for many ecological and epidemiological studies to generate time series and genealogies that capture the relatedness between individuals. Many software packages allow one to simulate phylogenetic trees but these tend to suffer from one or two major limitations. First, the underlying population dynamics model is often simplistic (e.g. constant population size or exponential growth). Second, the software packages are not appropriate to simulate a large number of trees. We introduce TiPS, an R package to generate trajectories and phylogenetic trees associated with a compartmental model. Trajectories are simulated using Gillespie's exact or approximate stochastic simulation algorithm, or a newly proposed mixed version of the two. Phylogenetic trees are simulated from a trajectory under a backwards‐in‐time approach (i.e. coalescent). TiPS is based on the Rcpp package, allowing to combine the flexibility of R for model definition and the speed of C++ for simulations execution. The model is defined in R with a set of reactions, which allow capturing heterogeneity in life cycles or any sort of population structure. TiPS converts the model into C++ code and compiles it into a simulator that is interfaced in R via a function. TiPS is flexible, easy‐to‐use and available on the CRAN at https://cran.r‐project.org/package=TiPS. Plus, benchmarking analyses show that TiPS is faster than existing packages. This package is particularly well suited for population genetics and phylodynamics studies that need to generate a large number of phylogenies used for population dynamics studies.

Insights into the Active Catalyst Formation from Dinuclear Palladium Acetate in Pd-Catalyzed Coupling Reactions: A DFT Study.

We studied the steps in the formation of active palladium catalyst species from palladium acetate dimer employing density functional theory calculations. We explored the possible pathways with an automated reaction search and studied the kinetics with stochastic simulation analysis. The dimeric form of palladium acetate is considered a resting state of the catalyst. Our reaction search starting from the dimeric form by sequential ligand addition resulted in experimentally observed monomeric species. We analyzed the bonding in the palladium acetate dimer and the role of Pd in the stability of the dimeric species. We implemented the Gillespie stochastic simulation algorithm to gain more insights into multichannel reaction paths and applied it to the degradation pathways. The analysis of the thermodynamic and kinetic data for these degradation pathways suggests that the dimeric form of the catalyst can be a potential catalytic precursor in the palladium acetate-catalyzed coupling reactions under the experimental reaction conditions.

Reversible bond kinetics from single-molecule force spectroscopy experiments close to equilibrium

Analysis of bond rupture data from single-molecule force spectroscopy experiments commonly relies on the strong assumption that the bond dissociation process is irreversible. However, with increased spatiotemporal resolution of instruments it is now possible to observe multiple unbinding-rebinding (or unfolding-refolding) events in a single pulling experiment. Here we augment the theory of force-induced unbinding by explicitly taking into account rebinding kinetics and provide approximate analytic solutions of the resulting rate equations. Furthermore, we use a short-time expansion of the exact kinetics to construct numerically efficient maximum likelihood estimators for the parameters of the force-dependent unbinding and rebinding rates, which pair well with and complement established methods, such as the analysis of rate maps. The estimators are independent of the applied force protocol and can therefore be used to globally analyze multiple force-extension traces for both force-ramp and force-clamp experiments. We provide an open-source implementation of the theory, evaluated for Bell-like rates, which we apply to synthetic data generated by a Gillespie stochastic simulation algorithm for time-dependent rates.

TemporalGSSA: A numerically robust R-wrapper to facilitate computation of a metabolite-specific and simulation time-dependent trajectory from stochastic simulation algorithm (SSA)-generated datasets

Whilst data on biochemical networks has increased several-fold, our comprehension of the underlying molecular biology is incomplete and inadequate. Simulation studies permit data collation from disparate time points and the imputed trajectories can provide valuable insights into the molecular biology of complex biochemical systems. Although, stochastic simulations are accurate, each run is an independent event and the data that is generated cannot be directly compared even with identical simulation times. This lack of robustness will preclude a biologically meaningful result for the metabolite(s) of concern and is a significant limitation of this approach. "TemporalGSSA" or temporal Gillespie Stochastic Simulation Algorithm is an R-wrapper which will collate and partition SSA-generated datasets with identical simulation times (trials) into finite sets of linear models (technical replicates). Each such model (time step of a single run, absolute number of molecules for a metabolite) computes several coefficients (slope, intercept, etc.). These coefficients are averaged (mean slope, mean intercept) across all trials of a technical replicate and along with an imputed time step (mean, median, random) is incorporated into a linear regression equation. The solution to this equation is the number of molecules of a metabolite which is used to compute the molar concentration of the metabolite per technical replicate. The summarized (mean, standard deviation) data of this vector of technical replicates is the outcome or numerical estimate of the molar concentration of a metabolite and is dependent on the duration of the simulation. If the SSA-generated dataset comprises runs with differing simulation times, "TemporalGSSA" can compute the time-dependent trajectory of a metabolite provided the trials-per technical replicate constraint is complied with. The algorithms deployed by "TemporalGSSA" are rigorous, have a sound theoretical basis and have contributed meaningfully to our comprehension of the mechanism(s) that drive complex biochemical systems. "TemporalGSSA", is robust, freely accessible and easy to use with several readily testable examples.

TauSSA

In this paper, we present TauSSA, a discrete-event simulation tool for stochastic queueing networks integrated in the LINE solver. TauSSA combines Gillespie's stochastic simulation algorithm with tau leaping, a methodology for optimistic simulation acceleration. Although tau leaping is frequently used in chemical reaction network simulation, it has so far found limited application in queueing theory. TauSSA offers one of the very first attempts to make this method broadly applicable to analyze extended queueing network models, which include class switching, fork-join, and non-exponential service and arrival distributions. We conceptualize various strategies for handling ordering and illegal states in tau leaping that arise specifically within queueing network models, and compare their performance through numerical experiments. Our main finding is that strategies that sort events based on the network topological order incur a better trade-off between speedup and approximation error.

Multiple epidemic waves as the outcome of stochastic SIR epidemics with behavioral responses: a hybrid modeling approach.

In the behavioral epidemiology (BE) of infectious diseases, little theoretical effort seems to have been devoted to understand the possible effects of individuals’ behavioral responses during an epidemic outbreak in small populations. To fill this gap, here we first build general, behavior implicit, SIR epidemic models including behavioral responses and set them within the framework of nonlinear feedback control theory. Second, we provide a thorough investigation of the effects of different types of agents’ behavioral responses for the dynamics of hybrid stochastic SIR outbreak models. In the proposed model, the stochastic discrete dynamics of infection spread is combined with a continuous model describing the agents’ delayed behavioral response. The delay reflects the memory mechanisms with which individuals enact protective behavior based on past data on the epidemic course. This results in a stochastic hybrid system with time-varying transition probabilities. To simulate such system, we extend Gillespie’s classic stochastic simulation algorithm by developing analytical formulas valid for our classes of models. The algorithm is used to simulate a number of stochastic behavioral models and to classify the effects of different types of agents’ behavioral responses. In particular this work focuses on the effects of the structure of the response function and of the form of the temporal distribution of such response. Among the various results, we stress the appearance of multiple, stochastic epidemic waves triggered by the delayed behavioral response of individuals.

Transition from growth to decay of an epidemic due to lockdown

We study the transition of an epidemic from growth phase to decay of the active infections in a population when lockdown health measures are introduced to reduce the probability of disease transmission. Although in the case of uniform lockdown, a simple compartmental model would indicate instantaneous transition to decay of the epidemic, this is not the case when partially isolated active clusters remain with the potential to create a series of small outbreaks. We model this using the Gillespie stochastic simulation algorithm based on a connected set of stochastic susceptible-infected-removed/recovered networks representing the locked-down majority population (in which the reproduction number is less than 1) weakly coupled to a large set of small clusters in which the infection may propagate. We find that the presence of such active clusters can lead to slower than expected decay of the epidemic and significantly delayed onset of the decay phase. We study the relative contributions of these changes, caused by the active clusters within the population, to the additional total infected population. We also demonstrate that limiting the size of the inevitable active clusters can be efficient in reducing their impact on the overall size of the epidemic outbreak. The deceleration of the decay phase becomes apparent when the active clusters form at least 5% of the population.

Heterogeneous responses to low level death receptor activation are explained by random molecular assembly of the Caspase-8 activation platform.

Ligand binding to death receptors activates apoptosis in cancer cells. Stimulation of death receptors results in the formation of intracellular multiprotein platforms that either activate the apoptotic initiator Caspase-8 to trigger cell death, or signal through kinases to initiate inflammatory and cell survival signalling. Two of these platforms, the Death-Inducing Signalling Complex (DISC) and the RIPoptosome, also initiate necroptosis by building filamentous scaffolds that lead to the activation of mixed lineage kinase domain-like pseudokinase. To explain cell decision making downstream of death receptor activation, we developed a semi-stochastic model of DISC/RIPoptosome formation. The model is a hybrid of a direct Gillespie stochastic simulation algorithm for slow assembly of the RIPoptosome and a deterministic model of downstream caspase activation. The model explains how alterations in the level of death receptor-ligand complexes, their clustering properties and intrinsic molecular fluctuations in RIPoptosome assembly drive heterogeneous dynamics of Caspase-8 activation. The model highlights how kinetic proofreading leads to heterogeneous cell responses and results in fractional cell killing at low levels of receptor stimulation. It reveals that the noise in Caspase-8 activation—exclusively caused by the stochastic molecular assembly of the DISC/RIPoptosome platform—has a key function in extrinsic apoptotic stimuli recognition.

DIFFpop: a stochastic computational approach to simulate differentiation hierarchies with single cell barcoding.

DIFFpop is an R package designed to simulate cellular differentiation hierarchies using either exponentially-expanding or fixed population sizes. The software includes functionalities to simulate clonal evolution due to the emergence of driver mutations under the infinite-allele assumption as well as options for simulation and analysis of single cell barcoding and labeling data. The software uses the Gillespie Stochastic Simulation Algorithm and a modification of expanding or fixed-size stochastic process models expanded to a large number of cell types and scenarios. DIFFpop is available as an R-package along with vignettes on Github (https://github.com/ferlicjl/diffpop). Supplementary data are available at Bioinformatics online.

SimInf: An R Package for Data-Driven Stochastic Disease Spread Simulations

We present the R package SimInf which provides an efficient and very flexible framework to conduct data-driven epidemiological modeling in realistic large scale disease spread simulations. The framework integrates infection dynamics in subpopulations as continuous-time Markov chains using the Gillespie stochastic simulation algorithm and incorporates available data such as births, deaths and movements as scheduled events at predefined time-points. Using C code for the numerical solvers and OpenMP to divide work over multiple processors ensures high performance when simulating a sample outcome. One of our design goal was to make SimInf extendable and enable usage of the numerical solvers from other R extension packages in order to facilitate complex epidemiological research. In this paper, we provide a technical description of the framework and demonstrate its use on some basic examples. We also discuss how to specify and extend the framework with user-defined models.

S-Leaping: An Adaptive, Accelerated Stochastic Simulation Algorithm, Bridging [Formula: see text]-Leaping and R-Leaping.

We propose the S-leaping algorithm for the acceleration of Gillespie's stochastic simulation algorithm that combines the advantages of the two main accelerated methods; the [Formula: see text]-leaping and R-leaping algorithms. These algorithms are known to be efficient under different conditions; the [Formula: see text]-leaping is efficient for non-stiff systems or systems with partial equilibrium, while the R-leaping performs better in stiff system thanks to an efficient sampling procedure. However, even a small change in a system's set up can critically affect the nature of the simulated system and thus reduce the efficiency of an accelerated algorithm. The proposed algorithm combines the efficient time step selection from the [Formula: see text]-leaping with the effective sampling procedure from the R-leaping algorithm. The S-leaping is shown to maintain its efficiency under different conditions and in the case of large and stiff systems or systems with fast dynamics, the S-leaping outperforms both methods. We demonstrate the performance and the accuracy of the S-leaping in comparison with the [Formula: see text]-leaping and R-leaping on a number of benchmark systems involving biological reaction networks.

Accuracy Analysis of Hybrid Stochastic Simulation Algorithm on Linear Chain Reaction Systems.

Noise in cellular systems is often modeled and simulated with Gillespie's stochastic simulation algorithm (SSA), but the low efficiency of the SSA limits its application to large biochemical networks. To improve the efficiency of stochastic simulations, Haseltine and Rawlings (HR) proposed a hybrid algorithm, which combines ordinary differential equations for traditional deterministic models and the SSA for stochastic models. In this paper, accuracy of the HR hybrid method is studied based on a linear chain reaction system. Mathematical analysis and numerical results both show that the HR hybrid method is accurate if either the quantity of reactant molecules in fast reactions is above a certain threshold, or the reaction rates of fast reactions are much larger than those of slow reactions. This analysis also shows that the HR hybrid method approximates the chemical master equation well for a much greater region in system parameter space than the slow-scale SSA and the stochastic quasi-steady-state assumption methods.

Selected-node stochastic simulation algorithm.

Stochastic simulations of biochemical networks are of vital importance for understanding complex dynamics in cells and tissues. However, existing methods to perform such simulations are associated with computational difficulties and addressing those remains a daunting challenge to the present. Here we introduce the selected-node stochastic simulation algorithm (snSSA), which allows us to exclusively simulate an arbitrary, selected subset of molecular species of a possibly large and complex reaction network. The algorithm is based on an analytical elimination of chemical species, thereby avoiding explicit simulation of the associated chemical events. These species are instead described continuously in terms of statistical moments derived from a stochastic filtering equation, resulting in a substantial speedup when compared to Gillespie's stochastic simulation algorithm (SSA). Moreover, we show that statistics obtained via snSSA profit from a variance reduction, which can significantly lower the number of Monte Carlo samples needed to achieve a certain performance. We demonstrate the algorithm using several biological case studies for which the simulation time could be reduced by orders of magnitude.

Probabilistic models of individual and collective animal behavior

Recent developments in automated tracking allow uninterrupted, high-resolution recording of animal trajectories, sometimes coupled with the identification of stereotyped changes of body pose or other behaviors of interest. Analysis and interpretation of such data represents a challenge: the timing of animal behaviors may be stochastic and modulated by kinematic variables, by the interaction with the environment or with the conspecifics within the animal group, and dependent on internal cognitive or behavioral state of the individual. Existing models for collective motion typically fail to incorporate the discrete, stochastic, and internal-state-dependent aspects of behavior, while models focusing on individual animal behavior typically ignore the spatial aspects of the problem. Here we propose a probabilistic modeling framework to address this gap. Each animal can switch stochastically between different behavioral states, with each state resulting in a possibly different law of motion through space. Switching rates for behavioral transitions can depend in a very general way, which we seek to identify from data, on the effects of the environment as well as the interaction between the animals. We represent the switching dynamics as a Generalized Linear Model and show that: (i) forward simulation of multiple interacting animals is possible using a variant of the Gillespie’s Stochastic Simulation Algorithm; (ii) formulated properly, the maximum likelihood inference of switching rate functions is tractably solvable by gradient descent; (iii) model selection can be used to identify factors that modulate behavioral state switching and to appropriately adjust model complexity to data. To illustrate our framework, we apply it to two synthetic models of animal motion and to real zebrafish tracking data.

The mean and noise of stochastic gene transcription with cell division.

Life growth and development are driven by continuous cell divisions. Cell division is a stochastic and complex process. In this paper, we study the impact of cell division on the mean and noise of mRNA numbers by using a two-state stochastic model of transcription. Our results show that the steady-state mRNA noise with symmetric cell division is less than that with binomial inheritance with probability 0.5, but the steady-state mean transcript level with symmetric division is always equal to that with binomial inheritance with probability 0.5. Cell division except random additive inheritance always decreases mean transcript level and increases transcription noise. Inversely, random additive inheritance always increases mean transcript level and decreases transcription noise. We also show that the steady-state mean transcript level (the steady-state mRNA noise) with symmetric cell division or binomial inheritance increases (decreases) with the average cell cycle duration. But the steady-state mean transcript level (the steady-state mRNA noise) with random additive inheritance decreases (increases) with the average cell cycle duration. Our results are confirmed by Gillespie stochastic simulation using plausible parameters.

Noise slows the rate of Michaelis–Menten reactions

Microscopic randomness and the small volumes of living cells combine to generate random fluctuations in molecule concentrations called “noise”. Here I investigate the effect of noise on biochemical reactions obeying Michaelis–Menten kinetics, concluding that substrate noise causes these reactions to slow. I derive a general expression for the time evolution of the joint probability density of chemical species in arbitrarily connected networks of non-linear chemical reactions in small volumes. This equation is a generalization of the chemical master equation (CME), a common tool for investigating stochastic chemical kinetics, extended to reaction networks occurring in small volumes, such as living cells. I apply this equation to a generalized Michaelis-Menten reaction in an open system, deriving the following general result: 〈p〉≤p¯ and 〈s〉≥s¯, where s¯ and p¯ denote the deterministic steady-state concentration of reactant and product species, respectively, and 〈s〉 and 〈p〉 denote the steady-state ensemble average over independent realizations of a stochastic reaction. Under biologically realistic conditions, namely when substrate is degraded or diluted by cell division, 〈p〉≤p¯. Consequently, noise slows the rate of in vivo Michaelis–Menten reactions. These predictions are validated by extensive stochastic simulations using Gillespie's exact stochastic simulation algorithm. I specify the conditions under which these effects occur and when they vanish, therefore reconciling discrepancies among previous theoretical investigations of stochastic biochemical reactions. Stochastic slowdown of reaction flux caused by molecular noise in living cells may have functional consequences, which the present theory may be used to quantify.

The impact of fluctuations and correlations in droplet growth by collision–coalescence revisited – Part 1: Numerical calculation of post-gel droplet size distribution

Abstract. The impact of stochastic fluctuations in cloud droplet growth is a matter of broad interest, since stochastic effects are one of the possible explanations of how cloud droplets cross the size gap and form the raindrop embryos that trigger warm rain development in cumulus clouds. Most theoretical studies on this topic rely on the use of the kinetic collection equation, or the Gillespie stochastic simulation algorithm. However, the kinetic collection equation is a deterministic equation with no stochastic fluctuations. Moreover, the traditional calculations using the kinetic collection equation are not valid when the system undergoes a transition from a continuous distribution to a distribution plus a runaway raindrop embryo (known as the sol–gel transition). On the other hand, the stochastic simulation algorithm, although intrinsically stochastic, fails to adequately reproduce the large end of the droplet size distribution due to the huge number of realizations required. Therefore, the full stochastic description of cloud droplet growth must be obtained from the solution of the master equation for stochastic coalescence. In this study the master equation is used to calculate the evolution of the droplet size distribution after the sol–gel transition. These calculations show that after the formation of the raindrop embryo, the expected droplet mass distribution strongly differs from the results obtained with the kinetic collection equation. Furthermore, the low-mass bins and bins from the gel fraction are strongly anticorrelated in the vicinity of the critical time, this being one of the possible explanations for the differences between the kinetic and stochastic approaches after the sol–gel transition. Calculations performed within the stochastic framework provide insight into the inability of explicit microphysics cloud models to explain the droplet spectral broadening observed in small, warm clouds.