Lecture notes. Written from the material of MIT 8.591J 2004, licensed CC BY-NC-SA 4.0. These are notes, not a transcript: the material has been reorganised and rewritten. This adaptation carries the same licence, and the original is linked above.
10. Stochastic Chemical Kinetics¶
What this covers¶
How should we describe a chemical reaction network when some of the reacting species are present in only tens or hundreds of copies — too few for the deterministic rate equations of ordinary chemical kinetics to be trusted? This chapter builds the probabilistic description from scratch: the master equation that replaces a rate equation for discrete molecule counts, the two limits in which it reduces to something familiar (the deterministic law for the mean, and a diffusion-like Fokker–Planck equation for large-but-not-infinite numbers), and the simulation algorithm and switching-time calculation that the master equation makes possible. It assumes the reader is comfortable with ordinary first-order chemical kinetics (\(dn/dt = k-\gamma n\)), with means and variances, and with solving simple ODEs.
Why molecule number cannot be treated as continuous¶
Take the simplest possible reaction network: a molecule \(X\) is created at a constant rate \(k\) and destroyed by a first-order reaction with rate \(\gamma\). If the number of molecules \(n\) is enormous — Avogadro's number, about \(10^{24}\), is the usual scale in a beaker — it costs nothing to pretend \(n\) is a continuous variable and write the ordinary rate equation
where \(f_n = k\) is the creation rate and \(g_n = \gamma n\) the destruction rate when there are \(n\) molecules present.
Inside a living cell this approximation is far less innocent: a typical protein is present in thousands of copies, a ribosome or an RNA polymerase in hundreds, an mRNA species in tens, and most genes in one or two copies. So take discreteness seriously and ask what the exact trajectory \(n(t)\) looks like.
One tempting guess is that \(n(t)\) still moves smoothly through the integers, with a creation event exactly every \(\Delta t = \Delta n/(k-\gamma n)\) for \(\Delta n = 1\) — a "clockwork" version of the same curve. This is not physically possible: it would require the system to keep track of the time elapsed since the last event, the way a clock counts ticks, but there is no clock here, only molecules colliding with each other. The system's future can depend only on its present state, never on how it got there. (No memory of the past beyond the present state is the defining property of a Markov process.)
What a memoryless system can do is have each possible reaction occur with some probability per unit time, proportional to its rate. Run the same experiment twice from the same initial condition and the two trajectories \(n(t)\) will differ — the system is genuinely stochastic. Run it many times and record \(n\) at some fixed time, and the outcomes form a distribution rather than a single value. The rest of the chapter is about that distribution.
The master equation¶
From a rate to a probability¶
If a given reaction fires at rate \(r\), then over a long stretch of time \(T\) it fires on average \(rT\) times. Split \(T\) into \(N\) equal slivers of length \(dt = T/N\); the chance the reaction happens in any one particular sliver is \(rT/N = r\,dt\). So in the limit of small \(dt\): the probability that a reaction with rate \(r\) occurs in a time interval \(dt\) is \(r\,dt\), and the chance of two occurrences within the same sliver is negligible.
Setting up an ensemble¶
Imagine many identical copies of the system, all started from the same initial condition, and let \(p_n(t)\) be the fraction of them that have exactly \(n\) molecules of \(X\) at time \(t\). Four things can change \(p_n\) in a small interval \(dt\):
- a system with \(n-1\) molecules gains one, at rate \(f_{n-1}\) — moving into the \(n\)-population;
- a system with \(n+1\) molecules loses one, at rate \(g_{n+1}\) — moving into the \(n\)-population;
- a system with \(n\) molecules gains one, at rate \(f_n\) — moving out of the \(n\)-population;
- a system with \(n\) molecules loses one, at rate \(g_n\) — moving out of the \(n\)-population.
Collecting all four fluxes gives the master equation:
This is really an infinite family of coupled linear ODEs, one for every \(n = 0,1,2,\dots\). Because it is linear, it does not matter whether \(p_n\) is a raw count of systems or a normalized probability with \(\sum_n p_n = 1\) — the second reading is the useful one: \(p_n(t)\) is the probability that a given system has exactly \(n\) molecules at time \(t\). Concretely, the ensemble could be a population of cells and \(p_n\) the fraction of cells carrying \(n\) copies of some protein.
The deterministic limit reappears in the mean¶
The master equation looks like a step backward — instead of one ODE for \(n(t)\) there are now infinitely many, for the whole distribution \(p_n(t)\). But its moments can be extracted without ever solving it. Take the mean number of molecules,
differentiate, substitute the master equation, and reindex the sums (using \(\sum_n h(n) = \sum_n h(n\pm1)\), valid whenever the sum runs over all integers). With the specific rates \(f_n = k\), \(g_n = \gamma n\) of the birth–death example, every term but a boundary piece cancels and what remains is
The mean molecule number obeys exactly the deterministic rate equation, even though individual systems fluctuate around it. This is not a special feature of this example: the mean always obeys the deterministic equation whenever the rates \(f_n, g_n\) are linear functions of \(n\).
The steady-state distribution: a Poisson¶
Set \(dp_n/dt = 0\) in the master equation for the birth–death system:
Write \(\bar n = k/\gamma\) (the deterministic steady state) and rearrange:
The right-hand side evaluated at \(n\) is exactly the left-hand side evaluated at \(n-1\), so this quantity is the same number for every \(n\) — call it \(C\). Normalizability forces \(C = 0\): a nonzero constant "current" between neighboring states cannot be sustained all the way out to \(n\to\infty\) by a distribution that sums to something finite. With \(C=0\) the recursion collapses:
and normalizing (using \(\sum_n \bar n^n/n! = e^{\bar n}\)) fixes \(p_0 = e^{-\bar n}\). The result is the Poisson distribution:
The argument that forced \(C=0\) will reappear, in a different guise, in the Fokker–Planck picture below: a genuine steady state cannot sustain a constant nonzero probability current.
How large is "large"?¶
For a Poisson distribution the mean and the variance coincide:
so the relative size of the fluctuations is
This is the precise sense in which molecule counts need to be "large" before the deterministic equation can be trusted: the relative deviation from the deterministic trajectory shrinks only as the inverse square root of the number of molecules. An ensemble of systems averaging 20 molecules shows about a 22% spread around that value; one averaging 500 molecules shows only about 4%. The two regimes look qualitatively different even though both obey the same equation for the mean.
The Fokker–Planck equation: a continuum approximation¶
Solving the master equation exactly means tracking every \(p_n\) for every \(n\) — an infinite system. There is a useful shortcut for intermediate molecule numbers: large enough that the step from \(n\) to \(n\pm1\) is a small perturbation, but not so large that fluctuations can be ignored altogether. Treat \(n\) as continuous and write \(h(n)\) for what was \(h_n\). Any smooth function Taylor-expands as
so expanding the two "neighbor" terms of the master equation around \(n\) (here \(\Delta n = \pm1\)),
Substituting into the master equation and keeping terms through second order gives the Fokker–Planck equation:
where \(J\) is a probability flux. The equation is a diffusion equation in disguise: think of each system in the ensemble as a particle whose position is its molecule number \(n\), and \(J\) as the current of particles across a given value of \(n\).
The zero-flux argument, and a potential¶
In steady state \(J\) must be constant in \(n\) — otherwise probability would pile up without bound somewhere. But the flux at \(n=0\) has to vanish, since no system can flow into negative molecule number, so \(J \equiv 0\) everywhere, not just at that one boundary. Zero flux means
Substituting \(q = (f+g)p\) turns this into a first-order linear ODE for \(q\):
so that
The steady-state probability is, up to the prefactor \(1/(f+g)\), a Boltzmann-like weight \(e^{-\varphi(n)}\) for an effective potential \(\varphi\) built entirely out of the reaction rates: wherever \(f>g\) (net production) \(\varphi\) decreases, and wherever \(g>f\) (net removal) it increases, so \(\varphi\) has minima exactly at the deterministic steady states.
Example: a bistable genetic switch¶
Take the autocatalytic gene-regulatory system from problem set 1:
There \(x\) was a concentration and the decay term came from dilution by cell growth. Here take the cell volume fixed and let the decay instead be an active degradation process, so that \(x\) can be reinterpreted directly as a molecule number. With \(\gamma=1\), \(K_1K_2=10^{-4}\), \(v_0=12.5\), \(v_1=200\), \((f-g)(x)\) crosses zero three times — two stable roots flanking one unstable root — so the system is bistable.
Feeding this \(f\) and \(g\) into \(\varphi(x) = -\int 2\,(f-g)/(f+g)\,dx'\) produces a double-well potential, with the two stable states sitting in the wells and the unstable state at the barrier between them. The corresponding steady-state distribution \(p(x) = \frac{A}{f+g}e^{-\varphi(x)}\) is bimodal: peaked in each well, and vanishingly small at the barrier.
This is the picture the rest of the chapter builds toward: a population of otherwise identical cells running this network splits into two subpopulations sitting in the two wells, with individual cells occasionally hopping between them.
How long between reactions? The waiting-time distribution¶
The master equation and its continuum relative describe the distribution \(p_n(t)\), but to actually simulate one system's trajectory — rather than the whole ensemble — it helps to know the distribution of the time until the next reaction, for a reaction firing at constant rate \(r\).
Let \(P(\tau)\,d\tau\) be the probability that the next occurrence falls between \(\tau\) and \(\tau+d\tau\). This factors as "does not happen before \(\tau\)" times "happens in the next \(d\tau\)":
\(Q\) satisfies its own recursion, since not occurring before \(\tau\) means not occurring before \(\tau - d\tau\) and not occurring in the last sliver:
(using \(Q(0)=1\)). So
the waiting time between successive firings of a constant-rate reaction is exponentially distributed, with mean \(\langle\tau\rangle = 1/r\) and variance \(\langle\delta\tau^2\rangle = 1/r^2\). This is the same memorylessness that made \(n(t)\) a Markov process in the first place: the exponential is the only distribution for which how long you have already waited carries no information about how much longer you will wait.
Simulating a reaction network: the Gillespie recipe¶
The exponential waiting-time law is exactly what is needed to simulate individual stochastic trajectories numerically, rather than solve for the whole distribution \(p_n(t)\).
Sampling trick. If \(u\) is drawn uniformly from \((0,1)\), then
is distributed exactly as the waiting time \(\tau\) above — exponential with rate \(r\). (This is the standard inverse-transform trick: it turns a uniform random-number generator into a generator for any distribution whose cumulative distribution function can be inverted.)
The algorithm. Suppose the network has \(m\) possible reaction channels — \(m=2\) for the birth–death system, creation and destruction — each currently firing at rate \(r_i\):
- Draw \(m\) independent putative waiting times \(\theta_i\) from the recipe above, one per channel.
- Let the channel with the smallest \(\theta_i\) be the one that actually fires next; advance the clock by that \(\theta_i\).
- Update the state accordingly — e.g. \(n \to n+1\) for a creation event, \(n \to n-1\) for a destruction event.
- Recompute all the rates \(r_i\) at the new state, and repeat from step 1.
Running this once generates one noisy trajectory \(n(t)\) of the kind sketched at the start of the chapter. Running it many times — 2,500 repeats, in the case worked through here — and recording only the final state at some fixed time reproduces the distribution of outcomes across the ensemble. This procedure is a slight simplification of the full Gillespie algorithm.
Spontaneous switching in a bistable system¶
Apply the simulation recipe to the bistable network above. A single cell does not sit permanently in one well: its trajectory \(x(t)\) occasionally jumps from one stable state to the other, driven by a run of fluctuations large enough to carry it over the barrier. Across a population, this produces exactly the bimodal histogram of the steady-state distribution derived from the Fokker–Planck potential — most cells sit in one well or the other, in proportions set by \(p(x)\), with rare individuals caught mid-transition.
The natural quantity to ask for is the lifetime of a state: the average time a cell spends in one well before escaping to the other. Simulation shows that this lifetime grows rapidly with the average number of \(X\) molecules in that state, and diverges as that number gets large — a second, independent way of seeing the crossover to deterministic behavior: with enough molecules, a state that was only metastable becomes, for all practical purposes, permanently stable.
This is not a purely academic question. Reading the switch (with time measured in units of the degradation time of \(X\), typically of order a cell lifetime) as a toy model of the phage-\(\lambda\) lysis/lysogeny decision, the induced state corresponds to lysogeny and an escape event to spontaneous lysis. Measured lysogens undergo spontaneous lysis at a rate of about \(10^{-8}\) per cell per generation — an extremely stable switch. Matching that stability in the simple model above requires roughly 1,000 molecules, consistent with measured repressor concentrations in the real system. Stability can also be bought by means other than raw molecule number: higher cooperativity in the regulatory interactions, or stabilizing feedback loops, both deepen the wells and steepen the barrier without needing more molecules.
Sources¶
All of this chapter comes from the lecture-outline handout for Lecture 12 of 7.81/8.591/9.531
Systems Biology (A. van Oudenaarden and M. Thattai, MIT, Fall 2004; MIT OpenCourseWare, CC
BY-NC-SA), sections 1–10, split across six converted files in lecture-outlines/12-outline/:
01-1-introduction.md— §1, the discreteness/Markov motivation (source Fig. 1).02-2-probabilistic-formulation-of-reaction-kinetics-the-master.md— §2 (master equation, source Fig. 2), §3 (deterministic limit for the mean), §4 (Poisson steady state), §5 (large-number limit, source Fig. 3).03-6-the-fokker-planck-equation.md— §6 (Fokker–Planck equation) and §7 (bistable-switch example, source Fig. 4).04-8-waiting-times-between-reaction-events.md— §8 (exponential waiting-time distribution).05-9-stochastic-simulation-of-chemical-reactions.md— §9 (Gillespie-style simulation recipe).06-10-spontaneous-switching-rates-in-a-bistable-system.md— §10 (escape times, source Fig. 5, phage-λ discussion).
No transcript, written notes or problem set were supplied for this lecture; the two diagrams in
this chapter (the master-equation flux diagram and the double-well potential) are original
renderings of what the source describes as its Figs. 2 and 4, not reproductions of the original
figures, which were not available. The bistable network's rate function \(f(x)\) is stated in the
source to come from "problem set 1" of the course, which was not among the supplied material and is
not reproduced here. The source files themselves are flagged as model-reconstructed from a PDF with
no text layer (sources/ocw-8591j-2004/lecture-outlines/12-outline.pdf), with every equation
marked unverified in that conversion — the equations above should be checked against the original
PDF before being cited as a primary source.
← 9. Genetic Toggle Switch and Oscillator · Contents · 11. The Turing Instability Condition →