No data available.
Please log in to see this content.
You have no subscription access to this content.
No metrics data to plot.
The attempt to load metrics for this article has failed.
The attempt to plot a graph for these metrics has failed.
Perspective: Stochastic algorithms for chemical kinetics
1. L. Wilhelmy, Ann. Phys. Chem. 81, 413 (1850);
1.L. Wilhelmy, Ann. Phys. Chem. 81, 499 (1850).
5. A. Arkin, J. Ross, and H. McAdams, Genetics 149, 1633 (1998).
8.If the observer knows the initial state x0 only as a random variable with some PDF Q0(x0), then the CME straightforwardly transforms, without change, to an equation for . But that function is not an “unconditioned” PDF P(x, t), because it is obviously conditioned on the initial distribution Q0. The essential role of an observer in the CME is also illustrated by the fact that in many cases, P(x, t|x0, t0) becomes independent of t as t → ∞ while X(t) does not; in that case, the only thing that eventually stops changing with time is best prediction of X that can be made by an observer who last observed X at time t0.
11.The derivation of Eq. (3a) that was given in Ref. 9 can be briefly summarized as follows: is the average “collision volume” that a randomly chosen S2 molecule sweeps out relative to the center of a randomly chosen S1 molecule in time dt. Dividing that collision volume by the system volume Ω gives, because the system is dilute and well-mixed, the probability that the center of the S1 molecule lies inside the collision volume, and hence the probability that the two molecules will collide in the next dt. That collision probability multiplied by qj gives the probability that the two molecules will react in the next dt. And finally, that single-pair reaction probability summed over all x1x2 distinct reactant pairs gives the probability defined in Eq. (2). If the two reactant molecules are of the same species, say S1, then the number of distinct reactant pairs will instead be x1(x1 − 1)/2. Determining the collision-conditioned reaction probability qj, which will always be some number between 0 and 1, requires additional physical reasoning. The best known example is for the model in which an Rj reaction will occur between the two colliding molecules if and only if their “collisional kinetic energy,” suitably defined, exceeds some threshold value Eth. In that case it can be shown [see for instance D. Gillespie, Physica A 188, 404 (1992) ] that qj = exp ( − Eth/kBT), which is the famous Arrhenius factor.
12. D. Gillespie, J. Chem. Phys. 131, 164109 (2009).
12.As discussed in Sec. VI of this reference, the analysis leading to the result (3b) can be regarded as a refined, corrected, and stochastically extended version of the analysis of F. Collins and G. Kimball, J. Colloid Sci. 4, 425 (1949).
12.A derivation of Eq. (3b) that is further improved over the one given in the 2009 paper can be found in Secs. 3.7 and 4.8 of D. Gillespie and E. Seitaridou, Simple Brownian Diffusion (Oxford University Press, 2012).
13. M. Smoluchowski, Z. Phys. Chem. 92, 129 (1917).
17.The inversion Monte Carlo method exploits the theorem that, if FY is the cumulative distribution function (CDF) of the random variable Y, then FY(Y) will be the unit-interval uniform random variable U. Thus, if u is a random sample of U, then solving (inverting) FY(y) = u will yield a random sample y of Y. Writing Eq. (4) as the product of and aj(x)/a0(x) reveals that τ and j are statistically independent random variables with those respective PDFs. Integrating (summing) those PDFs over their respective arguments gives their CDFs, and Eqs. (5a) and (5b) then follow from the inversion method. The generic forms of Eqs. (5a) and (5b), being such straightforward consequences of applying the inversion Monte Carlo method to the aforementioned PDFs of τ and j, were known long before Gillespie's 1976 paper;9 e.g., see page 36 of J. Hammersley and D. Handscomb, Monte Carlo Methods (Methuen, 1964). The main contributions of Gillespie's 1976 paper9 were: (i) proving from simple kinetic theory that bimolecular reactions in a dilute gas, like unimolecular reactions, are describable by propensity functions as defined in (2); and (ii) proving that propensity functions as defined in (2) imply that the “time to next reaction” and the “index of next reaction” are random variables distributed according to the joint PDF (4).
24. D. Anderson and T. Kurtz, “Continuous time Markov chain models for chemical reaction networks,” in Design and Analysis of Biomolecular Circuits: Engineering Approaches to Systems and Synthetic Biology, edited by H. Koeppl et al. (Springer, 2011).
28.An early conjecture was that the negativity problem in tau-leaping was caused by the unbounded Poisson random variables in the tau-leaping formula (7) occasionally allowing too many reaction firings in a single leap. That conjecture led to proposals to replace the Poisson random variables in Eq. (7) with binomial random variables, which are bounded. But it was subsequently determined that the principal causes of the negativity problem lay elsewhere. First was the flawed implementation of the leap condition (6) that was used in Ref. 27: the small-valued propensity functions that it failed to protect often have a reactant with a small molecular population which can easily be driven negative. Second, since the firing numbers of the individual reactions in the tau-leaping formula (7) are generated independently of each other, two or more reaction channels that decrease the population of a common species could inadvertently collude to overdraw that species. Neither of these two problems is fixed by the ad hoc substitution of binomial random variables for the Poisson random variables in Eq. (7). But both problems are effectively dealt with in the heavily revised tau-leaping procedure described in Ref. 29, which uses the theoretically appropriate Poisson random variables.
30. D. Gillespie, “Simulation methods in systems biology,” in Formal Methods for Computational Systems Biology, edited by M. Bernardo, P. Degano, and G. Zavattaro (Springer, 2008), pp. 125–167, Sec. 3 gives a tutorial on tau-leaping, Secs. 4 and 5 give a tutorial on the ssSSA.
34.The view we have taken here of dt as an independent real variable on the interval [0, ɛ), where ɛ is an arbitrarily small positive number, means that in Eq. (9b) is perfectly well defined. For an explanation of the connection between the factor in Eq. (9b) and “Gaussian white noise,” see D. Gillespie, Am. J. Phys. 64, 225 (1996);
35. D. Gillespie, J. Chem. Phys. 113, 297 (2000). The CLE (9b) can be shown to be mathematically equivalent to the Fokker-Planck equation that is obtained by first making a formal Taylor series expansion of the right side of the CME (1), which yields the so-called Kramers-Moyal equation, and then truncating that expansion after the second derivative term. However, that way of “obtaining” the CLE, which was known well before this 2000 paper, does not qualify as a derivation, because it does not make clear under what conditions the truncation will be accurate. In contrast, the derivation of the CLE given here provides a clear and testable criterion for the accuracy of its approximations, namely, the extent to which both leap conditions are satisfied. But see also the proviso in Ref. 36.
36.The approximation that was made in deriving the CLE from the tau-leaping formula (7), while accurate for likely values of those two random variables, is very inaccurate in the tails of their probability densities, even when m ≫ 1. Although both tails are “very near 0,” they differ by many orders of magnitude. Since rare events arise from the “unlikely” firing numbers under those tails, it follows that the CLE will not accurately describe the atypical behavior of a chemical system, even if both leap conditions are well satisfied.
38.The foregoing chain of reasoning can be summarized from the perspective of computational mathematics as follows: The definition (2) of the propensity function implies, for any time step τ that is small enough to satisfy the first leap condition (6), the tau-leaping formula (7). If the second leap condition (8) is also satisfied, the tau-leaping formula becomes Eq. (9a), which is the forward Euler formula for a stochastic differential equation. And that formula becomes in the thermodynamic limit, where its diffusion term will be negligibly small compared to its drift term, the forward Euler formula for an ordinary differential equation.
39.A common misconception is that, while the molecular population X is obviously a discrete variable, the molecular concentration Z ≡ X/Ω is a “continuous” variable. The error in that view becomes apparent when one realizes that simply by adopting a unit of length that gives Ω the value 1, Z becomes numerically equal to X. But even if that is not done, a sudden change in X, say from 10 to 9, will always result a discontinuous 10% decrease in Z. The molecular concentration Z is no less discrete, and no more continuous, than the molecular population X.
44.The last step in Eq. (11) can be heuristically understood as follows: The integral is essentially a sum of independent Poisson random variables, each indexed by t′ and having mean aj(X(t′))dt′. An established result in random variable theory is that the sum of K independent Poisson random variables with means m1, …, mK is a Poisson random variable with mean . Therefore, the integral in Eq. (11) is a Poisson random variable with mean . For an explanation of how that Poisson random variable can be viewed as a “unit-rate Poisson process” Yj that is “scaled” by , see Refs. 23 and 24. Note that all the integrals here can be written as finite algebraic sums, since X(t′) stays constant between successive reactions.
45. Y. Cao, D. Gillespie, and L. Petzold, J. Chem. Phys. 122, 014116 (2005);
45.A more concise tutorial presentation of the ssSSA, featuring illustrative applications to reactions (11) and the classical enzyme-substrate reactions E + S ⇌ ES → E + P, is given in Ref. 30. A more thorough examination of the relation between the ssSSA's treatment of the enzyme-substrate system and the classical Michaelis-Menten approximation is given in K. Sanft, D. Gillespie, and L. Petzold, IET Syst. Biol. 5, 58 (2011).
46.Under the condition c2 ≫ c3, R1 will be a “fast” reaction regardless of the size of c1. That is because the conversion of an S1 molecule into an S3 molecule necessarily takes exactly one more R1 reaction than R2 reaction. Thus, in most time intervals there will be at least as many R1 firings as R2 firings; hence, if R2 is fast, then R1 must be also. This illustrates the important fact that fast and slow reaction channels cannot always be identified solely on the basis of the magnitudes of their rate constants.
51. K. Sanft, Ph.D. thesis, University of California, Santa Barbara, 2012.
57. R. Rubinstein and D. Kroese, The Cross-Entropy Method (Springer, 2004).
67.See, for example, D. Gillespie and E. Seitaridou, Simple Brownian Diffusion (Oxford University Press, 2012), Chap. 5. That chapter's revised Sec. 5.6 (which can be downloaded from the book's webpage on the publisher's website) shows that below a certain value of h, no propensity function can give a physically accurate modeling of diffusive transfers of molecules between voxels.
68.Equation (14) also follows by applying the “centered finite difference method” in numerical analysis to the diffusion Eq. (13).
91. S. Isaacson
, “A convergent reaction-diffusion master equation
,” preprint arXiv:1211.6772v1
Article metrics loading...
Full text loading...
Most read this month