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.
f
Perspective: Stochastic algorithms for chemical kinetics
Rent:
Rent this article for
Abstract
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.
© 2013 AIP Publishing LLC
Received 05 January 2013
Accepted 25 March 2013
Published online 01 May 2013
Acknowledgments: The work of D.T.G. was funded by the University of California, Santa Barbara under professional services Agreement No. 130401A40, pursuant to National Institutes of Health (NIH) Award No. R01-EB014877-01. The work of A.H. and L.R.P. was funded by National Science Foundation (NSF) Award No. DMS-1001012, ICB Award No. W911NF-09-0001 from the U.S. Army Research Office, NIBIB of the NIH under Award No. R01-EB014877-01, and (U.S.) Department of Energy (DOE) Award No. DE-SC0008975. The content of this paper is solely the responsibility of the authors and does not necessarily represent the official views of these agencies.
Article outline:
I. INTRODUCTION
II. DILUTE WELL-MIXED CHEMICAL SYSTEMS
A. The chemical master equation and the propensity function
B. Physical justification for the propensity function
C. The stochastic simulation algorithm
D. Tau-leaping
E. Connection to the traditional ODE approach
F. Stiff systems and the slow-scale SSA
G. Rare events
H. Sensitivity analysis
III. BEYOND WELL-MIXED SYSTEMS
A. The reaction-diffusion master equation and simulation algorithm
B. Algorithms for spatial stochastic simulation
C. The RDME on small length scales
IV. ACCOMPLISHMENTS AND CHALLENGES
/content/aip/journal/jcp/138/17/10.1063/1.4801941
1.
1. L. Wilhelmy, Ann. Phys. Chem. 81, 413 (1850);
1.L. Wilhelmy, Ann. Phys. Chem. 81, 499 (1850).
5.
5. A. Arkin, J. Ross, and H. McAdams, Genetics 149, 1633 (1998).
8.
8.If the observer knows the initial state x_{0} only as a random variable with some PDF Q_{0}(x_{0}), 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 Q_{0}. The essential role of an observer in the CME is also illustrated by the fact that in many cases, P(x, t|x_{0}, t_{0}) 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 t_{0}.
11.
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 S_{2} molecule sweeps out relative to the center of a randomly chosen S_{1} 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 S_{1} 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 q_{j} gives the probability that the two molecules will react in the next dt. And finally, that single-pair reaction probability summed over all x_{1}x_{2} distinct reactant pairs gives the probability defined in Eq. (2). If the two reactant molecules are of the same species, say S_{1}, then the number of distinct reactant pairs will instead be x_{1}(x_{1} − 1)/2. Determining the collision-conditioned reaction probability q_{j}, 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 R_{j} reaction will occur between the two colliding molecules if and only if their “collisional kinetic energy,” suitably defined, exceeds some threshold value E_{th}. In that case it can be shown [see for instance D. Gillespie, Physica A 188, 404 (1992) ] that q_{j} = exp ( − E_{th}/k_{B}T), which is the famous Arrhenius factor.
http://dx.doi.org/10.1016/0378-4371(92)90283-V
12.
12. D. Gillespie, J. Chem. Phys. 131, 164109 (2009).
http://dx.doi.org/10.1063/1.3253798
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).
http://dx.doi.org/10.1016/0095-8522(49)90023-9
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.
13. M. Smoluchowski, Z. Phys. Chem. 92, 129 (1917).
17.
17.The inversion Monte Carlo method exploits the theorem that, if F_{Y} is the cumulative distribution function (CDF) of the random variable Y, then F_{Y}(Y) will be the unit-interval uniform random variable U. Thus, if u is a random sample of U, then solving (inverting) F_{Y}(y) = u will yield a random sample y of Y. Writing Eq. (4) as the product of and a_{j}(x)/a_{0}(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 paper^{9} 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.
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.
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.
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.
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);
http://dx.doi.org/10.1119/1.18210
35.
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.
http://dx.doi.org/10.1063/1.481811
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.
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.
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.
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 a_{j}(X(t′))dt′. An established result in random variable theory is that the sum of K independent Poisson random variables with means m_{1}, …, m_{K} 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” Y_{j} 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.
45. Y. Cao, D. Gillespie, and L. Petzold, J. Chem. Phys. 122, 014116 (2005);
http://dx.doi.org/10.1063/1.1824902
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).
http://dx.doi.org/10.1049/iet-syb.2009.0057
46.
46.Under the condition c_{2} ≫ c_{3}, R_{1} will be a “fast” reaction regardless of the size of c_{1}. That is because the conversion of an S_{1} molecule into an S_{3} molecule necessarily takes exactly one more R_{1} reaction than R_{2} reaction. Thus, in most time intervals there will be at least as many R_{1} firings as R_{2} firings; hence, if R_{2} is fast, then R_{1} 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.
51. K. Sanft, Ph.D. thesis, University of California, Santa Barbara, 2012.
57.
57. R. Rubinstein and D. Kroese, The Cross-Entropy Method (Springer, 2004).
67.
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.
68.Equation (14) also follows by applying the “centered finite difference method” in numerical analysis to the diffusion Eq. (13).
91.
91. S. Isaacson, “
A convergent reaction-diffusion master equation,” preprint
arXiv:1211.6772v1 (
2012).
http://aip.metastore.ingenta.com/content/aip/journal/jcp/138/17/10.1063/1.4801941
Article metrics loading...
/content/aip/journal/jcp/138/17/10.1063/1.4801941
2013-05-01
2014-10-24
Abstract
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...
/deliver/fulltext/aip/journal/jcp/138/17/1.4801941.html;jsessionid=7bcfnia3aq8pp.x-aip-live-06?itemId=/content/aip/journal/jcp/138/17/10.1063/1.4801941&mimeType=html&fmt=ahah&containerItemId=content/aip/journal/jcp
Most read this month
Article
content/aip/journal/jcp
Journal
5
3
true
true
Commenting has been disabled for this content