#### DMCA

## Stochastic simulation of chemical kinetics

### Cached

### Download Links

Venue: | Annu. Rev. Phys. Chem |

Citations: | 200 - 0 self |

### Citations

3105 | Numerical Recipes: The Art of Scientific Computing - Press, Teukolsky, et al. - 2007 |

427 | Efficient exact stochastic simulation of chemical systems with many species and many channels
- Gibson, Bruck
(Show Context)
Citation Context ...es. It can be proved (L. Lok, unpublished manuscript) that this procedure generates values for τ and j = (l, j ) in exact accord with the joint density function (Equation 8). It reduces to the direct method if all the reactions are taken to be members of one family, and it reduces to the first-reaction method if each reaction is taken to be a family unto itself. For intermediate partitionings, the method may afford bookkeeping advantages when M is large. A reformulation of the SSA that, for large N and M, significantly increases its speed as compared with the direct method is Gibson & Bruck’s (20) next-reaction method. 42 Gillespie A nn u. R ev . P hy s. C he m . 2 00 7. 58 :3 5- 55 . D ow nl oa de d fr om w w w .a nn ua lr ev ie w s. or g by U ni ve rs ity o f E di nb ur gh o n 01 /0 6/ 11 . F or p er so na l u se o nl y. ANRV308-PC58-02 ARI 21 February 2007 11:10 Essentially a heavily revised version of the first-reaction method, the next-reaction method saves the putative next firing times of all reaction channels in an indexed binary tree priority queue, which is constructed so that the firing time of each parent node is always earlier than the firing times of its two daughter node... |

224 |
Exact stochastic simulation of coupled chemical reactions.
- DT
- 1977
(Show Context)
Citation Context ...ven X(t) = x, that the next reaction in the system will occur in the infinitesimal time interval [t + τ, t + τ + dτ ), and will be an Rj reaction. (7) Formally, this function is the joint probability density function of the two random variables time to the next reaction (τ ) and index of the next reaction ( j ), given that the system is currently in state x. It is not difficult to derive an exact formula for p(τ, j |x, t) by applying the laws of probability to the fundamental premise (Equation 2). The result is (8–11) p(τ, j |x, t) = aj(x) exp(−a0(x) τ ), (8) where a0(x) = M∑ j ′=1 aj ′ (x). (9) Equation 8 is the mathematical basis for the stochastic simulation approach. It implies that τ is an exponential random variable with mean (and standard deviation) 1/a0(x), while j is a statistically independent integer random variable with point probabilities aj(x)/a0(x). There are several exact Monte Carlo procedures for generating samples of τ and j according to these distributions. Perhaps the simplest is the so-called direct method, which follows by applying the standard inversion generating method of Monte Carlo theory (11): We draw two random numbers r1 and r2 from the uniform www.annu... |

179 | Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics - Haseltine, Rawlings |

146 |
Stochastic chemical kinetics and the quasi-steady-state assumption: Application to the Gillespie algorithm
- Rao, Arkin
(Show Context)
Citation Context ... the last two terms in the brackets. Tests suggest that the trapezoidal formula often gives more accurate results than Equation 23 for a given value of τ . 6.2. The Slow-Scale Stochastic Simulation Algorithm A different approach to stochastically simulating stiff chemical systems is one that was inspired by the well-known Michaelis-Menten approximation in deterministic chemical kinetics (39). Several different ways of realizing this approach have been proposed (40–48), but all are rooted in the same basic idea of eliminating the fast entities through a kind of quasi-steady-state approximation (41). Arguably the clearest articulation of this approach, at least as regards its theoretical justification within stochastic chemical kinetics, is the ssSSA of References 42 and 45. The ssSSA proceeds in a series of steps, the first of which is to make a provisional partitioning of the reaction channels R = {R1, . . . , RM} into fast and slow subsets, Rf and Rs. Assigned to Rf are those reactions whose propensity functions tend to have the largest values. All the other reactions are assigned to Rs. If it is not obvious how to make this partitioning, then it may be that the system is not really s... |

