^{1,a)}, Andreas Hellander

^{2,b)}and Linda R. Petzold

^{2,c)}

### 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.

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.

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

### Key Topics

- Diffusion
- 26.0
- Chemical reactions
- 21.0
- Poisson's equation
- 18.0
- Chemical kinetics
- 14.0
- Langevin equation
- 9.0

_{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}.

_{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.

_{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).

_{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.

_{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.

Data & Media loading...

Article metrics loading...

### 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...

Commenting has been disabled for this content