NOTICE: Scitation Maintenance Sunday, March 1, 2015.

Scitation users may experience brief connectivity issues on Sunday, March 1, 2015 between 12:00 AM and 7:00 AM EST due to planned network maintenance.

Thank you for your patience during this process.

banner image
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.
The full text of this article is not currently available.
Perspective: Stochastic algorithms for chemical kinetics
Rent this article for
Access full text Article
1. L. Wilhelmy, Ann. Phys. Chem. 81, 413 (1850);
1.L. Wilhelmy, Ann. Phys. Chem. 81, 499 (1850).
2. M. Delbrück, J. Chem. Phys. 8, 120 (1940).
3. D. McQuarrie, J. Appl. Probab. 4, 413 (1967).
4. H. McAdams and A. Arkin, Proc. Natl. Acad. Sci. U.S.A. 94, 814 (1997).
5. A. Arkin, J. Ross, and H. McAdams, Genetics 149, 1633 (1998).
6. M. Elowitz, A. Levine, E. Siggia, and P. Swain, Science 297, 1183 (2002).
7. L. Weinberger, J. Burnett, J. Toettcher, A. Arkin, and D. Schaffer, Cell 122, 169 (2005).
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.
9. D. Gillespie, J. Comput. Phys. 22, 403 (1976).
10. D. Gillespie, J. Phys. Chem. 81, 2340 (1977).
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).
14. T. Nakanishi, J. Phys. Soc. Jpn. 32, 1313 (1972);
14.T. Nakanishi, J. Phys. Soc. Jpn. 40, 1232 (1976).
15. M. Šolc and I. Horsák, Collect. Czech. Chem. Commun. 37, 2994 (1972);
15.M. Šolc and I. Horsák, Collect. Czech. Chem. Commun. 38, 2200 (1973);
15.M. Šolc and I. Horsák, Collect. Czech. Chem. Commun. 40, 321 (1975).
16. D. Bunker, B. Garrett, T. Kleindienst, and G. Long III, Combust. Flame 23, 373 (1974).
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).
18. M. Gibson and J. Bruck, J. Phys. Chem. 104, 1876 (2000).
19. D. Gillespie, Annu. Rev. Phys. Chem. 58, 35 (2007).
20. Y. Cao, H. Li, and L. Petzold, J. Chem. Phys. 121, 4059 (2004).
21. J. McCollum, G. Peterson, C. Cox, M. Simpson, and N. Samatova, Comput. Biol. Chem. 30, 39 (2006).
22. D. Anderson, J. Chem. Phys. 127, 214107 (2007).
23. T. Kurtz, Ann. Probab. 8, 682 (1980);
23.S. Ethier and T. Kurtz, Markov Processes: Characterization and Convergence (Wiley, 1986).
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).
25. A. Slepoy, A. Thompson, and S. Plimpton, J. Chem. Phys. 128, 205101 (2008).
26. S. Mauch and M. Stalzer, IEEE/ACM Trans. Comput. Biol. Bioinf. 8, 27 (2011).
27. D. Gillespie, J. Chem. Phys. 115, 1716 (2001).
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.
29. Y. Cao, D. Gillespie, and L. Petzold, J. Chem. Phys. 124, 044109 (2006).
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. 125167, Sec. 3 gives a tutorial on tau-leaping, Secs. 4 and 5 give a tutorial on the ssSSA.
31. M. Rathinam, L. Petzold, Y. Cao, and D. Gillespie, J. Chem. Phys. 119, 12784 (2003);
31.M. Rathinam, L. Petzold, Y. Cao, and D. Gillespie, J. Chem. Phys. 121, 12169 (2004);
31.M. Rathinam, L. Petzold, Y. Cao, and D. Gillespie, Multiscale Model. Simul. 4, 867 (2005).
32. A. Auger, P. Chatelain, and P. Koumoutsakos, J. Chem. Phys. 125, 084103 (2006).
33. D. Anderson, J. Chem. Phys. 128, 054103 (2008).
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);
34.D. Gillespie, Am. J. Phys. 64, 1246 (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.
37. D. Gillespie, J. Phys. Chem. B 113, 1640 (2009).
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 ZX/Ω 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.
40. T. Kurtz, Stochastic Proc. Appl. 6, 223 (1978).
41. N. van Kampen, Adv. Chem. Phys. 34, 245 (1976);
41.Stochastic Processes in Physics and Chemistry (North-Holland, 1992).
42. E. Wallace, D. Gillespie, K. Sanft, and L. Petzold, IET Syst. Biol. 6, 102 (2012).
43. R. Grima, P. Thomas, and A. Straube, J. Chem. Phys. 135, 084103 (2011).
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.Y. Cao, D. Gillespie, and L. Petzold, J. Chem. Phys. 123, 144917 (2005).
45.A more concise tutorial presentation of the ssSSA, featuring illustrative applications to reactions (11) and the classical enzyme-substrate reactions E + SESE + 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 c2c3, 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.
47. Y. Cao, D. Gillespie, and L. Petzold, J. Comput. Phys. 206, 395 (2005).
48. L. Ferm, P. Lotstedt, and A. Hellander, J. Sci. Comput. 34, 127 (2008).
49. A. Singh and J. Hespanha, IEEE Trans. Autom. Control. 56, 414 (2011).
50. W. E, D. Liu, and E. Vanden-Eijnden, J. Chem. Phys. 123, 194107 (2005);
50.D. Gillespie, L. Petzold, and Y. Cao, J. Chem. Phys. 126, 137101 (2007).
51. K. Sanft, Ph.D. thesis, University of California, Santa Barbara, 2012.
52. H. Kuwahara and I. Mura, J. Chem. Phys. 129, 165101 (2008).
53. D. Gillespie, M. Roh, and L. Petzold, J. Chem. Phys. 130, 174103 (2009).
54. M. Roh, D. Gillespie, and L. Petzold, J. Chem. Phys. 133, 174106 (2010).
55. B. Daigle, Jr., M. Roh, D. Gillespie, and L. Petzold, J. Chem. Phys. 134, 044110 (2011).
56. M. Roh, B. Daigle, Jr., D. Gillespie, and L. Petzold, J. Chem. Phys. 135, 234108 (2011).
57. R. Rubinstein and D. Kroese, The Cross-Entropy Method (Springer, 2004).
58. M. Rathinam, P. Sheppard, and M. Khammash, J. Chem. Phys. 132, 034103 (2010).
59. D. Anderson, SIAM J. Numer. Anal. 50, 2237 (2012).
60. K. Takahashi, S. Tănase-Nicola, and P. ten Wolde, Proc. Natl. Acad. Sci. U.S.A. 107, 2473 (2010).
61. K. Doubrovinski and M. Howard, Proc. Natl. Acad. Sci. U.S.A. 102, 9808 (2005).
62. J. Elf and M. Ehrenberg, Syst. Biol. 1, 230 (2004).
63. R. Metzler, Phys. Rev. Lett. 87, 068103 (2001).
64. M. Sturrock, A. Hellander, M. Matzavinos, and M. Chaplain, J. R. Soc. Interface 10, 20120988 (2013).
65. D. Fange and J. Elf, PLoS Comput. Biol. 2(6), e80 (2006).
66. C. Gardiner, K. McNeil, D. Walls, and I. Matheson, J. Stat. Phys. 14, 307 (1976).
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).
69. T. Kurtz, J. Appl. Probab. 7, 49 (1970).
70. S. Isaacson and C. Peskin, SIAM J. Sci. Comput. 28, 47 (2006).
71. S. Isaacson, D. McQueen, and C. Peskin, Proc. Nat. Acad. Sci. U.S.A. 108, 3815 (2011).
72. S. Engblom, L. Ferm, A. Hellander, and P. Lötstedt, SIAM J. Sci. Comput. 31, 1774 (2009).
73. B. Drawert, S. Engblom, and A. Hellander, BMC Syst. Biol. 6, 76 (2012).
74. H.-W. Kang, L. Zheng, and H. Othmer, J. Math. Biol. 65, 1017 (2012).
75. D. Rossinelli, B. Bayati, and P. Koumoutsakos, Chem. Phys. Lett. 451, 136 (2008).
76. K. Iyengar, L. Harris, and P. Clancy, J. Chem. Phys. 132, 094101 (2010).
77. T. Marquez-Lago and K. Burrage, J. Chem. Phys. 127, 104101 (2007).
78. S. Lampoudi, D. Gillespie, and L. Petzold, J. Chem. Phys. 130, 094104 (2009).
79. W. Koh and K. Blackwell, J. Chem. Phys. 134, 154103 (2011).
80. B. Drawert, M. Lawson, L. Petzold, and M. Khammash, J. Chem. Phys. 132, 074101 (2010).
81. L. Ferm, A. Hellander, and P. Lötstedt, J. Comput. Phys. 229, 343 (2010).
82. D. Gillespie, S. Lampoudi, and L. Petzold, J. Chem. Phys. 126, 034302 (2007).
83. S. Lampoudi, D. Gillespie, and L. Petzold, J. Comput. Phys. 228, 3656 (2009).
84. R. Grima, J. Chem. Phys. 132, 185102 (2010).
85. J. Zon, S. van Zon, and P. Rein ten Wolde, Phys. Rev. Lett. 94, 128103 (2005).
86. R. Erban and J. Chapman, Phys. Biol. 6, 046001 (2009).
87. D. Fange, O. G. Berg, P. Sjöberg, and J. Elf, Proc. Natl. Acad. Sci. U.S.A. 107, 19820 (2010).
88. S. Isaacson, SIAM J. Appl. Math. 70, 77 (2009).
89. S. Isaacson and D. Isaacson, Phys. Rev. E 80, 066106 (2009).
90. S. Hellander, A. Hellander, and L. Petzold, Phys. Rev. E 85, 042901 (2012).
91. S. Isaacson, “A convergent reaction-diffusion master equation,” preprint arXiv:1211.6772v1 (2012).
92. M. Doi, J. Phys. A: Math. Gen. 9, 1465 (1976);
92.M. Doi, J. Phys. A: Math. Gen. 9, 1479 (1976).
93. A. Hellander, S. Hellander, and P. Lötstedt, Multiscale Model. Simul. 10, 585 (2012).
94. M. Flegg, S. Chapman, and R. Erban, J Roy. Soc. Interface 9, 859 (2012).
95. K. Sanft, S. Wu, M. Roh, J. Fu, R. Lim, and L. Petzold, Bioinformatics 27, 2457 (2011).

Data & Media loading...


Article metrics loading...



We outline our perspective on stochastic chemical kinetics, paying particular attention to numerical simulation algorithms. We first focus on dilute, well-mixed systems, whose description using ordinary differential equations has served as the basis for traditional chemical kinetics for the past 150 years. For such systems, we review the physical and mathematical rationale for a discrete-stochastic approach, and for the approximations that need to be made in order to regain the traditional continuous-deterministic description. We next take note of some of the more promising strategies for dealing stochastically with stiff systems, rare events, and sensitivity analysis. Finally, we review some recent efforts to adapt and extend the discrete-stochastic approach to systems that are not well-mixed. In that currently developing area, we focus mainly on the strategy of subdividing the system into well-mixed subvolumes, and then simulating diffusional transfers of reactant molecules between adjacent subvolumes together with chemical reactions inside the subvolumes.


Full text loading...

This is a required field
Please enter a valid email address
752b84549af89a08dbdd7fdb8b9568b5 journal.articlezxybnytfddd
Scitation: Perspective: Stochastic algorithms for chemical kinetics