124 | The slow-scale stochastic simulation algorithms
- Cao, Gillespie, et al.
(Show Context)
Citation Context ... A nn u. R ev . P hy s. C he m . 2 00 7. 58 :3 5- 55 . D ow nl oa de d fr om w w w .a nn ua lr ev ie w s. or g by U ni ve rs ity o f E di nb ur gh o n 01 /0 6/ 11 . F or p er so na l u se o nl y. ANRV308-PC58-02 ARI 21 February 2007 11:10 reactions is deemed acceptable; otherwise we must try another set of fast reactions, or else we must conclude that the system is not really stiff, and the ssSSA cannot be applied. With the stochastic stiffness conditions satisfied, we now invoke the slow-scale approximation—a result that can be mathematically derived from the fundamental premise (Equation 2) (42). The slow-scale approximation states, in essence, that we can ignore the fast reactions and simulate the system one slow reaction at a time, provided we replace the propensity function of each slow reaction by its average with respect to the asymptotic virtual fast process Xf(∞). More precisely, if P (yf, ∞ |xf, xs) is the probability that Xf(∞) = yf given that X(t) = (xf, xs), then the propensity function a sj (x f, xs) of each slow reaction Rsj at time t can be approximated on the timescale of the slow reactions by a sj (x f, xs) ≡ ∑ yf P (yf, ∞ |xf, xs) a sj (yf, xs). (24) The ssSSA t... |

115 |
Stochastic approach to chemical kinetics,
- McQuarrie
- 1967
(Show Context)
Citation Context ...uation 6). The catch is that the SSA is often very slow, essentially because it insists on simulating every individual reaction event. The mathematical reason for this slowness can be traced to the factor 1/a0(x) in Equation 10a, which will be small if any reactant population is large. To illustrate the foregoing ideas, consider the simple reaction S c−−→ Ø. (11) The propensity function for this reaction is a(x) = c x, and the state-change vector is ν = −1. The RRE (Equation 6) reads d X/dt = −c X, and the solution to this ODE for the initial condition X(0) = x0 is X(t) = x0e−ct(RREsolution). (12) The CME (Equation 4) reads ∂ dt P (x, t |x0, 0) = a(x + 1)P (x + 1, t |x0, 0) − a(x)P (x, t |x0, 0). 40 Gillespie A nn u. R ev . P hy s. C he m . 2 00 7. 58 :3 5- 55 . D ow nl oa de d fr om w w w .a nn ua lr ev ie w s. or g by U ni ve rs ity o f E di nb ur gh o n 01 /0 6/ 11 . F or p er so na l u se o nl y. ANRV308-PC58-02 ARI 21 February 2007 11:10 X(t) 100 80 S Ø 60 40 20 0 0 1 2 3 4 5 6 t c = 1, X (0) = 100 Figure 1 Simulating the simple isomerization reaction (Equation 11). The thin light blue line shows the solution (Equation 12) of the reaction-rate equation (RRE). The two dashed gray l... |

105 |
Improved leap-size selection for accelerated stochastic simulation
- Gillespie, Petzold
- 2003
(Show Context)
Citation Context ...n the exact SSA. But several practical issues need to be resolved to effectively implement this strategy: First, how can we estimate in advance the largest value of τ that satisfies the leap condition? Second, how can we ensure that the generated kj-values do not cause some Rj to fire so many times that the population of some reactant is driven negative? Finally, how can we arrange it so that tau-leaping segues efficiently to the SSA? The method originally suggested in Reference 26 for estimating the largest value of τ that satisfies the leap condition has undergone two successive refinements (31, 32). The latest τ -selection procedure (32) is not only more accurate than the earlier procedures, but also faster, especially if M is large. It computes the largest value of τ for which the estimated fractional change τ aj/aj in each propensity function during τ is bounded by a user-specified accuracy-control parameter ε (0 < ε 1). However, it does this in an indirect way: It chooses τ so that the estimated fractional change τ xi/xi in each reactant population is bounded by an amount εi = εi(ε, xi) (except that no xi is required to change by an amount less than one), where the functions 46 Gil... |

104 | Stiffness in stochastic chemically reacting systems: The implicit tau-leaping method
- Rathinam, Petzold, et al.
- 2003
(Show Context)
Citation Context ...ort is usually more than compensated by the ability to use much-larger values of t for the same degree of accuracy. The tau-leaping formula (Equation 17), where x = X(t), is obviously an explicit updating formula. To make it implicit by replacing x in the argument of the Poisson random variable with X(t+τ ), however, raises some serious questions in the context of Markov process theory, where updates are supposed to be past-forgetting; moreover, even if that replacement could be justified theoretically, there appears to be no way to solve the resulting equation for X(t + τ ). Rathinam et al. (37) have proposed a partial implicitization, in the following implicit tau-leaping formula: X(t + τ ) .= x + M∑ j = 1 [Pj(aj(x)τ ) − aj(x)τ + aj(X(t + τ ))τ ]νj. (23) In this formula, the mean of the Poisson random variable Pj(aj(x)τ ) is subtracted out and replaced by its value at the later time t + τ , but the variance has been left unchanged. The advantage of this formula is that, once the Poisson random numbers have been generated using the current state x, the equation can then be solved for X(t + τ ) using the same (deterministic) numerical techniques developed for implicit ODE solvers (36)... |

