Stochastic simulation methods for biochemical models with multi-state species

Markovian pure jump processes model a wide range of phenomena, including chemical reactions at the molecular level, dynamics of wireless communication networks, and the spread of epidemic diseases in small populations. There exist algorithms such as Gillespie's stochastic simulation algorithm (SSA) and Anderson's modified next reaction method (MNRM) that simulate a single path with the exact distribution of the process, but this can be time consuming when many reactions take place during a short time interval. Gillespie's approximated tau-leap method, on the other hand, can be used to reduce computational time, but it may lead to nonphysical values due to a positive one-step exit probability, and it also introduces a time discretization error. Here, we present a novel hybrid algorithm for simulating individual paths which adaptively switches between the SSA and the tau-leap method. The switching strategy is based on a comparison of the expected interarrival time of the SSA and an adaptive time step derive...

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.

Chapter 7 - Monte Carlo stochastic simulation of the master equation for unimolecular reaction systems

Fine-Grained Parallel and Distributed Spatial Stochastic Simulation of Biological Reactions

Kinetic Monte CarloA procedure for simulating time evolution of a system define by a specified state, with a given set of transition probabilities between states. Spatially-homogeneous: chemical reactions in the gas phase of well-mixed systems (Gillespie's Stochastic Simulation Algorithm) Non-uniform structures: e.g., catalytic surface kinetics, applications involving diffusion and reaction processes in condensed matterHomogeneous Models X i (t): number of species S i at time t (i = 1, . . ., N). X(t) = (X 1 (t), . . ., X N (t)) T -column-vector Markov stochastic process.-the probability that, at time t, there are k 1 units of species S 1 , k 2 units of species S 2 , ..., given X(t 0 ) = k 0 . a j : Z N + R + -propensity function for the respective reaction R j . a j (k)dt -given X(t) = k, the probability that reaction R j takes place in the next infinitesimal time interval [t, t + dt). v j -the jth column of stoichiometry matrix V .

Stochastic models for chemical reaction networks have become very popular in recent years. For such models, the estimation of parameter sensitivities is an important and challenging problem. Sensitivity values help in analyzing the network, understanding its robustness properties and also in identifying the key reactions for a given outcome. Most of the methods that exist in the literature for the estimation of parameter sensitivities, rely on Monte Carlo simulations using Gillespie's stochastic simulation algorithm or its variants. It is well-known that such simulation methods can be prohibitively expensive when the network contains reactions firing at different time-scales, which is a feature of many important biochemical networks. For such networks, it is often possible to exploit the time-scale separation and approximately capture the original dynamics by simulating a model, which is obtained by eliminating the fast reactions in a certain way. The aim of this paper is to tie these model reduction techniques with sensitivity analysis. We prove that under some conditions, the sensitivity values for the reduced model can be used to approximately recover the sensitivity values for the original model. Through an example we illustrate how our result can help in sharply reducing the computational costs for the estimation of parameter sensitivities for reaction networks with multiple time-scales. To prove our result, we use coupling arguments based on the random time change representation of Kurtz. We also exploit certain connections between the distributions of the occupation times of Markov chains and multi-dimensional wave equations.

Antagonism and bistability in protein interaction networks

Noise Treatment in Models of Genetic Switches

BackgroundBiochemical systems with relatively low numbers of components must be simulated stochastically in order to capture their inherent noise. Although there has recently been considerable work on discrete stochastic solvers, there is still a need for numerical methods that are both fast and accurate. The Bulirsch-Stoer method is an established method for solving ordinary differential equations that possesses both of these qualities.ResultsIn this paper, we present the Stochastic Bulirsch-Stoer method, a new numerical method for simulating discrete chemical reaction systems, inspired by its deterministic counterpart. It is able to achieve an excellent efficiency due to the fact that it is based on an approach with high deterministic order, allowing for larger stepsizes and leading to fast simulations. We compare it to the Euler τ-leap, as well as two more recent τ-leap methods, on a number of example problems, and find that as well as being very accurate, our method is the most robust, in terms of efficiency, of all the methods considered in this paper. The problems it is most suited for are those with increased populations that would be too slow to simulate using Gillespie’s stochastic simulation algorithm. For such problems, it is likely to achieve higher weak order in the moments.ConclusionsThe Stochastic Bulirsch-Stoer method is a novel stochastic solver that can be used for fast and accurate simulations. Crucially, compared to other similar methods, it better retains its high accuracy when the timesteps are increased. Thus the Stochastic Bulirsch-Stoer method is both computationally efficient and robust. These are key properties for any stochastic numerical method, as they must typically run many thousands of simulations.

The budding yeast cell cycle provides an excellent example of the need for modeling stochastic effects in mathematical modeling of biochemical reactions. The continuous deterministic approach using ordinary differential equations is adequate for understanding the average behavior of cells, while the discrete stochastic approach accurately captures noisy events in the growth-division cycle. This paper presents a stochastic approximation of the cell cycle for budding yeast using Gillespie's stochastic simulation algorithm. To compare the stochastic results with the average behavior, the simulation must be run thousands of times. A load balancing algorithm reduces the cost for making those runs by 14% when run on a parallel supercomputer.

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.

Stochastic simulation is an essential method for modeling complex geological structures in geosciences. Evaluating the uncertainty of the realizations of stochastic simulations can better describe real phenomena. However, uncertainty evaluation of stochastic simulation methods remains a challenge due to the limited data from geological surveys and the uncertainty in reliability estimation with stochastic simulation models. In addition, understanding the sensitivity of the parameters in stochastic simulation models is invaluable when exploring the parameters with a higher influence on the uncertainty associated with predictions generated from stochastic simulation. To facilitate uncertainty evaluation in stochastic simulation methods, we use the circular treemap as an interactive workflow to explore prediction uncertainty in and the parameter sensitivity of multiple-point geostatistical (MPS) stochastic simulation methods. In this work, we present a novel visualization framework for assessing the uncertainty in MPS stochastic simulation methods and exploring the parameter sensitivity of the MPS methods. We present a new indicator to integrate multiple metrics that characterize geospatial features and visualize these metrics to assist domain experts in making decisions. Parallel coordinates-scatter matrix plots and multi-dimensional scaling (MDS) plots are used to analyze the parametric sensitivity of MPS stochastic simulation methods. The realizations and parameters of two MPS stochastic simulation methods are used to test the applicability of the proposed visualization workflow and the visualization methods. The results demonstrate that our workflow and the visualization methods can assist experts in finding the model with less uncertainty and improve the efficiency of parameter adjustment using different MPS stochastic simulation methods.

The accepted stochastic descriptions of biochemical dynamics under well-mixed conditions are given by the Chemical Master Equation and the Stochastic Simulation Algorithm, which are equivalent. The latter is a Monte-Carlo method, which, despite enjoying broad availability in a large number of existing software packages, is computationally expensive due to the huge amounts of ensemble averaging required for obtaining accurate statistical information. The former is a set of coupled differential-difference equations for the probability of the system being in any one of the possible mesoscopic states; these equations are typically computationally intractable because of the inherently large state space. Here we introduce the software package intrinsic Noise Analyzer (iNA), which allows for systematic analysis of stochastic biochemical kinetics by means of van Kampen’s system size expansion of the Chemical Master Equation. iNA is platform independent and supports the popular SBML format natively. The present implementation is the first to adopt a complementary approach that combines state-of-the-art analysis tools using the computer algebra system Ginac with traditional methods of stochastic simulation. iNA integrates two approximation methods based on the system size expansion, the Linear Noise Approximation and effective mesoscopic rate equations, which to-date have not been available to non-expert users, into an easy-to-use graphical user interface. In particular, the present methods allow for quick approximate analysis of time-dependent mean concentrations, variances, covariances and correlations coefficients, which typically outperforms stochastic simulations. These analytical tools are complemented by automated multi-core stochastic simulations with direct statistical evaluation and visualization. We showcase iNA’s performance by using it to explore the stochastic properties of cooperative and non-cooperative enzyme kinetics and a gene network associated with circadian rhythms. The software iNA is freely available as executable binaries for Linux, MacOSX and Microsoft Windows, as well as the full source code under an open source license.

A belief network comprises a graphical representation of dependencies between variables of a domain and a set of conditional probabilities associated with each dependency. Unless rho =NP, an efficient, exact algorithm does not exist to compute probabilistic inference in belief networks. Stochastic simulation methods, which often improve run times, provide an alternative to exact inference algorithms. Such a stochastic simulation algorithm, D-BNRAS, which is a randomized approximation scheme is presented. To analyze the run time, belief networks are parameterized, by the dependence value D/sub xi /, which is a measure of the cumulative strengths of the belief network dependencies given background evidence xi . This parameterization defines the class of f-dependence networks. The run time of D-BNRAS is polynomial when f is a polynomial function. Thus, the results prove the existence of a class of belief networks for which inference approximation is polynomial and, hence, provably faster than any exact algorithm. >

Journal of Computational ChemistryVolume 34, Issue 5 p. iii-iv Cover ImageFree Access Inside Cover, Volume 34, Issue 5 First published: 23 January 2013 https://doi.org/10.1002/jcc.23240AboutPDF ToolsRequest permissionExport citationAdd to favoritesTrack citation ShareShare Give accessShare full text accessShare full-text accessPlease review our Terms and Conditions of Use and check box below to share full-text version of article.I have read and accept the Wiley Online Library Terms and Conditions of UseShareable LinkUse the link below to share a full-text version of this article with your friends and colleagues. Learn more.Copy URL Share a linkShare onFacebookTwitterLinked InRedditWechat Abstract Many approximate accelerated stochastic simulation methods exist that can provide a means of employing Kinetic Monte Carlo techniques to model chemical kinetics. The drawback in using such methods is that sometimes the population of the reactant species present during the simulations becomes negative, leading to unphysical results. This is due to the manner in which the change in time during the reactions is handled by the stochastic simulation methods in such systems. The work by Shantanu Kadam and Kumar Vanka on page 394 discusses a means of avoiding this problem by combining the representative reaction approach (RRA) based methodology with the stochastic simulation algorithm (SSA) and the binomial method. Volume34, Issue515 February 2013Pages iii-iv RelatedInformation