In this blog post I describe the advantages of taking a stochastic view of chemical systems based on the work of D. T. Gillespie and subsequent publications. Gillespie presented his formalism for considering stochastic chemical kinetics, now referred to as the Gillespie Algorithm, in two papers published in 1976 and 1977 (Gillespie, D. T. J. Comp. Phys. 1976, Gillespie D. T. J. Phys. Chem. 1977) – if you want to see the full derivation for the Gillespie Algorithm along with many examples I recommend giving them both a read.
The essential question of chemical kinetics as stated by Gillespie is:
“If a fixed volume V contains a spatially uniform mixture of N chemical species which can inter-react through M specified chemical reaction channels, then given the numbers of molecules of each species present at some initial time, what will these molecular population levels be at any later time?”
Parsing this rather technical language, Gillespie is essentially saying that chemical kinetics is about predicting the time evolution of a chemical system. If we mix some amounts of some different things that react with one another in a flask, and then come back some time later, can we predict the new composition of the flask? Answering these types of questions proves important across basic and applied chemistry. For example, characterizing the kinetics of protein/drug interactions and the lifetimes of drugs in human cells is an important consideration during drug design. Chemical kinetics provides a theoretical framework to understand the timescales of chemical reactions.
There are two different theoretical frameworks used to understand chemical kinetics, typically referred to as (1) deterministic and (2) stochastic chemical kinetics. The term classical chemical kinetics is used interchangeably with deterministic chemical kinetics. In the deterministic formulation (Figure 1, left), the rates of change of the quantities of chemical species {1, 2, …, N} are described by the time derivatives of their concentrations {X1, X2, …, XN}, giving a set of coupled differential equations for {dX1/dt, dX2/dt, dXN/dt}. By using differential equations to describe the time-dependence of molecular concentrations, we make two critical assumptions – that the concentrations {X1, X2, …, XN} are continuous functions of time, and that the time evolution of the system is a deterministic process.
Gillespie argues that these assumptions are actually quite poor for a chemical system; we know, a priori, that the quantities of the chemical species cannot be continuous, because the system is a collection of a large number of discrete molecules, and it is easily observable that chemical processes are inherently stochastic and not deterministic.
Seeking to do away with these assumptions, Gillespie proposed a method for stochastic simulation of the time evolution of a chemical system. Rather than write down a set of differential equations, we consider a collection of single molecules that can react by one or more pathways. There are two key questions – which of the available reaction pathways will be taken? how long will this reaction take? In the Gillespie Algorithm, we answer these questions by randomly selecting results from the joint probability density function P(τ, μ), which gives the probability that reaction μ will take place in the next time interval (t+τ+dτ). Generating functions that provide values of τ and μ for any given reaction scheme are shown in Figure 1 (right panel). In practice, all you need to generate values of τ and μ are one unit-interval random number each (r1 and r2). Note that the various parameters represented by an “a” in the right panel of Figure 1 correspond to the rates of reactions pathways µ. With modern programming it is trivial to generate such random numbers, making the implementation of the Gillespie Algorithm fast and easy.
While the generating functions for τ and μ in Figure 1 appear somewhat intimidating they have very simple interpretations. For τ, we are just randomly sampling a time from an exponential distribution with mean 1/a0. In the case of μ, we can visualize the selection process using a number line extending from 0 to 1. Each of the possible reaction pathways are assigned a portion of this number line corresponding to the rate of their reaction normalized by the sum of the possible reaction rates (denoted a0). If the random number r2 falls in the interval assigned to a particular reaction then that reaction pathway is taken. In a case with two possible reaction pathways out of the current state that have equivalent reaction rates, we will thus also have equivalent probabilities.
Gillespie claims that this view of chemical kinetics is always superior to the classical view because it dispenses with the assumptions of continuity and determinism. One negative aspect of the stochastic formulation of chemical kinetics is that many reactions must be simulated which can take more time than classical chemical kinetics; with classical kinetics, once the differential equations are solved the time-dependence of all concentrations can be easily computed. In practice, however, the sets of differential equations may have no known analytical solution, making numerical methods necessary regardless.
Personally, I have found the Gillespie Algorithm useful when checking results from other models. For example, imagine a situation in which you have extracted a first passage time distribution for the folding time of a protein using molecular dynamics simulations and you would like to understand more about its folding kinetics. Is this protein a two-state folder? Does it appear to be populating an intermediate? Using the Gillespie Algorithm, you can test various reaction schemes to determine which one gives the best agreement with the molecular dynamics data. Stochastic chemical kinetics also has a strong connection to the framework of statistical mechanics from thermal physics, with reaction pathways representing random transitions between accessible microstates of the system.