102 |
Efficient formulation of the stochastic simulation algorithm for chemically reacting systems
- Cao, Li, et al.
- 2004
(Show Context)
Citation Context ...eactions and species at the beginning of a simulation, and for the yeast system in particular, a just-in-time introduction of the reactions and species makes feasible an otherwise infeasible simulation. Lok also observes that this just-in-time strategy cannot be used with the RRE, which is noteworthy because the RRE too becomes unwieldy when enormous numbers of species and reaction channels are involved; thus we have yet another reason for using stochastic simulation instead of deterministic simulation on biological systems. For further discussion of these points, see Reference 22. Cao et al. (21) recently introduced a modified direct method that often makes the direct method competitive in speed with the next-reaction method. These authors observed that if the reaction channels are indexed so that reactions Rj with larger propensity functions are assigned lower index values j , then the average number of terms summed in the computation (Equation 10b) is minimized. The consequent gain in speed can be significant for systems with many reactions and a wide range of propensity function values, a common circumstance in biological models. The modified direct method starts off with a relativ... |

75 |
Binomial leap methods for simulating stochastic chemical kinetics.
- Tian, Burrage
- 2004
(Show Context)
Citation Context ... ev . P hy s. C he m . 2 00 7. 58 :3 5- 55 . D ow nl oa de d fr om w w w .a nn ua lr ev ie w s. or g by U ni ve rs ity o f E di nb ur gh o n 01 /0 6/ 11 . F or p er so na l u se o nl y. ANRV308-PC58-02 ARI 21 February 2007 11:10 Condition 21 is deemed to be adequately fulfilled if it is satisfied by both 〈τ xi〉 and√ var{τ xi}. The resulting set of inequalities yields an efficient, explicit formula for the largest permissible value of τ (32). As for keeping the generated random numbers kj from driving the Rj reactant populations negative, several strategies have been proposed. Tian & Burrage (33)— and, independently, Chatterjee et al. (34)—proposed approximating the unbounded Poisson random numbers kj with bounded binomial random numbers. But it turns out that it is usually not the unboundedness of the Poisson kj’s that produces negative populations, but rather the lack of coordination in tau-leaping between different reaction channels that separately decrease the population of a common species. Cao et al. (35) have proposed a different approach that resolves this difficulty and also establishes a smooth connection with the SSA. In their approach, we first identify as critical all tho... |

65 | Multiscale stochastic simulation algorithm with stochastic partial equilibrium assumption for chemically reacting systems. - Cao, Gillespie, et al. - 2005 |

59 |
Automatic generation of cellular reaction networks with Moleculizer 1.0. Nat Biotechnol
- Lok, Brent
- 2005
(Show Context)
Citation Context ...the next occurring reaction are therefore always available at the top node of the queue. The indexing scheme and the binary-tree structure of the queue facilitate updating the queue as a result of changes caused by occurring reactions. With added effort, we can even make the next-reaction method consume only one uniform random number per additional reaction event, in contrast to the two required by the direct method. Although the next-reaction method can be significantly faster than the direct method, it is much more challenging to code. For more details, see References 20 and 21. Lok & Brent (22) have developed a stochastic simulation software package called Moleculizer that uses a slightly simplified version of the next-reaction method, but with a unique twist: Reaction channels and species are introduced only when they are needed, and removed when they are not needed. The target application for Moleculizer is the simulation of the pheromone-response mechanism in yeast, a chemical system that entails a potentially enormous number of species and reaction channels. Lok showed that, at least when using the SSA, it is not necessary to introduce all the reactions and species at the beginn... |

