Stochastic Simulation Algorithm Research Articles

A multi-time-scale analysis of chemical reaction networks: II. Stochastic systems.

We consider stochastic descriptions of chemical reaction networks in which there are both fast and slow reactions, and for which the time scales are widely separated. We develop a computational algorithm that produces the generator of the full chemical master equation for arbitrary systems, and show how to obtain a reduced equation that governs the evolution on the slow time scale. This is done by applying a state space decomposition to the full equation that leads to the reduced dynamics in terms of certain projections and the invariant distributions of the fast system. The rates or propensities of the reduced system are shown to be the rates of the slow reactions conditioned on the expectations of fast steps. We also show that the generator of the reduced system is a Markov generator, and we present an efficient stochastic simulation algorithm for the slow time scale dynamics. We illustrate the numerical accuracy of the approximation by simulating several examples. Graph-theoretic techniques are used throughout to describe the structure of the reaction network and the state-space transitions accessible under the dynamics.

Brownian-motion based simulation of stochastic reaction–diffusion systems for affinity based sensors

In this work, we develop a 2D algorithm for stochastic reaction–diffusion systems describing the binding and unbinding of target molecules at the surfaces of affinity-based sensors. In particular, we simulate the detection of DNA oligomers using silicon-nanowire field-effect biosensors. Since these devices are uniform along the nanowire, two dimensions are sufficient to capture the kinetic effects features. The model combines a stochastic ordinary differential equation for the binding and unbinding of target molecules as well as a diffusion equation for their transport in the liquid. A Brownian-motion based algorithm simulates the diffusion process, which is linked to a stochastic-simulation algorithm for association at and dissociation from the surface. The simulation data show that the shape of the cross section of the sensor yields areas with significantly different target-molecule coverage. Different initial conditions are investigated as well in order to aid rational sensor design. A comparison of the association/hybridization behavior for different receptor densities allows optimization of the functionalization setup depending on the target-molecule density.

A Survey of Stochastic Simulation and Optimization Methods in Signal Processing

Modern signal processing (SP) methods rely very heavily on probability and statistics to solve challenging SP problems. SP methods are now expected to deal with ever more complex models, requiring ever more sophisticated computational inference techniques. This has driven the development of statistical SP methods based on stochastic simulation and optimization. Stochastic simulation and optimization algorithms are computationally intensive tools for performing statistical inference in models that are analytically intractable and beyond the scope of deterministic inference methods. They have been recently successfully applied to many difficult problems involving complex statistical models and sophisticated (often Bayesian) statistical inference techniques. This survey paper offers an introduction to stochastic simulation and optimization methods in signal and image processing. The paper addresses a variety of high-dimensional Markov chain Monte Carlo (MCMC) methods as well as deterministic surrogate methods, such as variational Bayes, the Bethe approach, belief and expectation propagation and approximate message passing algorithms. It also discusses a range of optimization methods that have been adopted to solve stochastic problems, as well as stochastic methods for deterministic optimization. Subsequently, areas of overlap between simulation and optimization, in particular optimization-within-MCMC and MCMC-driven optimization are discussed.

An adaptive tau-leaping method for stochastic simulations of reaction-diffusion systems

Stochastic modelling is critical for studying many biochemical processes in a cell, in particular when some reacting species have low population numbers. For many such cellular processes the spatial distribution of the molecular species plays a key role. The evolution of spatially heterogeneous biochemical systems with some species in low amounts is accurately described by the mesoscopic model of the Reaction-Diffusion Master Equation. The Inhomogeneous Stochastic Simulation Algorithm provides an exact strategy to numerically solve this model, but it is computationally very expensive on realistic applications. We propose a novel adaptive time-stepping scheme for the tau-leaping method for approximating the solution of the Reaction-Diffusion Master Equation. This technique combines effective strategies for variable time-stepping with path preservation to reduce the computational cost, while maintaining the desired accuracy. The numerical tests on various examples arising in applications show the improved efficiency achieved by the new adaptive method.

