Daniel T. Gillespie

Simulating the evolution of a chemically reacting system using the bimolecular propensity function, as is done by the stochastic simulation algorithm and its reaction-diffusion extension, entails making statistically inspired guesses as to where the reactant molecules are at any given time. Those guesses will be physically justified if the system is dilute and well-mixed in the reactant molecules. Otherwise, an accurate simulation will require the extra effort and expense of keeping track of the positions of the reactant molecules as the system evolves. One molecule-tracking algorithm that pays careful attention to the physics of molecular diffusion is the enhanced Green's function reaction dynamics (eGFRD) of Takahashi, Tănase-Nicola, and ten Wolde [Proc. Natl. Acad. Sci. U.S.A. 107, 2473 (2010)]. We introduce here a molecule-tracking algorithm that has the same theoretical underpinnings and strategic aims as eGFRD, but a different implementation procedure. Called the small-voxel tracking algorithm (SVTA), it combines the well known voxel-hopping method for simulating molecular diffusion with a novel procedure for rectifying the unphysical predictions of the diffusion equation on the small spatiotemporal scale of molecular collisions. Indications are that the SVTA might be more computationally efficient than eGFRD for the problematic class of non-dilute systems. A widely applicable, user-friendly software implementation of the SVTA has yet to be developed, but we exhibit some simple examples which show that the algorithm is computationally feasible and gives plausible results.

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.

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 discrete-stochastic model of diffusion introduced in the preceding chapter implies that the system state evolves in time as what is known as a jump Markov process. The time evolution of such a process is generally governed by what is called a master equation. In this case, the assumption that the solute molecules move about independently of each other gives rise to two master equations, one for a single solute molecule and one for the entire collection of solute molecules. In this chapter the chapter derives both of those master equations. The chapter obtains the exact time-dependent solution of the single-molecule master equation, and the exact time-independent (equilibrium) solution of the many-molecule master equation. The chapter also derives companion stochastic simulation algorithms for the two master equations, and the chapter applies those algorithms to the microfluidic diffusion experiment of Chapter 5. The chapter uses simulation to test the validity of the discrete-stochastic version of Fick's Law.

This chapter introduces a model of diffusion in which the position of the solute molecule is described by a discrete stochastic variable, in contrast to the continuous stochastic variable of the Einstein model. The aim of this new model is to provide a description of diffusion which is for most practical purposes equivalent to the Einstein model, but which affords new ways of describing it mathematically. The system volume is imagined to be subdivided into equal-size cells, and the current state of the system is given by the current numbers of solute molecules in the cells. The diffusion process is modeled by supposing that the solute molecules independently jump to adjacent cells according to a specific probability rule. The chapter establishes consistency between that cell-jumping rule and both the Fickian and the Einstein models of diffusion, and the chapter deduces a generalization of the cell-jumping rule for unequal cell sizes. The chapter shows that there must be some lower bound on the cell size, and it derives a generalization of the standard cell-jumping rule that relaxes that lower bound. The chapter examines in detail the discrete-stochastic model's version of Fick's Law. Finally, the chapter illustrates the model in the context of a laboratory experiment in which the solute molecules are tiny polystyrene beads diffusing in water inside a specially fabricated microfluidic chip.

Chapter 3 showed how Einstein's analysis of Brownian motion shifted the focus of the classical diffusion equation from the average behavior of many solute molecules to the probabilistic behavior of an individual solute molecule. This chapter deduces some further implications of the Einstein theory of diffusion, using not only analytical reasoning but also numerical simulation. This chapter shows how simulation can be used to construct plots of a single diffusing molecule's position probability density function for some simple boundary conditions. The chapter also shows how simulation can be used to construct “snapshots” of an unrestricted solute molecule's position at successive instants of time. But the chapter will discovers that such snapshots expose a serious physical limitation of the Einstein theory. The chapter rationalizes a quick fix that allows us to complete the derivation which has begun in Chapter 3 of a formula for the stochastic rate of a diffusion-controlled bimolecular chemical reaction. But from a broader view, it will become apparent that, while the Einstein theory of diffusion has a wide range of practical utility, a physically more accurate theory is needed.