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.     

Avoiding negative populations in explicit Poisson tau-leaping

The explicit tau-leaping procedure attempts to speed up the stochastic simulation of a chemically reacting system by approximating the number of firings of each reaction channel during a chosen time increment tau as a Poisson random variable. Since the Poisson random variable can have arbitrarily large sample values, there is always the possibility that this procedure will cause one or more reaction channels to fire so many times during tau that the population of some reactant species will be driven negative. Two recent papers have shown how that unacceptable occurrence can be avoided by replacing the Poisson random variables with binomial random variables, whose values are naturally bounded. This paper describes a modified Poisson tau-leaping procedure that also avoids negative populations, but is easier to implement than the binomial procedure. The new Poisson procedure also introduces a second control parameter, whose value essentially dials the procedure from the original Poisson tau-leaping at one extreme to the exact stochastic simulation algorithm at the other; therefore, the modified Poisson procedure will generally be more accurate than the original Poisson procedure.     

Efficient binomial leap method for simulating chemical kinetics

The binomial tau-leaping method of simulating the stochastic time evolution in a reaction system uses a binomial random number to approximate the number of reaction events. Theory implies that this method can avoid negative molecular numbers in stochastic simulations when a larger time step tau is used. Presented here is a modified binomial leap method for improving the accuracy and application range of the binomial leap method. The maximum existing population is first defined in this approach in order to determine a better bound of the number reactions. To derive a general leap procedure in chemically reacting systems, in this method a new sampling procedure based on the species is also designed for the maximum bound of consumed molecules of a reactant species in reaction channel. Numerical results indicate that the modified binomial leap method can be applied to a wider application range of chemically reacting systems with much better accuracy than the existing binomial leap method.     

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     

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.     

R-leaping: accelerating the stochastic simulation algorithm by reaction leaps

A novel algorithm is proposed for the acceleration of the exact stochastic simulation algorithm by a predefined number of reaction firings (R-leaping) that may occur across several reaction channels. In the present approach, the numbers of reaction firings are correlated binomial distributions and the sampling procedure is independent of any permutation of the reaction channels. This enables the algorithm to efficiently handle large systems with disparate rates, providing substantial computational savings in certain cases. Several mechanisms for controlling the accuracy and the appearance of negative species are described. The advantages and drawbacks of R-leaping are assessed by simulations on a number of benchmark problems and the results are discussed in comparison with established methods.     

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.     

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.     

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.     

Integral tau methods for stiff stochastic chemical systems

Tau leaping methods enable efficient simulation of discrete stochastic chemical systems. Stiff stochastic systems are particularly challenging since implicit methods, which are good for stiffness, result in noninteger states. The occurrence of negative states is also a common problem in tau leaping. In this paper, we introduce the implicit Minkowski-Weyl tau (IMW-τ) methods. Two updating schemes of the IMW-τ methods are presented: implicit Minkowski-Weyl sequential (IMW-S) and implicit Minkowski-Weyl parallel (IMW-P). The main desirable feature of these methods is that they are designed for stiff stochastic systems with molecular copy numbers ranging from small to large and that they produce integer states without rounding. This is accomplished by the use of a split step where the first part is implicit and computes the mean update while the second part is explicit and generates a random update with the mean computed in the first part. We illustrate the IMW-S and IMW-P methods by some numerical examples, and compare them with existing tau methods. For most cases, the IMW-S and IMW-P methods perform favorably.     

Developing Itô stochastic differential equation models for neuronal signal transduction pathways

Mathematical modeling and simulation of dynamic biochemical systems are receiving considerable attention due to the increasing availability of experimental knowledge of complex intracellular functions. In addition to deterministic approaches, several stochastic approaches have been developed for simulating the time-series behavior of biochemical systems. The problem with stochastic approaches, however, is the larger computational time compared to deterministic approaches. It is therefore necessary to study alternative ways to incorporate stochasticity and to seek approaches that reduce the computational time needed for simulations, yet preserve the characteristic behavior of the system in question. In this work, we develop a computational framework based on the It? stochastic differential equations for neuronal signal transduction networks. There are several different ways to incorporate stochasticity into deterministic differential equation models and to obtain It? stochastic differential equations. Two of the developed models are found most suitable for stochastic modeling of neuronal signal transduction. The best models give stable responses which means that the variances of the responses with time are not increasing and negative concentrations are avoided. We also make a comparative analysis of different kinds of stochastic approaches, that is the It? stochastic differential equations, the chemical Langevin equation, and the Gillespie stochastic simulation algorithm. Different kinds of stochastic approaches can be used to produce similar responses for the neuronal protein kinase C signal transduction pathway. The fine details of the responses vary slightly, depending on the approach and the parameter values. However, when simulating great numbers of chemical species, the Gillespie algorithm is computationally several orders of magnitude slower than the It? stochastic differential equations and the chemical Langevin equation. Furthermore, the chemical Langevin equation produces negative concentrations. The It? stochastic differential equations developed in this work are shown to overcome the problem of obtaining negative values.     

Multinomial tau-leaping method for stochastic kinetic simulations