53 |
Avoiding negative populations in explicit Poisson tau-leaping.
- Cao, DT, et al.
- 2005
(Show Context)
Citation Context ...issible value of τ (32). As for keeping the generated random numbers kj from driving the Rj reactant populations negative, several strategies have been proposed. Tian & Burrage (33)— and, independently, Chatterjee et al. (34)—proposed approximating the unbounded Poisson random numbers kj with bounded binomial random numbers. But it turns out that it is usually not the unboundedness of the Poisson kj’s that produces negative populations, but rather the lack of coordination in tau-leaping between different reaction channels that separately decrease the population of a common species. Cao et al. (35) have proposed a different approach that resolves this difficulty and also establishes a smooth connection with the SSA. In their approach, we first identify as critical all those reactions with nonzero propensities that are currently within nc firings of exhausting one of its reactants, nc being a user-specified integer. All other reactions are called noncritical. The noncritical reactions are handled by the regular Poisson tau-leaping method, and a maximum leap time τ ′ is computed for them using the procedure described in the preceding paragraph. For the critical reactions, we use the direc... |

41 |
Approximate accelerated stochastic simulation of chemically reacting systems.
- DT
- 2001
(Show Context)
Citation Context ...dom variable with mean (and variance) aj(x) τ . Therefore, to the degree that the leap condition is 44 Gillespie A nn u. R ev . P hy s. C he m . 2 00 7. 58 :3 5- 55 . D ow nl oa de d fr om w w w .a nn ua lr ev ie w s. or g by U ni ve rs ity o f E di nb ur gh o n 01 /0 6/ 11 . F or p er so na l u se o nl y. ANRV308-PC58-02 ARI 21 February 2007 11:10 Chemical Langevin equation (CLE): a differential equation driven by zero-mean Gaussian noise that describes tau-leaping when the reactant populations are sufficiently large satisfied, we can approximately leap the system ahead by a time τ by taking (25, 26) X(t + τ ) .= x + M∑ j = 1 Pj(aj(x)τ )νj, (17) where x = X(t), and Pj(mj) is a statistically independent Poisson random variable with mean (and variance) mj. Equation 17 is the basic tau-leaping formula. In the next section below I discuss how we can use it to perform faster stochastic simulations. But for now, let us suppose that τ is not only small enough to satisfy the leap condition, but also large enough that the expected number of firings of each reaction channel Rj during τ is 1: aj(x) τ 1 for all j = 1, . . . , M. (18) Then, denoting the normal (Gaussian) random variable with mean ... |

33 |
Accelerated stochastic simulation of the stiff enzyme-substrate reaction.
- Cao, DT, et al.
- 2005
(Show Context)
Citation Context ...bit the populations of the fast species whenever desired by Monte Carlo sampling the probability function P . The approaches of References 43, 44, 47 and 48 differ from the ssSSA approach described above in the way the averages (Equation 24) are computed and the way in which the fast-species populations are generated. All rely on making relatively short SSA runs of the virtual fast process between the slow reactions. Although the ssSSA can be challenging to implement, it has been successfully applied to a number of simple stiff systems (42), as well as the prototypical MichaelisMenten system (45) that is so ubiquitous in enzymatic reactions. These applications showed increases in simulation speed over the exact SSA of two to three orders of magnitude with no perceptible loss of simulation accuracy. 7. OUTLOOK The robustness and efficiency of both the SSA and the explicit tau-leaping algorithm have been considerably improved in the past few years, and those procedures seem to be nearing maturity. However, there is still room for improvement on the stiffness problem. Such improvement may come in the form of refinements to the implicit tau-leaping procedure and the ssSSA, and a clarifica... |

19 | Overcoming stiffness in stochastic simulation stemming from partial equilibrium: a multiscale Monte Carlo algorithm,” - Vlachos - 2005 |

17 | Nested stochastic simulation algorithm for chemical kinetic systems with disparate rates. - Weinan, Liu, et al. - 2005 |

16 |
Binomial distribution based τ -leap accelerated stochastic simulation.
- Chatterjee, Vlachos, et al.
- 2005
(Show Context)
Citation Context ... D ow nl oa de d fr om w w w .a nn ua lr ev ie w s. or g by U ni ve rs ity o f E di nb ur gh o n 01 /0 6/ 11 . F or p er so na l u se o nl y. ANRV308-PC58-02 ARI 21 February 2007 11:10 Condition 21 is deemed to be adequately fulfilled if it is satisfied by both 〈τ xi〉 and√ var{τ xi}. The resulting set of inequalities yields an efficient, explicit formula for the largest permissible value of τ (32). As for keeping the generated random numbers kj from driving the Rj reactant populations negative, several strategies have been proposed. Tian & Burrage (33)— and, independently, Chatterjee et al. (34)—proposed approximating the unbounded Poisson random numbers kj with bounded binomial random numbers. But it turns out that it is usually not the unboundedness of the Poisson kj’s that produces negative populations, but rather the lack of coordination in tau-leaping between different reaction channels that separately decrease the population of a common species. Cao et al. (35) have proposed a different approach that resolves this difficulty and also establishes a smooth connection with the SSA. In their approach, we first identify as critical all those reactions with nonzero propensities that ... |

