Perspective: Stochastic algorithms for chemical kinetics

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.

Figure 1

Stochastic chemical kinetics is premised on the definition 2 of the propensity function in the top box, a definition which must look to molecular physics for its justification. The two solid-outlined boxes in yellow denote mathematically exact consequences of that definition: the chemical master equation 1 and the stochastic simulation algorithm 4. Dashed-outlined boxes denote approximate consequences: tau-leaping 7, the chemical Langevin equation 9, the chemical Fokker-Planck equation (not discussed here but see Ref. 35), and the reaction rate equation 10. The bracketed condition by each dashed inference arrow is the condition enabling that approximation: reading from top to bottom, those conditions are the first leap condition, the second leap condition, and the thermodynamic limit. The rationale for viewing the linear noise approximation (LNA) as an intermediate result between the CLE and the RRE is detailed in Ref. . It has been shown, that for realistic propensity functions, getting “close enough” to the thermodynamic limit will ensure simultaneous satisfaction of the first and second leap conditions, at least for finite spans of time; therefore, the top-to-bottom progression indicated in the figure will inevitably occur as the molecular populations and the system volume become larger. But a given chemical system might be such that the largest value of τ that satisfies the first leap condition will not be large enough to satisfy the second leap condition; in that case, there will be no accurate description of the system below the discrete-stochastic level in the figure.

Figure 2

Parts of a Cartesian mesh (a) and an unstructured triangular mesh (b). Molecules are assumed to be well-mixed in the local volumes that make up the dual elements of the mesh (depicted in pink color). For the Cartesian grid (a), the dual is simply the staggered grid. The dual of the triangular mesh in (b) is obtained by connecting the midpoints of the edges and the centroids of the triangles. (c) shows how a model of a eukaryotic cell with a nucleus (green) can be discretized with a mesh made up of triangles and tetrahedra. The figure is adapted from Ref. , where a model of nuclear import was simulated on this domain using the URDME software.

Figure 3

Schematic representation of the RDME's behavior as a function of the voxel size h. For h < h*, no local correction to the conventional mesoscopic reaction rates exists that will make the RDME consistent with the Smoluchowski model for the simple problem of diffusion to a target. Figure adapted from Ref. .