Exact Stochastic Simulation Algorithm Research Articles

Stochastic Simulations of Casual Groups

Free-forming or casual groups are groups in which individuals are in face-to-face interactions and are free to maintain or terminate contact with one another, such as clusters of people at a cocktail party, play groups in a children’s playground or shopping groups in a mall. Stochastic models of casual groups assume that group sizes are the products of natural processes by which groups acquire and lose members. The size distributions predicted by these models have been the object of controversy since their derivation in the 1960s because of the neglect of fluctuations around the mean values of random variables that characterize a collection of groups. Here, we check the validity of these mean-field approximations using an exact stochastic simulation algorithm to study the processes of the acquisition and loss of group members. In addition, we consider the situation where the appeal of a group of size i to isolates is proportional to iα. We find that, for α≤1, the mean-field approximation fits the equilibrium simulation results very well, even for a relatively small population size N. However, for α>1, this approximation scheme fails to provide a coherent description of the distribution of group sizes. We find a discontinuous phase transition at αc>1 that separates the regime where the variance of the group size does not depend on N from the regime where it grows linearly with N. In the latter regime, the system is composed of a single large group that coexists with a large number of isolates. Hence, the same underlying acquisition-and-loss process can explain the existence of small, temporary casual groups and of large, stable social groups.

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.

A Gaussian jump process formulation of the reaction-diffusion master equation enables faster exact stochastic simulations.

We propose a Gaussian jump process model on a regular Cartesian lattice for the diffusion part of the Reaction-Diffusion Master Equation (RDME). We derive the resulting Gaussian RDME (GRDME) formulation from analogy with a kernel-based discretization scheme for continuous diffusion processes and quantify the limits of its validity relative to the classic RDME. We then present an exact stochastic simulation algorithm for the GRDME, showing that the accuracies of GRDME and RDME are comparable, but exact simulations of the GRDME require only a fraction of the computational cost of exact RDME simulations. We analyze the origin of this speedup and its scaling with problem dimension. The benchmarks suggest that the GRDME is a particularly beneficial model for diffusion-dominated systems in three dimensional spaces, often occurring in systems biology and cell biology.

Intrinsic stochastic resonance via set-point variation.

In the present paper, the possibility of invoking stochastic resonance (SR, periodic and aperiodic) by regulating the operating value of an appropriate parameter is explored. The operating values of these parameters are defined as the set point of the system throughout the present paper. Brusselator, a mathematical model [I. Prigogine and R. Lefever, J. Chem. Phys. 48, 1695 (1968)JCPSA60021-960610.1063/1.1668896] of nonlinear chemical reactions, is used for this purpose. We consider the effect of intrinsic noise in the Brusselator due to the Markovian nature of the chemical reactions. The stochastic time evolution is studied using the Gillespie algorithm [D. T. Gillespie, J. Comput. Phys. 22, 403 (1976)JCTPAH0021-999110.1016/0021-9991(76)90041-3], which is an exact stochastic simulation algorithm. We analyze the dependence of the resonance point on both the strength of the intrinsic noise as well as the distance from the bifurcation point. Subsequently, the phenomena of SR is explored using both periodic and aperiodic stimulus. It was found that, for a given system size, in both cases, SR is achieved by variation of the set point. Set-point variation can be achieved by regulating either the source concentration or the rate constants. Resonance is observed in both cases. However, this resonance occurs at different values of the set point, even with a fixed system size. This is clearly seen in the set-point versus system-size plane, where the resonance line has different slopes for the two scenarios. Our semianalytic treatment points to the fact that for a given system size intrinsic noise is affected differently for different methods involving the variation of the set point. This is explained by writing the corresponding chemical Langevin equation and comparing the various intrinsic noise sources.

Slow update stochastic simulation algorithms for modeling complex biochemical networks

The stochastic simulation algorithm (SSA) based modeling is a well recognized approach to predict the stochastic behavior of biological networks. The stochastic simulation of large complex biochemical networks is a challenge as it takes a large amount of time for simulation due to high update cost. In order to reduce the propensity update cost, we proposed two algorithms: slow update exact stochastic simulation algorithm (SUESSA) and slow update exact sorting stochastic simulation algorithm (SUESSSA). We applied cache-based linear search (CBLS) in these two algorithms for improving the search operation for finding reactions to be executed. Data structure used for incorporating CBLS is very simple and the cost of maintaining this during propensity update operation is very low. Hence, time taken during propensity updates, for simulating strongly coupled networks, is very fast; which leads to reduction of total simulation time. SUESSA and SUESSSA are not only restricted to elementary reactions, they support higher order reactions too.We used linear chain model and colloidal aggregation model to perform a comparative analysis of the performances of our methods with the existing algorithms. We also compared the performances of our methods with the existing ones, for large biochemical networks including B cell receptor and FcϵRI signaling networks.