We introduce the multinomial tau-leaping (MtauL) method for general reaction networks with multichannel reactant dependencies. The MtauL method is an extension of the binomial tau-leaping method where efficiency is improved in several ways. First, tau-leaping steps are determined simply and efficiently using a priori information and Poisson distribution-based estimates of expectation values for reaction numbers over a tentative tau-leaping step. Second, networks are partitioned into closed groups of reactions and corresponding reactants in which no group reactant set is found in any other group. Third, product formation is factored into upper-bound estimation of the number of times a particular reaction occurs. Together, these features allow larger time steps where the numbers of reactions occurring simultaneously in a multichannel manner are estimated accurately using a multinomial distribution. Furthermore, we develop a simple procedure that places a specific upper bound on the total reaction number to ensure non-negativity of species populations over a single multiple-reaction step. Using two disparate test case problems involving cellular processes--epidermal growth factor receptor signaling and a lactose operon model--we show that the tau-leaping based methods such as the MtauL algorithm can significantly reduce the number of simulation steps thus increasing the numerical efficiency over the exact stochastic simulation algorithm by orders of magnitude.     

Forward flux sampling-type schemes for simulating rare events: efficiency analysis

We analyze the efficiency of several simulation methods which we have recently proposed for calculating rate constants for rare events in stochastic dynamical systems in or out of equilibrium. We derive analytical expressions for the computational cost of using these methods and for the statistical error in the final estimate of the rate constant for a given computational cost. These expressions can be used to determine which method to use for a given problem, to optimize the choice of parameters, and to evaluate the significance of the results obtained. We apply the expressions to the two-dimensional nonequilibrium rare event problem proposed by Maier and Stein [Phys. Rev. E 48, 931 (1993)]. For this problem, our analysis gives accurate quantitative predictions for the computational efficiency of the three methods.     

Binomial tau-leap spatial stochastic simulation algorithm for applications in chemical kinetics

In cell biology, cell signaling pathway problems are often tackled with deterministic temporal models, well mixed stochastic simulators, and/or hybrid methods. But, in fact, three dimensional stochastic spatial modeling of reactions happening inside the cell is needed in order to fully understand these cell signaling pathways. This is because noise effects, low molecular concentrations, and spatial heterogeneity can all affect the cellular dynamics. However, there are ways in which important effects can be accounted without going to the extent of using highly resolved spatial simulators (such as single-particle software), hence reducing the overall computation time significantly. We present a new coarse grained modified version of the next subvolume method that allows the user to consider both diffusion and reaction events in relatively long simulation time spans as compared with the original method and other commonly used fully stochastic computational methods. Benchmarking of the simulation algorithm was performed through comparison with the next subvolume method and well mixed models (MATLAB), as well as stochastic particle reaction and transport simulations (CHEMCELL, Sandia National Laboratories). Additionally, we construct a model based on a set of chemical reactions in the epidermal growth factor receptor pathway. For this particular application and a bistable chemical system example, we analyze and outline the advantages of our presented binomial tau-leap spatial stochastic simulation algorithm, in terms of efficiency and accuracy, in scenarios of both molecular homogeneity and heterogeneity.     

Transition-state theory rate calculations with a recrossing-free moving dividing surface

Two different methods for transition-state theory (TST) rate calculations are presented that use the recently developed notions of the moving dividing surface and the associated moving separatrices: one is based on the flux-over-population approach and the other on the calculation of the reactive flux. The flux-over-population rate can be calculated in two ways by averaging the flux first over the noise and then over the initial conditions or vice versa. The former entails the calculation of reaction probabilities and is closely related to previous TST rate derivations. The latter results in an expression for the transmission factor as the noise average of a stochastic variable that is given explicitly as a function of the moving separatrices. Both the reactive-flux and flux-over-population methods suggest possible new ways of calculating approximate rates in anharmonic systems. In particular, numerical simulations of harmonic and anharmonic systems have been used to calculate reaction rates based on the reactive flux calculation using the fixed and moving dividing surfaces so as to illustrate the computational advantages of the latter.     

Efficient step size selection for the tau-leaping simulation method

The tau-leaping method of simulating the stochastic time evolution of a well-stirred chemically reacting system uses a Poisson approximation to take time steps that leap over many reaction events. Theory implies that tau leaping should be accurate so long as no propensity function changes its value "significantly" during any time step tau. Presented here is an improved procedure for estimating the largest value for tau that is consistent with this condition. This new tau-selection procedure is more accurate, easier to code, and faster to execute than the currently used procedure. The speedup in execution will be especially pronounced in systems that have many reaction channels.     

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.     

Two classes of quasi-steady-state model reductions for stochastic kinetics

The quasi-steady-state approximation (QSSA) is a model reduction technique used to remove highly reactive species from deterministic models of reaction mechanisms. In many reaction networks the highly reactive intermediates (QSSA species) have populations small enough to require a stochastic representation. In this work we apply singular perturbation analysis to remove the QSSA species from the chemical master equation for two classes of problems. The first class occurs in reaction networks where all the species have small populations and the QSSA species sample zero the majority of the time. The perturbation analysis provides a reduced master equation in which the highly reactive species can sample only zero, and are effectively removed from the model. The reduced master equation can be sampled with the Gillespie algorithm. This first stochastic QSSA reduction is applied to several example reaction mechanisms (including Michaelis-Menten kinetics) [Biochem. Z. 49, 333 (1913)]. A general framework for applying the first QSSA reduction technique to new reaction mechanisms is derived. The second class of QSSA model reductions is derived for reaction networks where non-QSSA species have large populations and QSSA species numbers are small and stochastic. We derive this second QSSA reduction from a combination of singular perturbation analysis and the Omega expansion. In some cases the reduced mechanisms and reaction rates from these two stochastic QSSA models and the classical deterministic QSSA reduction are equivalent; however, this is not usually the case.