Efficient Constant-Time Complexity Algorithm for Stochastic Simulation of Large Reaction Networks.

Exact stochastic simulation is an indispensable tool for a quantitative study of biochemical reaction networks. The simulation realizes the time evolution of the model by randomly choosing a reaction to fire and update the system state according to a probability that is proportional to the reaction propensity. Two computationally expensive tasks in simulating large biochemical networks are the selection of next reaction firings and the update of reaction propensities due to state changes. We present in this work a new exact algorithm to optimize both of these simulation bottlenecks. Our algorithm employs the composition-rejection on the propensity bounds of reactions to select the next reaction firing. The selection of next reaction firings is independent of the number reactions while the update of propensities is skipped and performed only when necessary. It therefore provides a favorable scaling for the computational complexity in simulating large reaction networks. We benchmark our new algorithm with the state of the art algorithms available in literature to demonstrate its applicability and efficiency.

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.

Stochastic Focusing and Defocusing in Biological Reaction Networks: Lessons Learned from Accurate Chemical Master Equation (ACME) Solutions

Biological reaction networks are stochastic due to random thermal fluctuations. Stochasticity plays important roles in regulation of biochemical reaction networks when the copy numbers of molecular species are small and often results in unexpected outcomes. For example, it has been shown that a basic enzymatic reaction system can display stochastic focusing (SF) by increasing the sensitivity of the network as a result of the increasing signal noise [1]. Although stochastic simulation algorithm has been widely used to study such systems, it is ineffective in examining rare events and this becomes a significant issue when the tails of probability distributions are relevant as is the case of SF. Here we use the ACME method for the exact solution of the discrete Chemical Master Equation and study the probability landscape of product molecules in the basic enzymatic reaction system used in the original SF study [1]. Examinations of the effects of signal molecules under different stochastic processes show that SF is at play as stochastic changes enhance the system sensitivity. However, we also observed that the noise in signaling under certain stochastic processes in the same reaction network lead to a decrease in the system sensitivities, thus the network experiences stochastic defocusing. We further show that signal molecules following certain stochastic processes in the same reaction network can give rise to noise-induced bistability in the distribution of product molecules. These results highlight the fundamental role of stochasticity in biological reaction networks and the need for exact computation of probability landscape of the molecules in the system. It also points to possible importance of positive and negative feedback loops in such networks for control of the intrinsic noise. [1] Paulsson et. al, 2000, PNAS.

CERENA: ChEmical REaction Network Analyzer--A Toolbox for the Simulation and Analysis of Stochastic Chemical Kinetics.

Gene expression, signal transduction and many other cellular processes are subject to stochastic fluctuations. The analysis of these stochastic chemical kinetics is important for understanding cell-to-cell variability and its functional implications, but it is also challenging. A multitude of exact and approximate descriptions of stochastic chemical kinetics have been developed, however, tools to automatically generate the descriptions and compare their accuracy and computational efficiency are missing. In this manuscript we introduced CERENA, a toolbox for the analysis of stochastic chemical kinetics using Approximations of the Chemical Master Equation solution statistics. CERENA implements stochastic simulation algorithms and the finite state projection for microscopic descriptions of processes, the system size expansion and moment equations for meso- and macroscopic descriptions, as well as the novel conditional moment equations for a hybrid description. This unique collection of descriptions in a single toolbox facilitates the selection of appropriate modeling approaches. Unlike other software packages, the implementation of CERENA is completely general and allows, e.g., for time-dependent propensities and non-mass action kinetics. By providing SBML import, symbolic model generation and simulation using MEX-files, CERENA is user-friendly and computationally efficient. The availability of forward and adjoint sensitivity analyses allows for further studies such as parameter estimation and uncertainty analysis. The MATLAB code implementing CERENA is freely available from http://cerenadevelopers.github.io/CERENA/.

SOLVING STOCHASTIC CHEMICAL KINETICS BY METROPOLIS-HASTINGS SAMPLING