15 |
Trapezoidal tau-leaping formula for the stochastic simulation of biochemical systems
- Cao, Petzold
(Show Context)
Citation Context ...tochastic Simulation of Chemical Kinetics 49 A nn u. R ev . P hy s. C he m . 2 00 7. 58 :3 5- 55 . D ow nl oa de d fr om w w w .a nn ua lr ev ie w s. or g by U ni ve rs ity o f E di nb ur gh o n 01 /0 6/ 11 . F or p er so na l u se o nl y. ANRV308-PC58-02 ARI 21 February 2007 11:10 simulations made with the exact SSA. But Equation 23 suppresses these fluctuations. However, we can restore the properly fluctuating fast variables whenever desired by taking a succession of much shorter explicit tau-leaps or SSA steps, a tactic called downshifting. For more details, see Reference 37. Cao & Petzold (38) subsequently proposed a trapezoidal implicit tau-leaping formula, which has the same form as Equation 23 except that a factor of 1/2 appears in front of each of the last two terms in the brackets. Tests suggest that the trapezoidal formula often gives more accurate results than Equation 23 for a given value of τ . 6.2. The Slow-Scale Stochastic Simulation Algorithm A different approach to stochastically simulating stiff chemical systems is one that was inspired by the well-known Michaelis-Menten approximation in deterministic chemical kinetics (39). Several different ways of realizing this ap... |