Stochastic simulation of biochemical reactions with partial-propensity and rejection-based approaches

We present in this paper a new exact algorithm for improving performance of exact stochastic simulation algorithm. The algorithm is developed on concepts of the partial-propensity and the rejection-based approaches. It factorizes the propensity bounds of reactions and groups factors by common reactant species for selecting next reaction firings. Our algorithm provides favorable computational advantages for simulating of biochemical reaction networks by reducing the cost for selecting the next reaction firing to scale with the number of chemical species and avoiding expensive propensity updates during the simulation. We present the details of our new algorithm and benchmark it on concrete biological models to demonstrate its applicability and efficiency.

APPROXIMATE BAYESIAN PARAMETER INFERENCE FOR DYNAMICAL SYSTEMS IN SYSTEMS BIOLOGY

A b s t r a c t: This paper proposes to use approximate instead of exact stochastic simulation algorithms for approximate Bayesian parameter inference of dynamical systems in systems biology. It first presents the mathematical framework for the description of systems biology models, especially from the aspect of a stochastic formulation as opposed to deterministic model formulations based on the law of mass action. In contrast to maximum likelihood methods for parameter inference, approximate inference methodsare presented which are based on sampling parameters from a known prior probability distribution, which gradually evolves tward a posterior distribution, through the comparison of simulated data from the model to a given data set of measurements. The paper then discusses the simulation process, where an overview is given of the different exact and approximate methods for stochastic simulation and their improvements that we propose. The exact and approximate simulators are implemented and used within approximate Bayesian parameter inference methods. Our evaluation of these methods on two tasks of parameter estimation in two different models shows that equally good results are obtained much faster when using approximate simulation as compared to using exact simulation.

Using Conserved Cycles in Exact Stochastic Simulation Algorithms

Biyokimyasal reaksiyon sistemleri farkli reaksiyonlar araciligiyla etkilesime giren bircok farkli turu icerir. Sistem icerisinde yer alan turlerin sayilari ve miktarlari cok yuksek oldugunda, diferansiyel denklemlere dayanan saf modelleme yaklasimlari cok boyutluluktan muzdarip olurlar. Eger bir sistem korunumlu donguler icerirse, bazi turlerin miktarlari cebirsel baglantilar yoluyla elde edilebilir bu da sistemin dinamiklerini temsil eden diferansiyel denklemlerin boyutunu dusurur. Bu calismada, biyokimyasal reaksiyon sistemlerinde yer alan korunumlu donguleri elde etmek icin Gauss-Jordan metodunu kullanan bir numerik algoritma oneriyoruz. Algoritmayi stokastik modelleme yaklasiminda konum vektorunun tam realizasyonlarini elde eden Direk Metod (DM), Ilk Reaksiyon Metodu (FRM) ve Sonraki Reaksiyon Metodu (FRM) icerisinde verdik. Bu uc algoritmayi korunum bagintilarini icerecek/icermecek sekilde farkli boyutlardaki biyokimyasal sistemlere uyguladik ve her tam lagoritmanin farkli iki versiyonunun hesaplama miktarlari kiyasladik.

Multi-level methods and approximating distribution functions

Biochemical reaction networks are often modelled using discrete-state, continuous-time Markov chains. System statistics of these Markov chains usually cannot be calculated analytically and therefore estimates must be generated via simulation techniques. There is a well documented class of simulation techniques known as exact stochastic simulation algorithms, an example of which is Gillespie’s direct method. These algorithms often come with high computational costs, therefore approximate stochastic simulation algorithms such as the tau-leap method are used. However, in order to minimise the bias in the estimates generated using them, a relatively small value of tau is needed, rendering the computational costs comparable to Gillespie’s direct method. The multi-level Monte Carlo method (Anderson and Higham, Multiscale Model. Simul. 10:146–179, 2012) provides a reduction in computational costs whilst minimising or even eliminating the bias in the estimates of system statistics. This is achieved by first crudely approximating required statistics with many sample paths of low accuracy. Then correction terms are added until a required level of accuracy is reached. Recent literature has primarily focussed on implementing the multi-level method efficiently to estimate a single system statistic. However, it is clearly also of interest to be able to approximate entire probability distributions of species counts. We present two novel methods that combine known techniques for distribution reconstruction with the multi-level method. We demonstrate the potential of our methods using a number of examples.

HRSSA – Efficient hybrid stochastic simulation for spatially homogeneous biochemical reaction networks

This paper introduces HRSSA (Hybrid Rejection-based Stochastic Simulation Algorithm), a new efficient hybrid stochastic simulation algorithm for spatially homogeneous biochemical reaction networks. HRSSA is built on top of RSSA, an exact stochastic simulation algorithm which relies on propensity bounds to select next reaction firings and to reduce the average number of reaction propensity updates needed during the simulation. HRSSA exploits the computational advantage of propensity bounds to manage time-varying transition propensities and to apply dynamic partitioning of reactions, which constitute the two most significant bottlenecks of hybrid simulation. A comprehensive set of simulation benchmarks is provided for evaluating performance and accuracy of HRSSA against other state of the art algorithms.