This study considers using Metropolis-Hastings algorithm for stochastic simulation of chemical reactions. The proposed method uses SSA (Stochastic Simulation Algorithm) distribution which is a standard method for solving well-stirred chemically reacting systems as a desired distribution. A new numerical solvers based on exponential form of exact and approximate solutions of CME (Chemical Master Equation) is employed for obtaining target and proposal distributions in Metropolis-Hastings algorithm to accelerate the accuracy of the tau-leap method. Samples generated by this technique have the same distribution as SSA and the histogram of samples show it's convergence to SSA.

Inflow performance relationships for layered solution-gas drive reservoir

Inflow performance relationship (IPR) is a very important tool to forecast well performance. Existing IPR models are idealistic since they are developed for homogeneous reservoirs; therefore, they are inappropriate for layered systems. Consequently, there is a need for IPR models that efficiently describe layered reservoir performance. This study investigates the effects of reservoir heterogeneity on IPR for layered solution-gas drive reservoirs. Multiphase flow in both two and multilayer reservoirs was simulated. Both fluid cross flow and no fluid cross flow among layers were considered. A stochastic simulation algorithm was used to generate various permeability realisations among layers. Three geostatistical models using uniform, Gaussian, and bimodal probability distributions were used to grasp optimum match between real reservoir behaviour and simulated data. The generated data were scrutinised to develop two accurate IPR equations. The first equation describes the well behaviour under current flowing conditions, whereas the second equation can be used to forecast future well performance.

Optimization for Service Routes of Pallet Service Center Based on the Pallet Pool Mode.

Service routes optimization (SRO) of pallet service center should meet customers' demand firstly and then, through the reasonable method of lines organization, realize the shortest path of vehicle driving. The routes optimization of pallet service center is similar to the distribution problems of vehicle routing problem (VRP) and Chinese postman problem (CPP), but it has its own characteristics. Based on the relevant research results, the conditions of determining the number of vehicles, the one way of the route, the constraints of loading, and time windows are fully considered, and a chance constrained programming model with stochastic constraints is constructed taking the shortest path of all vehicles for a delivering (recycling) operation as an objective. For the characteristics of the model, a hybrid intelligent algorithm including stochastic simulation, neural network, and immune clonal algorithm is designed to solve the model. Finally, the validity and rationality of the optimization model and algorithm are verified by the case.

Stochastic Simulation of a CumSum Knock Controller

Single time-history simulations of knock control strategies cannot provide repeatable measures of closed loop performance because such responses are strongly dependent on the random arrival of knock events during the specific simulation or experiment. Stochastic simulation algorithms address this issue by computing the probability distribution of the closed loop response as a function of time from any given initial distribution. In this paper, a new algorithm is presented for stochastic simulation of a Cumulative Summation-based knock control strategy. The stochastic simulator is then used to assess the performance of this strategy in a scenario in which feedback sensitivity is also enhanced through the use of an optimized knock threshold. The results are compared to those for a traditional knock controller operating using a more conventional knock threshold level.

Fock space, symbolic algebra, and analytical solutions for small stochastic systems.

Randomness is ubiquitous in nature. From single-molecule biochemical reactions to macroscale biological systems, stochasticity permeates individual interactions and often regulates emergent properties of the system. While such systems are regularly studied from a modeling viewpoint using stochastic simulation algorithms, numerous potential analytical tools can be inherited from statistical and quantum physics, replacing randomness due to quantum fluctuations with low-copy-number stochasticity. Nevertheless, classical studies remained limited to the abstract level, demonstrating a more general applicability and equivalence between systems in physics and biology rather than exploiting the physics tools to study biological systems. Here the Fock space representation, used in quantum mechanics, is combined with the symbolic algebra of creation and annihilation operators to consider explicit solutions for the chemical master equations describing small, well-mixed, biochemical, or biological systems. This is illustrated with an exact solution for a Michaelis-Menten single enzyme interacting with limited substrate, including a consideration of very short time scales, which emphasizes when stiffness is present even for small copy numbers. Furthermore, we present a general matrix representation for Michaelis-Menten kinetics with an arbitrary number of enzymes and substrates that, following diagonalization, leads to the solution of this ubiquitous, nonlinear enzyme kinetics problem. For this, a flexible symbolic maple code is provided, demonstrating the prospective advantages of this framework compared to stochastic simulation algorithms. This further highlights the possibilities for analytically based studies of stochastic systems in biology and chemistry using tools from theoretical quantum physics.

