Stochastic simulations of chemical reaction kinetics

Attempts to develop a numerical solution of the master equation invariably face insurmountable difficulties. In principle, the probability distribution can be expressed in terms of its moments, and the master equation be equivalently written as a set of ordinary differential equations involving these probability moments. In this scheme, unless severe approximations are made, higher moments are present in the equation for each moment. A closure scheme then is not feasible and one is left with an infinite number of coupled equations. Numerous, approximate closure schemes have been proposed in the literature, hut none has proven adequate. 

One can assume a probability distribution type, e.g., Gaussian, and eliminate high-order moments. In such a case the master equation reduces to a Fokker-Planck equation, as discussed in Chapter 13, for which a numerical solution is often available. Fokker-Planck equations, however, cannot capture the behavior of systems that escape the confines of normal distributions. 

Practically, one is left with methods that sample the probability distribution and its changes in time. In this chapter, we discuss elements of numerical methods to simulate stochastic processes. We first present 

---

Figure 18.1 Regimes of the problem space for multiscale stochastic simulations of chemical reaction kinetics. The r-axis represents the number of molecules of reacting species, x, and the ) -axis measures the frequency of reaction events, A. The threshold variables demarcate the partitions of modeling formalisms. In area I, the number of molecules is so small and the reaction events are so infrequent that a discrete-stochastic simulation algorithm, like the SSA, is needed. In contrast, in area V, which extends to infinity, the thermodynamic limit assumption becomes vahd and a continuous-deterministic modehng formalism becomes valid. Other areas admit different modehng formalisms, such as ones based on chemical Langevin equations, or probabilistic steady-state assumptions.

For near and supercritical conditions, combustion gas-phase data are often used as the point of reference to assess solvent effects. The gas-phase values of kig, available for temperatures 800-2500 K, show the activation energy 90 kJ mol In condensed phase, stabilization of H2O molecules via H-bonding may increase the activation barrier, but on the other hand the reaction can be promoted by the solvent cage effect. Diffusion-kinetic modelling and stochastic simulation of chemical reactions in radiation tracks have shown that the occurrence of reaction (15.19) is consistent with the anomalous increase in H2 yield observed in water radiolysis at temperatures above 523 K, if kig is of the order of 1-2x10 s (4-8x10 s ) at 573 K. Considering the two... 

Reaction kinetic models can be simulated not only by solving the kinetic system of differential equations but also via simulating the equivalent stochastic models. Computer codes are available that solve the stochastic kinetic equations. One of these is the Chemical Kinetics Simulator (CKS) program that was developed at IBM s Almaden Research Centre. It provides a rapid, interactive method for the accurate simulation of chemical reactions. CKS is a good tool for teaching the principles of stochastic reaction kinetics to students and trainees. 

This paper has focused on two recent computer methods for discrete simulation of chemical kinetics. Beginning with the realization that truly microscopic computer experiments are not at all feasible, I have tried to motivate the development of a hierarchy of simulations in studies of a class of chemical problems which best illustrate the absolute necessity for simulation at levels above molecular dynamics. It is anticipated (optimistically ) that the parallel development of discrete event simulations at different levels of description may ultimately provide a practical interface between microscopic physics and macroscopic chemistry in complex physicochemical systems. With the addition to microscopic molecular dynamics of successively higher-level simulations intermediate between molecular dynamics at one extreme and differential equations at the other, it should be possible to examine explicitly the validity of assumptions invoked at each stage in passing from the molecular level to the stochastic description and finally to the macroscopic formulation of chemical reaction kinetics. 

Although the collision and transition state theories represent two important methods of attacking the theoretical calculation of reaction rates, they are not the only approaches available. Alternative methods include theories based on nonequilibrium statistical mechanics, stochastic theories, and Monte Carlo simulations of chemical dynamics. Consult the texts by Johnson (62), Laidler (60), and Benson (59) and the review by Wayne (63) for a further introduction to the theoretical aspects of reaction kinetics. 

Ray Kapral came to Toronto from the United States in 1969. His research interests center on theories of rate processes both in systems close to equilibrium, where the goal is the development of a microscopic theory of condensed phase reaction rates,89 and in systems far from chemical equilibrium, where descriptions of the complex spatial and temporal reactive dynamics that these systems exhibit have been developed.90 He and his collaborators have carried out research on the dynamics of phase transitions and critical phenomena, the dynamics of colloidal suspensions, the kinetic theory of chemical reactions in liquids, nonequilibrium statistical mechanics of liquids and mode coupling theory, mechanisms for the onset of chaos in nonlinear dynamical systems, the stochastic theory of chemical rate processes, studies of pattern formation in chemically reacting systems, and the development of molecular dynamics simulation methods for activated chemical rate processes. His recent research activities center on the theory of quantum and classical rate processes in the condensed phase91 and in clusters, and studies of chemical waves and patterns in reacting systems at both the macroscopic and mesoscopic levels. 