Comparison of Different Dynamic Monte Carlo Methods for the Simulation of Olefin Polymerization

SummaryIn this work, Monte Carlo methods were used to simulate olefin polymerization with coordination catalysts: the Direct method (DM), the First Reaction method (FRM), the Next Reaction method (NRM), and the τ‐Leaping method. The first three methods are exact stochastic simulation algorithms (SSA), while the τ‐Leaping is an approximate method with faster computation times. It is shown that all four methods predict similar polymer microstructures, but require significantly different computation times. The τ‐Leaping method is the fastest, being recommended when complex polymerization mechanisms are being investigated. The NRM, because of its intelligent data storage and handling approach, is the best among the SSA.

Fast and accurate representations of stochastic ion channel fluctuations

Random ion channel gating is an important source of noise at the single neuron level. For mesoscale ion channel population sizes, fast, accurate representation of channel noise fluctuations remains an important challenge. We highlight recent progress in three areas. In [1] we present an exact stochastic simulation algorithm, which takes into account the time dependence of the voltage sensitive transitions due to rapid voltage changes during action potentials. The exact algorithm is similar to a widely used approximate stochastic simulation algorithm, in which transition rates are held fixed during the intervals between channel state transitions. We compare the algorithms and show that they are inequivalent in a strong sense, meaning that sample paths diverge when driven with identical Poisson processes, leading to different precise firing times. But the stationary histograms produced are practically indistinguishable for modest channel numbers (circa N>100), indicating weak equivalence. For channel numbers in the hundreds or higher, numerical stochastic differential equations (SDE) algorithms based on the system size expansion can be significantly faster than simulations based on discrete Markov chain simulations, while retaining reasonable accuracy. For large networks, however, even SDE based simulations become costly, particularly as more complex channel gating schemes are introduced. Schmandt and Galan introduced a stochastic shielding approximation as a fast, accurate way of simulating stochastic ion channel kinetics [2]. In the SDE representation, each edge in the graph generates both a mean population flux between adjacent nodes, and a fluctuation about the mean. Only the fluctuations arising from edges connecting functionally distinct states directly affect fluctuations in the observed behavior of the cell. In [3] we analyze stochastic shielding both for the HH sodium and potassium gating models, and for an ensemble of random graphs. We derive a quantitative measure of edge importance related to the eigenvalue/eigenvector decomposition of the graph Laplacian matrix. Channel noise makes a regularly spiking neuron a stochastic oscillator. The classical definition for the asymptotic phase of an oscillator breaks down when stochasticity is taken into account. Alternative definitions of the ``phase'' have been based on the mean first passage time property of a system of isochronal surfaces [4] and in terms of the spectral decomposition of a the adjoint Kolmogorov operator [5,6]. In the vanishing noise limit, the spectrally derived isochrons approach those of the underlying mean field system, if the latter has a finite period limit cycle. Together, these results expand our analytic, numerical, and conceptual tools for understanding the effects of random ion channel gating in conductance based neural models.

Constant-complexity stochastic simulation algorithm with optimal binning.

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.

Microdomain calcium fluctuations as a colored noise process.

Calcium ions play a key role in subcellular signaling as localized transients of the intracellular calcium concentration modify the activity of ion channels, enzymes and transcription factors, among others. The intracellular calcium concentration is inherently noisy, as diffusion, the transient binding to and dissociation from buffer molecules and stochastically gating calcium channels contribute to the fluctuations of the local copy number of Ca2+ ions. We study the properties of the fluctuating calcium concentration in sub-femtoliter volumes using an exact stochastic simulation algorithm and approximations to the exact stochastic solution. It is shown that the time course of the local calcium concentration represents a colored noise process whose autocorrelation time is a function of buffer kinetics and diffusion constants. Using the chemical Langevin description and the excess buffer approximation of the process, fast approximative algorithms and theoretical connections to the Ornstein-Uhlenbeck process are obtained. In a generic example, we show how calcium noise can couple to the dynamics of a single variable moving in a double-well potential, leading to a colored noise induced transition. Our work shows how a multitude of intracellular signaling pathways may be influenced by the inherent stochasticity of calcium signals, a key messenger in virtually any cell type, and how the calcium signal can be implemented efficiently in cellular signaling models.

Efficient rejection-based simulation of biochemical reactions with stochastic noise and delays.