Hybrid stochastic simulation of reaction-diffusion systems with slow and fast dynamics.

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

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.

Simulation of AMPA and NMDA contribution to postsynaptic response

The synaptic response depends on the dynamics of the released neurotransmitter, as well as on the ability of ionic receptors (AMPA and NMDA) to activate and deactivate in time. Simulation of these receptors allow to study interesting patterns as neuronal plasticity, which modulates the efficiency of transmission in two different manners: facilitation (incremented response to weak stimuli), and depression (decreased response to strong stimuli). In this work, we have simulated the postsynaptic response due to AMPA and NMDA activation. Our model considers the postsynaptic membrane as a disk of 0.1 microns of radius where receptors lie with a ratio of 4:1 (80 AMPA and 20 NMDA) [1]. They are uniformly distributed in the whole area of the disk, assuming they have the same probability of being activated as they come in contact with neurotransmitter molecules which are dispersed homogeneously over the postsynaptic membrane. We use the kinetic models proposed by [2] and [3] for AMPA and NMDA receptors, respectively. The model was solved with a stochastic simulation algorithm, and plots show the mean values calculated from 500 runs. Figures ​Figures1A1A and ​and1B1B show the results obtained when the neurotransmitter comes in the form of two pulses of different concentrations in a short period of time. Figure ​Figure1A1A corresponds to the case of low concentration pulses (40 and 160 molecules for the first and second pulses, respectively). In this simulation, it can be seen how AMPA receptors increase their response to the second stimulus, exhibiting facilitation. Figure ​Figure1B1B shows the case of high concentration pulses (1207 and 4787 molecules, respectively). In this simulation, a desensitized response is observed, that is, a depression effect which reduces the number of open receptors for the second stimulus. In both cases it is observed that the contribution made by the NMDA receptors is identical, whereas AMPA receptors have a contribution that follows the rapid dynamics of the neurotransmitter released by the presynaptic neuron. Figure 1 Simulations A Low-concentration pulses. B High-concentration pulses. Above input trains (red) and below AMPA (green), NMDA (magenta) and Total receptors (black) responses. Other simulations for different trains of neurotransmitter were performed, and similar results were obtained. We also performed simulations for different stimulation frequencies in larger timescales and again, AMPA receptors followed rapid neurotransmitter dynamics while NMDA receptors exhibited changes only in time activation. Then, we conclude that the contribution of AMPA receptors is the most important one in short timescales, whereas NMDA receptors participate mainly in larger timescale responses. These results are the basis of our current work on proposing a postsynaptic model intended to work with our presynaptic model [4].

Abstract A03: A brain cancer cell migration simulator based on a motor-clutch model

Abstract WHO Grade IV astrocytoma, also known as glioblastoma, is an invasive disease where cell migration drives dispersion of tumor cells throughout the brain. To predict cell migration speeds in complex mechanochemical environments, we developed a cell migration simulator based on our published motor-clutch framework [1-3]. One of the key predictions of our previous modeling is that cell traction force is biphasic with respect to environmental stiffness, i.e. there is an optimal stiffness at which cell traction force is maximal. To test whether such optimality exists within a whole cell migration simulator, we extended our model to allow for multiple cellular protrusions in 1D, 2D, or 3D environments. The cell migration simulator imposes that forces balance and that masses of key elements- actin, adhesions, and myosin motors- are conserved. Using the Gillespie stochastic simulation algorithm implemented in matlab, we simulated cell traction force, F-actin retrograde flow, cell aspect ratio, and cell migration as a function of cell stiffness. We found that all four variables depends biphasically on simulated environmental stiffness with approximately overlapping optima. In addition, the model predicts that the optimal stiffness can be shifted by coordinate shifting of motor and adhesive clutch numbers, a prediction that we are now testing experimentally in vitro. Overall, our modeling results predict that effective anti-migratory treatment of patients will require knowledge of the cell migration speed landscape within parameter space. In particular, we predict that simultaneously inhibiting motors and clutches will be beneficial for patients whose tumors operate near the optimum, largely ineffective for those below the optimum, and potentially detrimental to those above the optimum. Thus, our cell migration simulator provides dynamic physical biomarkers that can potentially be used to classify glioblastoma patients for more precise treatment.

