Binomial distribution based tau-leap accelerated stochastic simulation

Recently, Gillespie introduced the tau-leap approximate, accelerated stochastic Monte Carlo method for well-mixed reacting systems [J. Chem. Phys. 115, 1716 (2001)]. In each time increment of that method, one executes a number of reaction events, selected randomly from a Poisson distribution, to enable simulation of long times. Here we introduce a binomial distribution tau-leap algorithm (abbreviated as BD-tau method). This method combines the bounded nature of the binomial distribution variable with the limiting reactant and constrained firing concepts to avoid negative populations encountered in the original tau-leap method of Gillespie for large time increments, and thus conserve mass. Simulations using prototype reaction networks show that the BD-tau method is more accurate than the original method for comparable coarse-graining in time.     

Binomial leap methods for simulating stochastic chemical kinetics

This paper discusses efficient simulation methods for stochastic chemical kinetics. Based on the tau-leap and midpoint tau-leap methods of Gillespie [D. T. Gillespie, J. Chem. Phys. 115, 1716 (2001)], binomial random variables are used in these leap methods rather than Poisson random variables. The motivation for this approach is to improve the efficiency of the Poisson leap methods by using larger stepsizes. Unlike Poisson random variables whose range of sample values is from zero to infinity, binomial random variables have a finite range of sample values. This probabilistic property has been used to restrict possible reaction numbers and to avoid negative molecular numbers in stochastic simulations when larger stepsize is used. In this approach a binomial random variable is defined for a single reaction channel in order to keep the reaction number of this channel below the numbers of molecules that undergo this reaction channel. A sampling technique is also designed for the total reaction number of a reactant species that undergoes two or more reaction channels. Samples for the total reaction number are not greater than the molecular number of this species. In addition, probability properties of the binomial random variables provide stepsize conditions for restricting reaction numbers in a chosen time interval. These stepsize conditions are important properties of robust leap control strategies. Numerical results indicate that the proposed binomial leap methods can be applied to a wide range of chemical reaction systems with very good accuracy and significant improvement on efficiency over existing approaches.     

Efficient stochastic simulation of simultaneous reaction and diffusion in a gas-liquid interface

A stochastic simulation of simultaneous reaction and diffusion is proposed for the gas-liquid interface formed in the surface of a gas bubble within a liquid. The interface between a carbon dioxide bubble and an aqueous solution of calcium hydroxide was simulated as an application example, taken from the integrated production of calcium carbonate. First Gillespie’s stochastic simulation algorithm was applied in separate reaction and diffusion simulations. The results from these simulations were consistent with deterministic solutions based on differential equations. However it was observed that stochastic diffusion simulations are extremely slow. The sampling of diffusion events was accelerated applying a group molecule transfer scheme based on the binomial distribution function. Simulations of the reaction-diffusion in the gas-liquid interface based on the standard Gillespie’s stochastic algorithm were also slow. However the application of the binomial distribution function scheme allowed to compute the concentration profiles in the gas-liquid interface in a fraction of the time required with the standard Gillespie’s stochastic algorithm.     

Look before you leap: a confidence-based method for selecting species criticality while avoiding negative populations in τ-leaping

The stochastic simulation algorithm was introduced by Gillespie and in a different form by Kurtz. There have been many attempts at accelerating the algorithm without deviating from the behavior of the simulated system. The crux of the explicit τ-leaping procedure is the use of Poisson random variables to approximate the number of occurrences of each type of reaction event during a carefully selected time period, τ. This method is acceptable providing the leap condition, that no propensity function changes "significantly" during any time-step, is met. Using this method there is a possibility that species numbers can, artificially, become negative. Several recent papers have demonstrated methods that avoid this situation. One such method classifies, as critical, those reactions in danger of sending species populations negative. At most, one of these critical reactions is allowed to occur in the next time-step. We argue that the criticality of a reactant species and its dependent reaction channels should be related to the probability of the species number becoming negative. This way only reactions that, if fired, produce a high probability of driving a reactant population negative are labeled critical. The number of firings of more reaction channels can be approximated using Poisson random variables thus speeding up the simulation while maintaining the accuracy. In implementing this revised method of criticality selection we make use of the probability distribution from which the random variable describing the change in species number is drawn. We give several numerical examples to demonstrate the effectiveness of our new method.     