We propose a new exact stochastic rejection-based simulation algorithm for biochemical reactions and extend it to systems with delays. Our algorithm accelerates the simulation by pre-computing reaction propensity bounds to select the next reaction to perform. Exploiting such bounds, we are able to avoid recomputing propensities every time a (delayed) reaction is initiated or finished, as is typically necessary in standard approaches. Propensity updates in our approach are still performed, but only infrequently and limited for a small number of reactions, saving computation time and without sacrificing exactness. We evaluate the performance improvement of our algorithm by experimenting with concrete biological models.

The N-leap method for stochastic simulation of coupled chemical reactions

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.

Stochastic Modeling in Systems Biology

Many cellular behaviors are regulated by gene regulation networks, kinetics of which is one of the main subjects in the study of systems biology. Because of the low number molecules in these reacting systems, stochastic effects are significant. In recent years, stochasticity in modeling the kinetics of gene regulation networks have been drawing the attention of many researchers. This paper is a self contained review trying to provide an overview of stochastic modeling. I will introduce the derivation of the main equations for modeling the biochemical systems with intrinsic noise (chemical master equation, Fokker-Plan equation, reaction rate equation, chemical Langevin equation), and will discuss the relations between these formulations. The mathematical formulations for systems with fluctuations in kinetic parameters are also discussed. Finally, I will introduce the exact stochastic simulation algorithm and the approximate explicit tau-leaping method for making numerical simulations.

A new approximate method for the stochastic simulation of chemical systems: The representative reaction approach

We have developed two new approximate methods for stochastically simulating chemical systems. The methods are based on the idea of representing all the reactions in the chemical system by a single reaction, i.e., by the "representative reaction approach" (RRA). Discussed in the article are the concepts underlying the new methods along with flowchart with all the steps required for their implementation. It is shown that the two RRA methods {with the reaction 2A −> B as the representative reaction (RR)} perform creditably with regard to accuracy and computational efficiency, in comparison to the exact stochastic simulation algorithm (SSA) developed by Gillespie and are able to successfully reproduce at least the first two moments of the probability distribution of each species in the systems studied. As such, the RRA methods represent a promising new approach for stochastically simulating chemical systems.

Intrinsic noise induced resonance in presence of sub-threshold signal in Brusselator

In a system of non-linear chemical reactions called the Brusselator, we show that intrinsic noise can be regulated to drive it to exhibit resonance in the presence of a sub-threshold signal. The phenomena of periodic stochastic resonance and aperiodic stochastic resonance, hitherto studied mostly with extrinsic noise, is demonstrated here to occur with inherent systemic noise using exact stochastic simulation algorithm due to Gillespie. The role of intrinsic noise in a couple of other phenomena is also discussed.

AESS: Accelerated Exact Stochastic Simulation

The Stochastic Simulation Algorithm (SSA) developed by Gillespie provides a powerful mechanism for exploring the behavior of chemical systems with small species populations or with important noise contributions. Gene circuit simulations for systems biology commonly employ the SSA method, as do ecological applications. This algorithm tends to be computationally expensive, so researchers seek an efficient implementation of SSA. In this program package, the Accelerated Exact Stochastic Simulation Algorithm (AESS) contains optimized implementations of Gillespieʼs SSA that improve the performance of individual simulation runs or ensembles of simulations used for sweeping parameters or to provide statistically significant results. Program summary Program title: AESS Catalogue identifier: AEJW_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEJW_v1_0.html Program obtainable from: CPC Program Library, Queenʼs University, Belfast, N. Ireland Licensing provisions: University of Tennessee copyright agreement No. of lines in distributed program, including test data, etc.: 10 861 No. of bytes in distributed program, including test data, etc.: 394 631 Distribution format: tar.gz Programming language: C for processors, CUDA for NVIDIA GPUs Computer: Developed and tested on various x86 computers and NVIDIA C1060 Tesla and GTX 480 Fermi GPUs. The system targets x86 workstations, optionally with multicore processors or NVIDIA GPUs as accelerators. Operating system: Tested under Ubuntu Linux OS and CentOS 5.5 Linux OS Classification: 3, 16.12 Nature of problem: Simulation of chemical systems, particularly with low species populations, can be accurately performed using Gillespieʼs method of stochastic simulation. Numerous variations on the original stochastic simulation algorithm have been developed, including approaches that produce results with statistics that exactly match the chemical master equation (CME) as well as other approaches that approximate the CME. Solution method: The Accelerated Exact Stochastic Simulation (AESS) tool provides implementations of a wide variety of popular variations on the Gillespie method. Users can select the specific algorithm considered most appropriate. Comparisons between the methods and with other available implementations indicate that AESS provides the fastest known implementation of Gillespieʼs method for a variety of test models. Users may wish to execute ensembles of simulations to sweep parameters or to obtain better statistical results, so AESS supports acceleration of ensembles of simulation using parallel processing with MPI, SSE vector units on x86 processors, and/or using NVIDIA GPUs with CUDA.