The small-voxel tracking algorithm for simulating chemical reactions among diffusing molecules.

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, Tanase-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.

Validity conditions for stochastic chemical kinetics in diffusion-limited systems.

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.

Perspective: Stochastic algorithms for chemical kinetics.

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.

A diffusional bimolecular propensity function.

We derive an explicit formula for the propensity function (stochastic reaction rate) of a generic bimolecular chemical reaction in which the reactant molecules move about by diffusion, as solute molecules in a bath of much smaller and more numerous solvent molecules. Our derivation assumes that the solution is macroscopically well stirred and dilute in the solute molecules. It effectively extends the physical rationale for the chemical master equation and the stochastic simulation algorithm from well-stirred dilute gases to well-stirred dilute solutions, with the former becoming a limiting case of the latter. This extension is important for cellular systems, where the solvent molecules are typically water and the solute (reactant) molecules are much larger organic structures, whose relatively low populations often require a discrete-stochastic formalism. In the course of our derivation, we illuminate some limitations on the ability of the classical diffusion equation to accurately describe how a diffusing molecule moves on spatial and temporal scales that are relevant to collision-induced chemical reactions.

Deterministic limit of stochastic chemical kinetics.

An analysis is presented of the approximating assumptions that underlie a recently proposed derivation of the traditional deterministic reaction rate equation from a discrete-stochastic formulation of chemical kinetics. It is shown that if the system is close enough to the thermodynamic limit, in which the molecular populations and the containing volume all approach infinity in such a way that the molecular concentrations remain finite, then the required approximating assumptions will be justified for practically all spatially homogeneous systems that one is likely to encounter.

Algorithms and software for stochastic simulation of biochemical reacting systems.

Traditional deterministic approaches for simulation of chemically reacting systems fail to capture the randomness inherent in such systems at scales common in intracellular biochemical processes. In this manuscript, we briefly review the state of the art in discrete stochastic and multiscale algorithms for simulation of biochemical systems and we present the StochKit software toolkit.

Adaptive explicit-implicit tau-leaping method with automatic tau selection.

The existing tau-selection strategy, which was designed for explicit tau leaping, is here modified to apply to implicit tau leaping, allowing for longer steps when the system is stiff. Further, an adaptive strategy that identifies stiffness and automatically chooses between the explicit and the (new) implicit tau-selection methods to achieve better efficiency is proposed. Numerical testing demonstrates the advantages of the adaptive method for stiff systems.

Stochastic simulation of chemical kinetics.

Stochastic chemical kinetics describes the time evolution of a well-stirred chemically reacting system in a way that takes into account the fact that molecules come in whole numbers and exhibit some degree of randomness in their dynamical behavior. Researchers are increasingly using this approach to chemical kinetics in the analysis of cellular systems in biology, where the small molecular populations of only a few reactant species can lead to deviations from the predictions of the deterministic differential equations of classical chemical kinetics. After reviewing the supporting theory of stochastic chemical kinetics, I discuss some recent advances in methods for using that theory to make numerical simulations. These include improvements to the exact stochastic simulation algorithm (SSA) and the approximate explicit tau-leaping procedure, as well as the development of two approximate strategies for simulating systems that are dynamically stiff: implicit tau-leaping and the slow-scale SSA.

Efficient step size selection for the tau-leaping simulation method.

The tau-leaping method of simulating the stochastic time evolution of a well-stirred chemically reacting system uses a Poisson approximation to take time steps that leap over many reaction events. Theory implies that tau leaping should be accurate so long as no propensity function changes its value "significantly" during any time step tau. Presented here is an improved procedure for estimating the largest value for tau that is consistent with this condition. This new tau-selection procedure is more accurate, easier to code, and faster to execute than the currently used procedure. The speedup in execution will be especially pronounced in systems that have many reaction channels.

Accelerated stochastic simulation of the stiff enzyme-substrate reaction.