A modified next reaction method for simulating chemical systems with time dependent propensities and delays

Chemical reaction systems with a low to moderate number of molecules are typically modeled as discrete jump Markov processes. These systems are oftentimes simulated with methods that produce statistically exact sample paths such as the Gillespie algorithm or the next reaction method. In this paper we make explicit use of the fact that the initiation times of the reactions can be represented as the firing times of independent, unit rate Poisson processes with internal times given by integrated propensity functions. Using this representation we derive a modified next reaction method and, in a way that achieves efficiency over existing approaches for exact simulation, extend it to systems with time dependent propensities as well as to systems with delays.     

Hybrid stochastic simulations of intracellular reaction–diffusion systems

With the observation that stochasticity is important in biological systems, chemical kinetics have begun to receive wider interest. While the use of Monte Carlo discrete event simulations most accurately capture the variability of molecular species, they become computationally costly for complex reaction–diffusion systems with large populations of molecules. On the other hand, continuous time models are computationally efficient but they fail to capture any variability in the molecular species. In this study a hybrid stochastic approach is introduced for simulating reaction–diffusion systems. We developed an adaptive partitioning strategy in which processes with high frequency are simulated with deterministic rate-based equations, and those with low frequency using the exact stochastic algorithm of Gillespie. Therefore the stochastic behavior of cellular pathways is preserved while being able to apply it to large populations of molecules. We describe our method and demonstrate its accuracy and efficiency compared with the Gillespie algorithm for two different systems. First, a model of intracellular viral kinetics with two steady states and second, a compartmental model of the postsynaptic spine head for studying the dynamics of Ca+2 and NMDA receptors.     

Unbiased tau-leap methods for stochastic simulation of chemically reacting systems

The tau-leap method first developed by Gillespie [D. T. Gillespie, J. Chem. Phys. 115, 1716 (2001)] can significantly speed up stochastic simulation of certain chemically reacting systems with acceptable losses in accuracy. Recently, several improved tau-leap methods, including the binomial, multinomial, and modified tau-leap methods, have been developed. However, in all these tau-leap methods, the mean of the number of times, K(m), that the mth reaction channel fires during a leap is not equal to the true mean. Therefore, all existing tau-leap methods produce biased simulation results, which limit the simulation accuracy and speed. In this paper, we analyze the mean of K(m) based on the chemical master equation. Using this analytical result, we develop unbiased Poisson and binomial tau-leap methods. Moreover, we analyze the variance of K(m), and then develop an unbiased Poisson/Gaussian/binomial tau-leap method to correct the errors in both the mean and variance of K(m). Simulation results demonstrate that our unbiased tau-leap method can significantly improve simulation accuracy without sacrificing speed.     

A kinetic Monte Carlo simulation study of inositol 1,4,5-trisphosphate receptor (IP3R) calcium release channel

Most of the previously theoretical studies about the stochastic nature of the IP3R calcium release channel gating use the chemical master equation (CME) approach. Because of the limitations of this approach we have used a stochastic simulation algorithm (SSA) presented by Gillespie. A single subunit of De Young-Keizer (DYK) model was simulated using Gillespie algorithm. The model has been considered in its complete form with eight states. We investigate the conditions which affect the open state of the model. Calcium concentrations were the subject of fluctuation in the previous works while in this study the population of the states is the subject of stochastic fluctuations. We found out that decreasing open probability is a function of Ca(2+) concentration in fast time domain, while in slow time domain it is a function of IP3 concentration. Studying the population of each state shows a time dependent reaction pattern in fast and medium time domains (10(-4) and 10(-3)s). In this pattern the state of X(010) has a determinative role in selecting the open state path. Also, intensity and frequency of fluctuations and Ca(2+) inhibitions have been studied. The results indicate that Gillespie algorithm can be a better choice for studying such systems, without using any approximation or elimination while having acceptable accuracy. In comparison with the chemical master equation, Gillespie algorithm is also provides a wide area for studying biological systems from other points of view.     

Entropy production of a mechanically driven single oligomeric enzyme: a consequence of fluctuation theorem