13 |
Markov Processes: An Introduction for Physical Scientists.
- DT
- 1992
(Show Context)
Citation Context ...t) = aj(x) exp(−a0(x) τ ), (8) where a0(x) = M∑ j ′=1 aj ′ (x). (9) Equation 8 is the mathematical basis for the stochastic simulation approach. It implies that τ is an exponential random variable with mean (and standard deviation) 1/a0(x), while j is a statistically independent integer random variable with point probabilities aj(x)/a0(x). There are several exact Monte Carlo procedures for generating samples of τ and j according to these distributions. Perhaps the simplest is the so-called direct method, which follows by applying the standard inversion generating method of Monte Carlo theory (11): We draw two random numbers r1 and r2 from the uniform www.annualreviews.org • Stochastic Simulation of Chemical Kinetics 39 A nn u. R ev . P hy s. C he m . 2 00 7. 58 :3 5- 55 . D ow nl oa de d fr om w w w .a nn ua lr ev ie w s. or g by U ni ve rs ity o f E di nb ur gh o n 01 /0 6/ 11 . F or p er so na l u se o nl y. ANRV308-PC58-02 ARI 21 February 2007 11:10 Stochastic simulation algorithm (SSA): a Monte Carlo procedure for numerically generating time trajectories of the molecular populations in exact accordance with the CME distribution in the unit interval, and take τ = 1 a0(x) ln ( 1 r1 ... |

13 |
Efficient stepsize selection for the tau-leaping simulation method”.
- Cao, Gillespie, et al.
- 2006
(Show Context)
Citation Context ...n the exact SSA. But several practical issues need to be resolved to effectively implement this strategy: First, how can we estimate in advance the largest value of τ that satisfies the leap condition? Second, how can we ensure that the generated kj-values do not cause some Rj to fire so many times that the population of some reactant is driven negative? Finally, how can we arrange it so that tau-leaping segues efficiently to the SSA? The method originally suggested in Reference 26 for estimating the largest value of τ that satisfies the leap condition has undergone two successive refinements (31, 32). The latest τ -selection procedure (32) is not only more accurate than the earlier procedures, but also faster, especially if M is large. It computes the largest value of τ for which the estimated fractional change τ aj/aj in each propensity function during τ is bounded by a user-specified accuracy-control parameter ε (0 < ε 1). However, it does this in an indirect way: It chooses τ so that the estimated fractional change τ xi/xi in each reactant population is bounded by an amount εi = εi(ε, xi) (except that no xi is required to change by an amount less than one), where the functions 46 Gil... |

11 |
Computer Methods for Ordinary Differential Equations and Differential Algebraic Equations.
- UM, LR
- 1998
(Show Context)
Citation Context ...wly, and off which the state point 48 Gillespie A nn u. R ev . P hy s. C he m . 2 00 7. 58 :3 5- 55 . D ow nl oa de d fr om w w w .a nn ua lr ev ie w s. or g by U ni ve rs ity o f E di nb ur gh o n 01 /0 6/ 11 . F or p er so na l u se o nl y. ANRV308-PC58-02 ARI 21 February 2007 11:10 Slow-scale SSA (ssSSA): an approximation to the SSA for systems with fast and slow reactions in which only the latter are explicitly simulated moves rapidly toward the slow manifold. Researchers have devoted much effort over the years to understanding and overcoming the computational problems posed by stiff ODEs (36) because such equations arise in many practical contexts. Stiff RREs in particular are quite commonplace. Stiffness is just as computationally troublesome in the stochastic context. When the SSA simulates a stiff system, it moves along as usual, one reaction at a time, oblivious to the stiffness. But because the great majority of the reactions are the usually uninteresting fast ones, the simulation proceeds slowly from a practical point of view. The explicit tau-leaping algorithm also performs as advertised on stiff systems. But because the τ -selection procedure that keeps the algorithm accur... |

9 | Equation-free probabilistic steady-state approximation: dynamic application to the stochastic simulation of biochemical reaction networks. - Salis, Kaznessis - 2005 |

5 |
Elements of the Theory of Markov Processes and Their Applications.
- AT
- 1960
(Show Context)
Citation Context ...ne shows the solution (Equation 12) of the reaction-rate equation (RRE). The two dashed gray lines show the one-standard-deviation envelope 〈X(t)〉 ± sdev {X(t)} of Equations 14a,b as predicted by the solution of the chemical master equation (CME) (Equation 13). The red and blue jagged curves show the trajectories produced by two separate runs of the stochastic simulation algorithm (SSA). Because P (x0 + 1, t |x0, 0) ≡ 0, we can solve this equation exactly by successively putting x = x0, x0 − 1, . . . , 0. The result is P (x, t |x0, 0) = x0!x!(x0 − x)! e −cxt(1 − e−ct)x0−x (x = 0, . . . , x0), (13) which we recognize as the probability density function for the binomial random variable with mean and standard deviation 〈X(t)〉 = x0e−ct, (14a) sdev{X(t)} = √ x0e−ct(1 − e−ct). (14b) Note that 〈X(t)〉 is identical to the solution (Equation 12) of the RRE; this will always be so if all the propensity functions are linear in the populations, but not otherwise. The SSA for this reaction is simple: In state x at time t, we draw a unit-interval uniform random number r , increase t by τ = (1/ax) ln(1/r), decrease x by 1, and then repeat. Figure 1 shows numerical results for c = 1 and x0 = 100. 3. EL... |

3 |
The multivariate Langevin and Fokker-Planck equations.
- DT
- 1996
(Show Context)
Citation Context ...∑ j = 1 νj √ aj(x)Nj(0, 1) √ τ , (19) where x = X(t), and each Nj(0, 1) is a statistically independent normal random variable with mean 0 and variance 1. Again, this equation is valid only to the extent that during τ , no propensity function changes its value significantly, yet every reaction channel fires many more times than once. It is usually possible to find a τ that satisfies these opposing conditions if all the reactant populations are sufficiently large. In the theory of continuous Markov processes, it can be shown that the CLE (Equation 19) can also be written in the white-noise form (25, 27) dX(t) dt .= M∑ j = 1 νj aj(X(t)) + M∑ j = 1 νj √ aj(X(t)) j(t) . (20) Here the j(t) are statistically independent Gaussian white-noise processes, satisfying 〈j(t) j ′ (t ′)〉 = δjj ′ δ(t−t ′), where the first delta function is Kronecker’s and the second is Dirac’s. Equation 20 is just another way of writing Equation 19; the two equations are mathematically equivalent. Equations of the form of Equation 20, with the right side the sum of a deterministic drift term and a stochastic diffusion term proportional to Gaussian white noise, are known as Langevin equations or stochastic differential ... |

3 | On the origins of approximations for stochastic chemical kinetics. - EL, JB - 2005 |

2 | A rigorous derivation of the chemical master equation. - DT - 1992 |

2 | Stochastic analysis of an oscillating chemical reaction. - Nakanishi - 1972 |

2 | The chemical Langevin and Fokker-Planck equations for the reversible isomerization reaction. - DT - 2002 |

2 |
Binomial distribution based % -leap accelerated stochastic simulation.
- Chatterjee, Vlachos, et al.
- 2005
(Show Context)
Citation Context ... nl oa de d fro m a rjo ur na ls. an nu al re vi ew s.o rg by U ni ve rs ite d e M on tre al o n 09 /2 6/ 07 . F or p er so na l u se o nl y. ANRV308-PC58-02 ARI 21 February 2007 11:10 Condition 21 is deemed to be adequately fulfilled if it is satisfied by both $#% xi% and) var{#% xi}. The resulting set of inequalities yields an efficient, explicit formula for the largest permissible value of % (32). As for keeping the generated random numbers kj from driving the Rj reactant populations negative, several strategies have been proposed. Tian & Burrage (33)— and, independently, Chatterjee et al. (34)—proposed approximating the unbounded Poisson random numbers kj with bounded binomial random numbers. But it turns out that it is usually not the unboundedness of the Poisson kj’s that produces negative populations, but rather the lack of coordination in tau-leaping between different reaction channels that separately decrease the population of a common species. Cao et al. (35) have proposed a different approach that resolves this difficulty and also establishes a smooth connection with the SSA. In their approach, we first identify as critical all those reactions with nonzero propensities that ... |

1 | Simulation of first-order chemical reaction as a stochastic process on a digital computer. - Solc, Horsak - 1972 |

1 | Simulation study of a stochastic model of reversible first-order reaction equilibrium. - Horsak, Solc - 1973 |

1 |
Discrete simulation methods in combustion kinetics.
- DL, Garret, et al.
- 1974
(Show Context)
Citation Context ...erefore, to the degree that the leap condition is 44 Gillespie A nn u. R ev . P hy s. C he m . 2 00 7. 58 :3 5- 55 . D ow nl oa de d fr om w w w .a nn ua lr ev ie w s. or g by U ni ve rs ity o f E di nb ur gh o n 01 /0 6/ 11 . F or p er so na l u se o nl y. ANRV308-PC58-02 ARI 21 February 2007 11:10 Chemical Langevin equation (CLE): a differential equation driven by zero-mean Gaussian noise that describes tau-leaping when the reactant populations are sufficiently large satisfied, we can approximately leap the system ahead by a time τ by taking (25, 26) X(t + τ ) .= x + M∑ j = 1 Pj(aj(x)τ )νj, (17) where x = X(t), and Pj(mj) is a statistically independent Poisson random variable with mean (and variance) mj. Equation 17 is the basic tau-leaping formula. In the next section below I discuss how we can use it to perform faster stochastic simulations. But for now, let us suppose that τ is not only small enough to satisfy the leap condition, but also large enough that the expected number of firings of each reaction channel Rj during τ is 1: aj(x) τ 1 for all j = 1, . . . , M. (18) Then, denoting the normal (Gaussian) random variable with mean m and variance σ 2 by N (m, σ 2), and invoking... |

1 |
Simulation study of a stochastic models of higher order reactions.
- Horsak, Solc
- 1975
(Show Context)
Citation Context ...n approximately leap the system ahead by a time τ by taking (25, 26) X(t + τ ) .= x + M∑ j = 1 Pj(aj(x)τ )νj, (17) where x = X(t), and Pj(mj) is a statistically independent Poisson random variable with mean (and variance) mj. Equation 17 is the basic tau-leaping formula. In the next section below I discuss how we can use it to perform faster stochastic simulations. But for now, let us suppose that τ is not only small enough to satisfy the leap condition, but also large enough that the expected number of firings of each reaction channel Rj during τ is 1: aj(x) τ 1 for all j = 1, . . . , M. (18) Then, denoting the normal (Gaussian) random variable with mean m and variance σ 2 by N (m, σ 2), and invoking the mathematical fact that a Poisson random variable with a mean and variance that is 1 can be approximated as a normal random variable with that same mean and variance, we can further approximate Equation 17 as X(t + τ ) .= x + M∑ j = 1 Nj(aj(x)τ, aj(x)τ ) νj = x + M∑ j = 1 [ aj(x)τ + √ aj(x)τNj(0, 1) ] νj. The last step here invokes the well-known property of the normal random variable that N (m, σ 2) = m + σN (0, 1). Collecting terms gives us what is known as the chemical Langevin... |

1 |
Time dependent fluctuation in a nonlinear chemical system.
- Nakanishi
- 1976
(Show Context)
Citation Context ...atical fact that a Poisson random variable with a mean and variance that is 1 can be approximated as a normal random variable with that same mean and variance, we can further approximate Equation 17 as X(t + τ ) .= x + M∑ j = 1 Nj(aj(x)τ, aj(x)τ ) νj = x + M∑ j = 1 [ aj(x)τ + √ aj(x)τNj(0, 1) ] νj. The last step here invokes the well-known property of the normal random variable that N (m, σ 2) = m + σN (0, 1). Collecting terms gives us what is known as the chemical Langevin equation (CLE) or Langevin leaping formula (25), X(t + τ ) .= x + M∑ j = 1 νjaj(x)τ + M∑ j = 1 νj √ aj(x)Nj(0, 1) √ τ , (19) where x = X(t), and each Nj(0, 1) is a statistically independent normal random variable with mean 0 and variance 1. Again, this equation is valid only to the extent that during τ , no propensity function changes its value significantly, yet every reaction channel fires many more times than once. It is usually possible to find a τ that satisfies these opposing conditions if all the reactant populations are sufficiently large. In the theory of continuous Markov processes, it can be shown that the CLE (Equation 19) can also be written in the white-noise form (25, 27) dX(t) dt .= M∑ j = 1 νj aj(X... |

1 |
The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior.
- JM, GD, et al.
- 2006
(Show Context)
Citation Context ...ned lower index values j , then the average number of terms summed in the computation (Equation 10b) is minimized. The consequent gain in speed can be significant for systems with many reactions and a wide range of propensity function values, a common circumstance in biological models. The modified direct method starts off with a relatively short prerun using the direct method in which the average sizes of the propensity functions are assessed. Then it reindexes the reactions accordingly and resumes the simulation, at a usually much-greater speed. See Reference 21 for details. McCollum et al. (23) have recently proposed a further improvement on the direct method in what they call the sorting direct method. Similar to the modified direct method, the sorting direct method seeks to index the reaction channels in order of decreasing values of their propensity functions so as to optimize the search in Equation 10b. But the sorting method does this dynamically, and without the need for a prerun, by using a simple bubble-up tactic: Whenever a reaction channel fires, the index of the firing channel is interchanged with the index of the next lower indexed channel www.annualreviews.org • Stochas... |

1 |
Stochastic chemical kinetics.
- DT
- 2005
(Show Context)
Citation Context ...rm on the left side of the CLE (Equation 20) and the first term on the right side both grow like the system size, whereas the second term on the right grows more slowly as the square root of the system size. In the full thermodynamic limit, the last term becomes negligibly small compared with the other terms, and the CLE (Equation 20) reduces to the RRE (Equation 6). Thus we have derived the RRE from the fundamental stochastic premise (Equation 2). The approximations made in this derivation are schematized in Figure 2, which summarizes the theoretical structure of stochastic chemical kinetics (29). 5. THE EXPLICIT TAU-LEAPING SIMULATION ALGORITHM The basic tau-leaping formula (Equation 17) suggests an obvious strategy for approximately doing stochastic simulations: In the current state x, we first choose a value for τ that satisfies the leap condition. Next, we generate for each j a sample kj of the Poisson random variable with mean aj(x) τ , by using, for example, the numerical procedure described in Reference 30. (Because the Poisson random numbers k1, . . . , kM are statistically independent, they could be generated simultaneously on M parallel processors, which would result in a su... |

1 |
Enzyme Kinetics.
- KM
- 1971
(Show Context)
Citation Context ...ore details, see Reference 37. Cao & Petzold (38) subsequently proposed a trapezoidal implicit tau-leaping formula, which has the same form as Equation 23 except that a factor of 1/2 appears in front of each of the last two terms in the brackets. Tests suggest that the trapezoidal formula often gives more accurate results than Equation 23 for a given value of τ . 6.2. The Slow-Scale Stochastic Simulation Algorithm A different approach to stochastically simulating stiff chemical systems is one that was inspired by the well-known Michaelis-Menten approximation in deterministic chemical kinetics (39). Several different ways of realizing this approach have been proposed (40–48), but all are rooted in the same basic idea of eliminating the fast entities through a kind of quasi-steady-state approximation (41). Arguably the clearest articulation of this approach, at least as regards its theoretical justification within stochastic chemical kinetics, is the ssSSA of References 42 and 45. The ssSSA proceeds in a series of steps, the first of which is to make a provisional partitioning of the reaction channels R = {R1, . . . , RM} into fast and slow subsets, Rf and Rs. Assigned to Rf are those re... |