Sipos, T., Toth, J. Erdi, P. (1974a). Stochastic simulation of complex chemical reactions by digital computer, I. The model. React. Kinet. Catal. Lett., 1, 113-17. Sipos, T., T6th, J. Erdi, P. (1974b). Stochastic simulation of complex chemical reactions by digital computer, II. Applications. React. Kinet. Catal. Lett., 1, 209-13. 

In a continuous-flow chemical reactor, the concern is not only with probabilistic transitions among chemical species but also with probabilistic liansitions of each chemical species between the interior and exterior of the reactor. Pippel and Philipp [8] used Markov chains for simulating the dynamics of a chemical system. In their approach, the kinetics of a chemical reaction are treated deterministically and the flow through the system are treated stochastically by means of a Markov chain. Shinnar et al. [9] superimposed the kinetics of the first order chemical reactions on a stochastically modeled mixing process to characterize the performance of a continuous-flow reactor and compared it with that of the corresponding batch reactor. Most stochastic approaches to analysis and modeling of chemical reactions in a flow system have combined deterministic chemical kinetics and stochastic flows. 

The irradiation of water is immediately followed by a period of fast chemistry, whose short-time kinetics reflects the competition between the relaxation of the nonhomogeneous spatial distributions of the radiation-induced reactants and their reactions. A variety of gamma and energetic electron experiments are available in the literature. Stochastic simulation methods have been used to model the observed short-time radiation chemical kinetics of water and the radiation chemistry of aqueous solutions of scavengers for the hydrated electron and the hydroxyl radical to provide fundamental information for use in the elucidation of more complex, complicated chemical, and biological systems found in real-world scenarios. 

It is often stated that MC methods lack real time and results are usually reported in MC events or steps. While this is immaterial as far as equilibrium is concerned, following real dynamics is essential for comparison to solutions of partial differential equations and/or experimental data. It turns out that MC simulations follow the stochastic dynamics of a master equation, and with appropriate parameterization of the transition probabilities per unit time, they provide continuous time information as well. For example, Gillespie has laid down the time foundations of MC for chemical reactions in a spatially homogeneous system.f His approach is easily extendable to arbitrarily complex computational systems when individual events have a prescribed transition probability per unit time, and is often referred to as the kinetic Monte Carlo or dynamic Monte Carlo (DMC) method. The microscopic processes along with their corresponding transition probabilities per unit time can be obtained via either experiments such as field emission or fast scanning tunneling microscopy or shorter time scale DFT/MD simulations discussed earlier. The creation of a database/lookup table of transition... 

Finally, an entirely different approach to simulating gelation is the Dynamic Monte Carlo (DMC) method, in which chemical reactions are modeled by stochastic integration of phenomenological kinetic rate laws [23]. This has been used successfully to understand the onset of gel formation, first-shell substitution effects, and the influence of cyclization in silicon alkoxide systems [24—26]. However, this approach has not so far been extended to include the instantaneous positions and diffusion of each oUgomer, which would be necessary in order for the calculation to generate an actual model of an aerogel that could be used in subsequent simulations. 

E. L. Haseltine and J. B. Rawlings, Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics, J. Chem. Phys., 117, 6959-6969, (2002). 

In this section, we will present results of microldnetics simulations based on elementary reaction energy schemes deduced from quantum chemical studies. We use an adapted scheme to enable analysis of the results in terms of the values of elementary rate constants selected. For the same reason, we ignore surface concentration dependence of adsorption energies, whereas this can be readily implemented in the simulations. We are interested in general trends and especially the temperature dependence of overall reaction rates. The simulations will also provide us with information on surface concentrations. In the simulations to be presented here, we exclude product readsorption effects. Microldnetics simulations are attractive since they do not require an assumption of rate-controlling steps or equilibration. Solutions for overall rates are found by solving the complete set of PDFs with proper initial conditions. While in kinetic Monte Carlo simulations these expressions are solved using stochastic techniques, which enable formation... 

Gillespie s algorithm numerically reproduces the solution of the chemical master equation, simulating the individual occurrences of reactions. This type of description is called a jump Markov process, a type of stochastic process. A jump Markov process describes a system that has a probability of discontinuously transitioning from one state to another. This type of algorithm is also known as kinetic Monte Carlo. An ensemble of simulation trajectories in state space is required to accurately capture the probabilistic nature of the transient behavior of the system. 

---