In this work we have shown how an applied mechanical force affects an oligomeric enzyme kinetics in a chemiostatic condition where the statistical characteristics of random walk of the substrate molecules over a finite number of active sites of the enzyme plays important contributing factors in governing the overall rate and nonequilibrium thermodynamic properties. The analytical results are supported by the simulation of single trajectory based approach of entropy production using Gillespie’s stochastic algorithm. This microscopic numerical approach not only gives the macroscopic entropy production from the mean of the distribution of entropy production which depends on the force but also a broadening of the distribution by the applied mechanical force, a kind of power broadening. In the nonequilibrium steady state (NESS), both the mean and the variance of the distribution increases and then saturates with the rise in applied force corresponding to the situation when the net rate of product formation reaches a limiting value with an activationless transition. The effect of the system-size and force on the entropy production distribution is shown to be constrained by the detailed fluctuation theorem.     

Simulating cumulative distributions under transient conditions in well-mixed,continuous and batch polymerization reactors

Instantaneous property methods are applied to the simulation of cumulative distributions of polymer chain architecture in well-mixed, continuous and batch reactors. We show how cumulative distributions can be simulated under transient conditions using a general convolution integral and suitable output from a dynamic reactor simulation package such as POLYRED^TM. We also show how the state equations that describe a polymerizing system can be augmented for direct simulation of these distributions. Examples of the impact of catalyst kinetics, reactor operating policies, and process upsets on the molecular weight distribution of bimodal high molecular weight polyethylene produced in well-mixed, stirred-bed gas-phase reactors are simulated for a binary transition-metal catalyst system.     

Discrete-time stochastic modeling and simulation of biochemical networks

Since inherent randomness in chemically reacting systems is evident, stochastic modeling and simulation are exceedingly important for investigating complex biological networks. Within the most common stochastic approach a network is modeled by a continuous-time Markov chain governed by the chemical master equation. We show how the continuous-time Markov chain can be converted to a stochastically identical discrete-time Markov chain and obtain a discrete-time version of the chemical master equation. Simulating the discrete-time Markov chain is equivalent to the Gillespie algorithm but requires less effort in that it eliminates the generation of exponential random variables. Thus, exactness as possessed by the Gillespie algorithm is preserved while the simulation can be performed more efficiently.     

Brownian dynamics of mixed surfactant micelles

We investigate micelle formation in a system containing two or more different amphiphiles with different geometries using a stochastic molecular-dynamics (MD) simulation method. For a binary system containing two amphiphiles, we calculate the critical micelle concentration (CMC) and cluster distribution for the mixture at several mole fractions and compare the simulation results with those predicted by analytic theories in the dilute limit and with experiments. We find that the CMC obtained from molecular mean-field theory agrees well with our simulation results. Motivated by the industrial use of mixed surfactant systems, we then extend our studies to a system containing six different chain lengths drawn from a Poisson distribution. We find that unlike a binary mixture of amphiphiles, the different species cancel the effects of each other so that the cluster distribution for the mixture has a shape of a system consisted entirely of amphiphiles of length equal to the mean chain length of the Poisson distribution.     

An accelerated algorithm for discrete stochastic simulation of reaction-diffusion systems using gradient-based diffusion and tau-leaping

Stochastic simulation of reaction-diffusion systems enables the investigation of stochastic events arising from the small numbers and heterogeneous distribution of molecular species in biological cells. Stochastic variations in intracellular microdomains and in diffusional gradients play a significant part in the spatiotemporal activity and behavior of cells. Although an exact stochastic simulation that simulates every individual reaction and diffusion event gives a most accurate trajectory of the system's state over time, it can be too slow for many practical applications. We present an accelerated algorithm for discrete stochastic simulation of reaction-diffusion systems designed to improve the speed of simulation by reducing the number of time-steps required to complete a simulation run. This method is unique in that it employs two strategies that have not been incorporated in existing spatial stochastic simulation algorithms. First, diffusive transfers between neighboring subvolumes are based on concentration gradients. This treatment necessitates sampling of only the net or observed diffusion events from higher to lower concentration gradients rather than sampling all diffusion events regardless of local concentration gradients. Second, we extend the non-negative Poisson tau-leaping method that was originally developed for speeding up nonspatial or homogeneous stochastic simulation algorithms. This method calculates each leap time in a unified step for both reaction and diffusion processes while satisfying the leap condition that the propensities do not change appreciably during the leap and ensuring that leaping does not cause molecular populations to become negative. Numerical results are presented that illustrate the improvement in simulation speed achieved by incorporating these two new strategies.     

