The chemical master equation (CME) and the mathematically equivalent stochastic simulation algorithm (SSA) assume that the reactant molecules in a chemically reacting system are “dilute” and “well-mixed” throughout the containing volume. Here we clarify what those two conditions mean, and we show why their satisfaction is necessary in order for bimolecular reactions to physically occur in the manner assumed by the CME and the SSA. We prove that these conditions are closely connected, in that a system will stay well-mixed if and only if it is dilute. We explore the implications of these validity conditions for the reaction-diffusion (or spatially inhomogeneous) extensions of the CME and the SSA to systems whose containing volumes are not necessarily well-mixed, but can be partitioned into cubical subvolumes (voxels) that are. We show that the validity conditions, together with an additional condition that is needed to ensure the physical validity of the diffusion-induced jump probability rates of molecules between voxels, require the voxel edge length to have a strictly positive lower bound. We prove that if the voxel edge length is steadily decreased in a way that respects that lower bound, the average rate at which bimolecular reactions occur in the reaction-diffusion CME and SSA will remain constant, while the average rate of diffusive transfer reactions will increase as the inverse square of the voxel edge length. We conclude that even though the reaction-diffusion CME and SSA are inherently approximate, and cannot be made exact by shrinking the voxel size to zero, they should nevertheless be useful in many practical situations.
REFERENCES
1.
Derivations of the CME and the SSA from the definitions of the propensity functions and the state-change vectors, derivations which are mathematically rigorous in that they invoke only the laws of probability, can be found in the following two publications: Secs. 2 and 3 of
D.
Gillespie
, “Stochastic chemical kinetics
” in Handbook of Materials Modeling
, edited by S.
Yip
(Springer
, 2005
), pp. 1735
–1752
;and Secs. 2.2 and 2.3 of
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
.2.
This fact arises from a general result in quantum mechanics which is widely known as “Fermi's Golden Rule” (although it is actually due to Dirac): For a wide class of energy-conserving transitions in atomic and molecular physics which are induced by a weak perturbation, the probability that the transition will occur in a time interval Δt will, to a good approximation, be linear in Δt, provided Δt is neither too small nor too large. The lower bound on Δt implied by the proviso is much smaller than what would be considered infinitesimally small on time scales typical of chemical reaction events. This seminal result also describes radioactive decay. But the quantum mechanical proof of this result is far from trivial, see, e.g., F. Mandl, Quantum Mechanics (Wiley, 1992), Chap. 9.
3.
D.
Gillespie
, J. Comput. Phys.
22
, 403
(1976
);D.
Gillespie
, J. Phys. Chem.
81
, 2340
(1977
);D.
Gillespie
, Physica A
188
, 404
(1992
).4.
F.
Collins
and G.
Kimball
, J. Colloid Sci.
4
, 425
(1949
).5.
An equivalent but logically neater derivation of Eq. (5) is given in D. Gillespie and E. Seitaridou, Simple Brownian Diffusion: An Introduction to the Standard Theories (Oxford University Press, 2012), Secs. 3.7 and 4.8. A version of Eq. (5) that appears more frequently in current literature can be obtained by multiplying the numerator and denominator of Eq. (5) by πσ12 and then defining |$k \equiv \pi \sigma _{12}^2 \bar v_{12} q_j $|: The quantity in parentheses on the right becomes, in the thermodynamic limit of infinitely large Ω and molecular populations, the “rate constant” for reaction Rj in traditional deterministic chemical kinetics. It was originally obtained in Ref. 4 by imposing on the standard diffusion equation a “radiation boundary condition” with an ad hoc parameter k. The advantage of the more recent derivation of Eq. (5) cited above is that it uses physical reasoning which makes it unnecessary to overtly postulate a radiation boundary condition, and it provides an explicit value for k. Note that k, which is sometimes called the “microscopic association rate,” is in fact just the reaction rate constant associated with the dilute gas propensity function in Eq. (3).
\begin{equation*}a_j (x_1,x_2 ) = \left( {\frac{{4\pi \sigma _{12} D_{12} k}}{{4\pi \sigma _{12} D_{12} + k}}} \right) \cdot | \Omega |^{ - 1} x_1 x_2.\end{equation*}
A rigorous derivation of the radiation boundary condition from the Kramers equation in the Langevin model of Brownian motion has been given by
D.
Bicout
, A.
Berezhkovskii
, and A.
Szabo
, J. Chem. Phys.
114
, 2293
(2001
). The 2009 derivation of Eq. (5) is based on the same underlying physics as the Bicout et al. derivation (i.e., ballistic motion on the smallest spatiotemporal scales), but it is arguably mathematically simpler and physically more transparent. 6.
7.
The overall logic in the following argument is not new; see, e.g.,
S. A.
Isaacson
, SIAM J. Appl. Math.
70
, 77
(2009
);R.
Grima
and S.
Schnell
, Essays Biochem.
45
, 41
(2008
). However, the specific conclusions we will draw from those arguments, namely, Eqs. (8) and (9) are so far as we can tell new.
[PubMed]
8.
In the ideal gas regime, we have |$4D_{12} \gg \sigma _{12} \bar v_{12} q_j $|, and the bimolecular propensity function (5) reduces to (3). From (3), it follows that the mean time to the next reaction of an S1 molecule with any one of x2 S2 the molecules in Ω is |$\tau _{\rm r} = ( {\pi \sigma _{12}^2 \bar v_{12} q_j | \Omega |^{ - 1} x_2 } )^{ - 1} $|. Since the molecules are now moving ballistically, the average time it takes the S1 molecule to explore the volume Ω is given, not by (6), but rather by |$\tau _{\rm b} = | \Omega |^{1/3} /\bar v_1 $|, where |$\bar v_1 $| is the average speed of an S1 molecule. The requirement for the reaction to occur under well-mixed conditions is thus τb ≪ τr. With the foregoing formulas for τb and τr, this becomes the requirement |$\sqrt {\pi ( {\bar v_{12} /\bar v_1 } )q_j x_2 } \cdot \sigma _{12} \ll | \Omega |^{1/3} $|. This condition can always be satisfied if the probability qj, that an S1−S2 collision will produce a reaction, is sufficiently close to 0. But if qj is not ≪1, then since |$\bar v_{12} /\bar v_1 \ge 1$|, the condition will be satisfied only if |$\sqrt {\pi x_2 } \sigma _{12} \ll | \Omega |^{1/3} $|, i.e., only if the reactant molecules are dilute inside Ω.
9.
C.
Gardiner
, K.
McNeil
, D.
Walls
, and I.
Matheson
, J. Stat. Phys.
14
, 307
(1976
). It is interesting to note that this seminal paper uses the terminology “diffusion-reaction master equation.”10.
D.
Gillespie
and E.
Seitaridou
, Simple Brownian Diffusion: An Introduction to the Standard Theories
(Oxford University Press
, 2012
), see Secs. 4.2–4.6.11.
Two different derivations of condition (14) can be found in Secs. 5.6 and 9.4 of the book cited in Ref. 10, as revised and corrected in the downloadable errata at http://ukcatalogue.oup.com/product/9780199664504.do#. These (revised) sections also make it clear that if condition (14) is not satisfied, then the movement of reactant molecules between voxels cannot be accurately described by any propensity function.
12.
See the paper by S. A. Isaacson cited in Ref. 7.
13.
The fact that dt here is an infinitesimal ensures that the occurrence of more than one reaction firing in dt will be so rare that those firings are, for all practical purposes, “mutually exclusive.” That is crucial for invoking the addition law of probability.
14.
D.
Fange
, O.
Berg
, P.
Sjöberg
, and J.
Elf
, Proc. Natl. Acad. Sci. U.S.A.
107
, 19820
(2010
).15.
R.
Erban
and S.
Chapman
, Phys. Biol.
6
, 046001
(2009
).16.
S.
Hellander
, A.
Hellander
, and L.
Petzold
, Phys. Rev. E
85
, 042901
(2012
).17.
The stepping algorithm implied by the Langevin model of molecular diffusion is derived in the book cited in Ref. 9, and is given explicitly in that book's Eq. (9.22). A generic version of this algorithm was first presented in a more general mathematical context in
D.
Gillespie
, Phys. Rev. E
54
, 2084
(1996
).The Einstein algorithm (16) evidently updates the position of the solute molecule from time t to time t + Δt through a single formula of the form, here, g is a simple function of D and Δt, and n is a “unit normal” random number—i.e., a sample of the normal random variable with mean 0 and variance 1. In contrast, the Langevin stepping algorithm updates the position and the velocity of the solute molecule using two coupled formulas of the form Here, τ = mD/(kBT), fi, and gi are explicit (and in some cases rather complicated) functions of their respective arguments; and nα and nβ are two independent unit normal random numbers. Importantly, the Langevin updating formulas are physically accurate for all Δt > 0, whereas the Einstein updating formula is physically accurate only if Δt ≫ τ. For Δt ≪ τ, the Langevin updating formula for x reduces (approximately) to xt + Δt = xt + vt · Δt, and in that case there is no question about what the value of xt′ was for any t′ ∈ (t, t + Δt); it was xt′ = xt + vt · (t′ − t).
\begin{equation*}x_{t+\Delta t}=x_t +g(D, \Delta t)\cdot n;\end{equation*}
\begin{equation*}\left\{ {\begin{array}{l} {v_{t + \Delta t} = f_1 (\tau,\Delta t) \cdot v_t + f_2 (D,\tau,\Delta t) \cdot n_\alpha {\rm , }} \\[6pt] {x_{t + \Delta t} = x_t + g_1 (\tau,\Delta t) \cdot v_t + g_2 (D,\tau,\Delta t) \cdot n_\alpha + g_3 (D,\tau,\Delta t) \cdot n_\beta.} \\\end{array}} \right.\end{equation*}
18.
K.
Takahashi
, S.
Tănase-Nicola
and P. R.
ten Wolde
, Proc. Natl. Acad. Sci. U.S.A.
107
, 2473
(2010
), see its Supplemental Information section.© 2014 AIP Publishing LLC.
2014
AIP Publishing LLC
You do not currently have access to this content.