The enzyme-catalyzed conversion of a substrate into a product is a common reaction motif in cellular chemical systems. In the three reactions that comprise this process, the intermediate enzyme-substrate complex is usually much more likely to decay into its original constituents than to produce a product molecule. This condition makes the reaction set mathematically "stiff." We show here how the simulation of this stiff reaction set can be dramatically speeded up relative to the standard stochastic simulation algorithm (SSA) by using a recently introduced procedure called the slow-scale SSA. The speedup occurs because the slow-scale SSA explicitly simulates only the relatively rare conversion reactions, skipping over occurrences of the other two less interesting but much more frequent reactions. We describe, explain, and illustrate this simulation procedure for the isolated enzyme-substrate reaction set, and then we show how the procedure extends to the more typical case in which the enzyme-substrate reactions occur together with other reactions and species. Finally, we explain the connection between this slow-scale SSA approach and the Michaelis-Menten [Biochem. Z. 49, 333 (1913)] formula, which has long been used in deterministic chemical kinetics to describe the enzyme-substrate reaction.

Avoiding negative populations in explicit Poisson tau-leaping.

The explicit tau-leaping procedure attempts to speed up the stochastic simulation of a chemically reacting system by approximating the number of firings of each reaction channel during a chosen time increment tau as a Poisson random variable. Since the Poisson random variable can have arbitrarily large sample values, there is always the possibility that this procedure will cause one or more reaction channels to fire so many times during tau that the population of some reactant species will be driven negative. Two recent papers have shown how that unacceptable occurrence can be avoided by replacing the Poisson random variables with binomial random variables, whose values are naturally bounded. This paper describes a modified Poisson tau-leaping procedure that also avoids negative populations, but is easier to implement than the binomial procedure. The new Poisson procedure also introduces a second control parameter, whose value essentially dials the procedure from the original Poisson tau-leaping at one extreme to the exact stochastic simulation algorithm at the other; therefore, the modified Poisson procedure will generally be more accurate than the original Poisson procedure.

The slow-scale stochastic simulation algorithm.

Reactions in real chemical systems often take place on vastly different time scales, with "fast" reaction channels firing very much more frequently than "slow" ones. These firings will be interdependent if, as is usually the case, the fast and slow reactions involve some of the same species. An exact stochastic simulation of such a system will necessarily spend most of its time simulating the more numerous fast reaction events. This is a frustratingly inefficient allocation of computational effort when dynamical stiffness is present, since in that case a fast reaction event will be of much less importance to the system's evolution than will a slow reaction event. For such situations, this paper develops a systematic approximate theory that allows one to stochastically advance the system in time by simulating the firings of only the slow reaction events. Developing an effective strategy to implement this theory poses some challenges, but as is illustrated here for two simple systems, when those challenges can be overcome, very substantial increases in simulation speed can be realized.

The numerical stability of leaping methods for stochastic simulation of chemically reacting systems.

Tau-leaping methods have recently been proposed for the acceleration of discrete stochastic simulation of chemically reacting systems. This paper considers the numerical stability of these methods. The concept of stochastic absolute stability is defined, discussed, and applied to the following leaping methods: the explicit tau, implicit tau, and trapezoidal tau.

Modeling quantum measurement probability as a classical stochastic process.

The time-dependent measurement probabilities for the simple two-state quantum oscillator seem to invite description as a classical two-state stochastic process. It has been shown that such a description cannot be achieved using a Markov process. Constructing a more general non-Markov process is a challenging task, requiring as it does the proper generalizations of the Markovian Chapman-Kolmogorov and master equations. Here we describe those non-Markovian generalizations in some detail, and we then apply them to the two-state quantum oscillator. We devise two non-Markovian processes that correctly model the measurement statistics of the oscillator, we clarify a third modeling process that was proposed earlier by others, and we exhibit numerical simulations of all three processes. Our results illuminate some interesting though widely unappreciated points in the theory of non-Markovian stochastic processes. But since quantum theory does not tell us which one of these quite different modeling processes "really" describes the behavior of the oscillator, and also since none of these processes says anything about the dynamics of other (noncommuting) oscillator observables, we can see no justification for regarding any of these processes as being fundamentally descriptive of quantum dynamics. (c) 2001 American Institute of Physics.