K-leap method for accelerating stochastic simulation of coupled chemical reactions

Leap methods are very promising for accelerating stochastic simulation of a well stirred chemically reacting system, while providing acceptable simulation accuracy. In Gillespie's tau-leap method [D. Gillespie, J. Phys. Chem. 115, 1716 (2001)], the number of firings of each reaction channel during a leap is a Poisson random variable, whose sample values are unbounded. This may cause large changes in the populations of certain molecular species during a leap, thereby violating the leap condition. In this paper, we develop an alternative leap method called the K-leap method, in which we constrain the total number of reactions occurring during a leap to be a number K calculated from the leap condition. As the number of firings of each reaction channel during a leap is upper bounded by a properly chosen number, our K-leap method can better satisfy the leap condition, thereby improving simulation accuracy. Since the exact stochastic simulation algorithm (SSA) is a special case of our K-leap method when K=1, our K-leap method can naturally change from the exact SSA to an approximate leap method during simulation, whenever the leap condition allows to do so.     

Replica-exchange extensions of simulated tempering method

In this paper we consider combinations of two well-known generalized-ensemble algorithms, namely, simulated tempering and replica-exchange method. We discuss two examples of such combinations. One is the replica-exchange simulated tempering and the other is the simulated tempering replica-exchange method. In the former method, a short replica-exchange simulation is first performed and the simulated tempering weight factor is obtained by the multiple-histogram reweighting techniques. This process of simulated tempering weight factor determination is faster and simpler than that in the usual iterative process. A long simulated tempering production run is then performed with this weight factor. The latter method is a further extension of the former in which a simulated tempering replica-exchange simulation is performed with a small number of replicas. These algorithms are particularly useful for studying frustrated systems with rough energy landscape. We give the formulations of these two methods in detail and demonstrate their effectiveness taking the example of the system of a 17-residue helical peptide.     

Further on creating normal density functions in Microsoft<Superscript>®</Superscript> Excel

Simulations play an increasing role in metrology and analytical chemistry, particularly the use of simulated normal distributions. Microsoft Excel is available to almost everybody and normal distributions can be simulated using either an innate function or other algorithms. We explore the success of five different algorithms to simulate normal distributions and compare the outcome with a simulation coded in R. The normality of 10⁶ data points simulated by the chosen algorithms was tested by different recognized statistical procedures and Q–Q plots. They all failed the statistical procedures, but evaluating the Q–Q plots revealed different types of deviations. It also showed that the deviations, although statistically significant, most likely do not have any practical relevance in laboratory work.     

激光溅射下原子团簇生长的非平衡动力学

在“半球模型”的基础上,在考虑了因等离子体膨胀和热辐射等引起的温度答体积的变化的情况下,进一步考虑了碳簇的链状、单环、双环、多环及富勒烯等五种结构及中性粒子与中性粒子、中性粒子与离子两种反应,计算了碳簇中性分子和离子的形成动力学及其尺寸分布。     

Dynamic Monte Carlo simulation of chain growth polymerization and its concentration effect

Polymerization kinetics and chain configuration are two of lasting hotspots in Polymer Science. As these problems are concerned, computer experiment affords a very useful method besides theoretical and experi-mental investigations. The conventional Monte …     

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 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. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2012     

Numerical Study of Pore Sizes Distribution in Gels

We introduce a new numerical technique for the calculation of the pore size distribution in two-dimensional disordered systems. Our method is based on a triangulation technique which allows a closer measurement of pores surface without any morphological hypothesis.In this work, we focus our calculations in simulated gels. Such materials are modeled in two different conditions: by means of the Diffusion-Limited and Reaction-Limited Cluster-cluster Aggregation algorithms, DLCA and RLCA, respectively. In both situations, when the particles concentration decreases, the average pores size increases. The more compact cluster in RLCA, compared with DLCA, is consistent with the pore size distribution we have calculated. The simulated mean pore size is quantitatively in agreement with experimental data from literature.