Comparison of different moment-closure approximations for stochastic chemical kinetics.

In recent years, moment-closure approximations (MAs) of the chemical master equation have become a popular method for the study of stochastic effects in chemical reaction systems. Several different MA methods have been proposed and applied in the literature, but it remains unclear how they perform with respect to each other. In this paper, we study the normal, Poisson, log-normal, and central-moment-neglect MAs by applying them to understand the stochastic properties of chemical systems whose deterministic rate equations show the properties of bistability, ultrasensitivity, and oscillatory behaviour. Our results suggest that the normal MA is favourable over the other studied MAs. In particular, we found that (i) the size of the region of parameter space where a closure gives physically meaningful results, e.g., positive mean and variance, is considerably larger for the normal closure than for the other three closures, (ii) the accuracy of the predictions of the four closures (relative to simulations using the stochastic simulation algorithm) is comparable in those regions of parameter space where all closures give physically meaningful results, and (iii) the Poisson and log-normal MAs are not uniquely defined for systems involving conservation laws in molecule numbers. We also describe the new software package MOCA which enables the automated numerical analysis of various MA methods in a graphical user interface and which was used to perform the comparative analysis presented in this paper. MOCA allows the user to develop novel closure methods and can treat polynomial, non-polynomial, as well as time-dependent propensity functions, thus being applicable to virtually any chemical reaction system.

A fast exact simulation method for a class of Markov jump processes.

A new method of the stochastic simulation algorithm (SSA), named the Hashing-Leaping method (HLM), for exact simulations of a class of Markov jump processes, is presented in this paper. The HLM has a conditional constant computational cost per event, which is independent of the number of exponential clocks in the Markov process. The main idea of the HLM is to repeatedly implement a hash-table-like bucket sort algorithm for all times of occurrence covered by a time step with length τ. This paper serves as an introduction to this new SSA method. We introduce the method, demonstrate its implementation, analyze its properties, and compare its performance with three other commonly used SSA methods in four examples. Our performance tests and CPU operation statistics show certain advantages of the HLM for large scale problems.

Relative frequencies of constrained events in stochastic processes: An analytical approach.

The stochastic simulation algorithm (SSA) and the corresponding Monte Carlo (MC) method are among the most common approaches for studying stochastic processes. They relies on knowledge of interevent probability density functions (PDFs) and on information about dependencies between all possible events. Analytical representations of a PDF are difficult to specify in advance, in many real life applications. Knowing the shapes of PDFs, and using experimental data, different optimization schemes can be applied in order to evaluate probability density functions and, therefore, the properties of the studied system. Such methods, however, are computationally demanding, and often not feasible. We show that, in the case where experimentally accessed properties are directly related to the frequencies of events involved, it may be possible to replace the heavy Monte Carlo core of optimization schemes with an analytical solution. Such a replacement not only provides a more accurate estimation of the properties of the process, but also reduces the simulation time by a factor of order of the sample size (at least ≈10(4)). The proposed analytical approach is valid for any choice of PDF. The accuracy, computational efficiency, and advantages of the method over MC procedures are demonstrated in the exactly solvable case and in the evaluation of branching fractions in controlled radical polymerization (CRP) of acrylic monomers. This polymerization can be modeled by a constrained stochastic process. Constrained systems are quite common, and this makes the method useful for various applications.