Stochastic Processes in Biology

⚡ Key Takeaways

Stochastic processes in biology model biological systems as collections of random variables evolving in time, capturing the inherent noise of molecular, cellular, ecological, and evolutionary dynamics. Key result: when molecule numbers are small, deterministic rate equations fail, and one must use the Chemical Master Equation (CME) or its diffusion approximations (Fokker–Planck, Langevin) to predict fluctuations, switching, and extinction. Applies to: gene expression, ion channels, population genetics, neural spiking, epidemics, enzyme kinetics. Who/when: rooted in Markov (1906), Wright–Fisher (1930s), Delbrück (1940), Gillespie (1976), and modern single-cell biology (2000s–).

💡 The Core Insight

When the number of interacting entities (molecules, individuals) is small, the intrinsic randomness of individual reaction or birth/death events produces fluctuations of relative magnitude $\sim 1/\sqrt{N}$ . These fluctuations are not averaged away and qualitatively change behavior: they create finite extinction times, multimodal phenotype distributions, stochastic switching between attractors, and noise-induced transitions absent from the deterministic description. The correct level of description is a Markov jump process on integer states, whose probability distribution obeys the Chemical Master Equation.

DefinitionMarkov Process and Transition Semigroup

A stochastic process $\{X(t)\}_{t\ge 0}$ taking values in a state space $S$ is a Markov process if for all $t_1<\dots<t_n<t$ and states,

1

\quad P\big(X(t)\in A \mid X(t_n),\dots,X(t_1)\big)=P\big(X(t)\in A \mid X(t_n)\big).

Derivation hint: This is the definition of the Markov property — the future depends on the present only, not the path. Conditioning collapses to the most recent observation. For a continuous-time Markov chain on a countable state $S$ , the dynamics are encoded by the transition rate matrix (generator) $Q=(q_{ij})$ with $q_{ij}\ge 0$ ( $i\ne j$ ) and $q_{ii}=-\sum_{j\ne i}q_{ij}$ . The transition probabilities $P_{ij}(t)=P(X(t)=j\mid X(0)=i)$ satisfy

2

\quad \frac{d}{dt}P(t)=P(t)\,Q=Q\,P(t),\qquad P(t)=e^{tQ}.

Derivation hint: Combine the Chapman–Kolmogorov equation $P(t+s)=P(t)P(s)$ with $P(h)\approx I+hQ$ for small $h$ , then take $h\to 0$ to get the forward/backward Kolmogorov equations.

DefinitionChemical Master Equation

Consider $N$ chemical species with copy-number vector $\mathbf{n}\in\mathbb{Z}_{\ge0}^N$ undergoing $M$ reactions. Reaction $r$ changes the state by stoichiometric vector $\boldsymbol{\nu}_r$ and fires with propensity $a_r(\mathbf{n})$ , where $a_r(\mathbf{n})\,dt$ is the probability of one firing in $[t,t+dt)$ . The probability $P(\mathbf{n},t)$ obeys the Chemical Master Equation (CME):

3

\quad \frac{\partial P(\mathbf{n},t)}{\partial t}=\sum_{r=1}^{M}\Big[a_r(\mathbf{n}-\boldsymbol{\nu}_r)\,P(\mathbf{n}-\boldsymbol{\nu}_r,t)-a_r(\mathbf{n})\,P(\mathbf{n},t)\Big].

Derivation hint: Balance probability flux into state $\mathbf{n}$ (from $\mathbf{n}-\boldsymbol{\nu}_r$ via reaction $r$ ) against flux out of $\mathbf{n}$ . This is a gain–loss bookkeeping over an infinitesimal interval. For mass-action kinetics in volume $\Omega$ , a bimolecular reaction $S_i+S_j\to\dots$ has propensity $a_r=\frac{c_r}{\Omega}n_i n_j$ , while a unimolecular reaction has $a_r=c_r n_i$ , where $c_r$ is the stochastic rate constant.

TheoremSystem-Size Expansion and the Linear Noise Approximation

Write $\mathbf{n}=\Omega\,\boldsymbol{\phi}(t)+\sqrt{\Omega}\,\boldsymbol{\xi}$ , where $\boldsymbol{\phi}$ is the macroscopic concentration. Van Kampen's expansion of the CME in powers of $\Omega^{-1/2}$ yields, at leading order, the deterministic rate equation

4

\quad \frac{d\boldsymbol{\phi}}{dt}=\sum_r \boldsymbol{\nu}_r\, \tilde a_r(\boldsymbol{\phi}),

Derivation hint: Substitute the ansatz into the CME, Taylor-expand the step operators $E^{\pm\nu}\approx 1\pm\Omega^{-1/2}\partial_\xi+\dots$ , and collect the $\Omega^{1/2}$ terms. and at next order a linear Fokker–Planck (Ornstein–Uhlenbeck) equation for $\boldsymbol\xi$ with drift matrix (Jacobian) $A_{ij}=\partial_{\phi_j}\!\sum_r\nu_{ri}\tilde a_r$ and diffusion matrix $B=\sum_r \boldsymbol{\nu}_r\boldsymbol{\nu}_r^{\!\top}\tilde a_r(\boldsymbol\phi)$ . The stationary covariance $\Sigma$ solves the Lyapunov equation

5

\quad A\Sigma+\Sigma A^{\top}+B=0.

Derivation hint: Collect the $\Omega^0$ terms to get a Fokker–Planck equation with linear drift $A\xi$ and constant diffusion $B$ ; its Gaussian stationary solution has covariance solving the continuous Lyapunov equation.

ProofFrom CME to the Chemical Langevin Equation

We sketch Gillespie's derivation of the diffusion approximation. 1. Counting reactions. Let $K_r(t,t+\tau)$ be the number of firings of reaction $r$ in $[t,t+\tau)$ . Then exactly $\mathbf{n}(t+\tau)=\mathbf{n}(t)+\sum_r \boldsymbol{\nu}_r K_r$ . 2. Poisson structure. If $\tau$ is small enough that propensities are approximately constant, $K_r\sim\text{Poisson}\big(a_r(\mathbf{n})\tau\big)$ , independent across $r$ . This is the first timescale condition. 3. Macroscopic infinitesimal regime. If additionally $a_r(\mathbf{n})\tau\gg 1$ (many firings, the second condition — large copy numbers), the Poisson is well-approximated by a Gaussian: $K_r\approx a_r\tau+\sqrt{a_r\tau}\,\mathcal N_r(0,1)$ . 4. Langevin form. Substituting and writing $\sqrt{\tau}\,\mathcal N_r\to dW_r$ ,

6

\quad d\mathbf{n}=\sum_r \boldsymbol{\nu}_r\, a_r(\mathbf{n})\,dt+\sum_r \boldsymbol{\nu}_r\sqrt{a_r(\mathbf{n})}\,dW_r.

Derivation hint: Replace the exact integer Poisson increments by Gaussian increments with matched mean and variance; the $\sqrt{a_r}$ multiplicative noise is the hallmark of intrinsic fluctuations. 5. Consistency. Equation [6] is the Chemical Langevin Equation; its Fokker–Planck equation has drift $\sum_r\boldsymbol\nu_r a_r$ and diffusion $\sum_r\boldsymbol\nu_r\boldsymbol\nu_r^\top a_r$ , matching the LNA when linearized. The two conditions of steps 2–3 carve out the regime where both deterministic ( $\Omega\to\infty$ ) and discrete-stochastic descriptions reconcile. $\blacksquare$

∎

MD

Max Delbrück

1906–1981

Trained as a physicist under Bohr, Delbrück brought the language of fluctuations into molecular biology. His 1940 paper "Statistical Fluctuations in Autocatalytic Reactions" was among the first to treat a biochemical reaction as a stochastic birth process and ask how molecular noise shapes outcomes. With Salvador Luria, the 1943 Luria–Delbrück experiment used the distribution of mutant counts (not just the mean) to prove that bacterial mutations arise spontaneously before selection — a landmark demonstration that the variance of a stochastic process carries decisive biological information. He shared the 1969 Nobel Prize in Physiology or Medicine.

🧠 Intuition

Think of a single gene that switches randomly between OFF and ON, producing mRNA only when ON, with mRNA degrading at a fixed rate. Deterministically, you would predict a single steady mRNA level. Stochastically, three things happen. First, intrinsic noise: even at fixed promoter state, mRNA arrives in random bursts and decays in random steps, so the count fluctuates with Fano factor $\ge 1$ . Second, promoter switching (extrinsic-like) noise: slow ON/OFF transitions imprint long-lived correlations, producing bursty expression and a heavy-tailed, sometimes bimodal, distribution. Third, timescale separation matters: if switching is fast relative to mRNA decay, the gene looks "averaged" (Poisson-like, Fano $\approx 1$ ); if switching is slow, the population splits into high- and low-expressing subpopulations even though every cell is genetically identical. The single most important number is the Fano factor $F=\mathrm{Var}(n)/\langle n\rangle$ : $F=1$ flags a pure Poisson (constitutive) source, $F>1$ flags bursting, and $F<1$ flags negative feedback or sub-Poissonian control.

Figure

Gillespie trajectory (stochastic, jagged) versus deterministic mean for a birth–death process; fluctuations of relative size 1/√N around the steady state.

```svg

```

ExampleBirth–Death Process: mRNA Statistics

Given: mRNA is produced at constant rate $k=20\ \text{min}^{-1}$ and degraded at per-molecule rate $\gamma=0.5\ \text{min}^{-1}$ . Reactions: $\varnothing\xrightarrow{k} m$ , $m\xrightarrow{\gamma} \varnothing$ , propensities $a_1=k$ , $a_2=\gamma n$ .

Step 1 — Mean equation. $\frac{d\langle n\rangle}{dt}=k-\gamma\langle n\rangle \Rightarrow \langle n\rangle_{ss}=k/\gamma = 20/0.5 = 40$ .

Step 2 — Stationary distribution. The CME stationary solution is Poisson with parameter $\lambda=k/\gamma$ :

7

\quad P_{ss}(n)=\frac{\lambda^n e^{-\lambda}}{n!},\qquad \lambda=\frac{k}{\gamma}.

Derivation hint: Detailed balance for the birth–death chain gives $P(n+1)/P(n)=a_1/a_2(n+1)=k/[\gamma(n+1)]$ ; iterating yields the Poisson form.

Step 3 — Variance and Fano factor. Poisson $\Rightarrow \mathrm{Var}(n)=\lambda=40$ , so $F=\mathrm{Var}/\langle n\rangle=1$ .

Step 4 — Coefficient of variation. $\mathrm{CV}=\sigma/\langle n\rangle = 1/\sqrt{\lambda}=1/\sqrt{40}\approx 0.158$ .

\boxed{\langle n\rangle_{ss}=40\ \text{molecules},\quad \mathrm{Var}=40,\quad F=1,\quad \mathrm{CV}\approx 15.8\%}

ExampleWright–Fisher / Moran Fixation Probability

Given: A haploid population of $N=1000$ individuals; a single new neutral mutant ( $i=1$ ) appears. We want its probability of eventual fixation, and the expected time, under the neutral Moran model (one birth and one death per step, sampling proportional to type frequency).

Step 1 — Fixation probability (neutral). By the martingale property of allele frequency $X_t$ (it is a bounded martingale; $\mathbb E[X_\infty]=X_0$ ), fixation probability equals initial frequency:

8

\quad \pi_{\text{fix}}(i)=\frac{i}{N}.

Derivation hint: $X_t=i/N$ is a martingale under neutrality; optional stopping at the absorbing boundaries $\{0,1\}$ gives $\mathbb E[X_\infty]=\pi_{\text{fix}}\cdot 1 + (1-\pi_{\text{fix}})\cdot 0 = X_0$ . For one mutant: $\pi_{\text{fix}}=1/N = 1/1000 = 0.001$ .

Step 2 — Selection correction. With selection coefficient $s$ (Moran with relative fitness $1+s$ ), the diffusion approximation gives

9

\quad \pi_{\text{fix}}(p)=\frac{1-e^{-2Nse\,p}}{1-e^{-2Nse}},

with $N_e$ the effective size; for one beneficial mutant of small $s$ , $\pi_{\text{fix}}\approx \frac{1-e^{-2s}}{1-e^{-2Ns}}\approx 2s$ when $1\ll 2Ns$ . If $s=0.01$ : $\pi_{\text{fix}}\approx 2s = 0.02$ — twenty-fold above neutral.

Step 3 — Expected neutral fixation time (for a mutation that does fix) scales as $\approx 2N$ generations (Moran) or $\approx 4N$ (diploid Wright–Fisher).

\boxed{\pi_{\text{neutral}}=0.001;\quad \pi_{\text{sel}}(s=0.01)\approx 0.02;\quad T_{\text{fix}}\sim \mathcal O(N)\ \text{generations}}

⚠️ Common mistake

"Adding Gaussian white noise to the ODE always gives the right stochastic model." False in general. Intrinsic biochemical noise is multiplicative and state-dependent ( $\sqrt{a_r(\mathbf n)}$ in [6]), not additive constant noise. Moreover the diffusion approximation requires large copy numbers; near extinction ( $n\sim\mathcal O(1)$ ) the Gaussian Langevin breaks down because it allows negative and non-integer counts and badly mis-estimates rare events like extinction times. Use the exact CME / Gillespie there. Also: Itô vs. Stratonovich matters — the CME-derived Langevin equation is unambiguously Itô because reaction firings are forward-in-time Poisson events.

Counterexample

Noise that creates an extra "phenotype" with no deterministic counterpart. Consider a positive-feedback gene with a single stable deterministic fixed point. One might assume the stationary distribution is unimodal around that point. But if the feedback is sharp and copy numbers are low, the stochastic system can be bimodal — a substantial fraction of cells sit at a low-expression state never visited by the ODE. Here the macroscopic potential and the stochastic potential (from $-\ln P_{ss}$ ) have qualitatively different shapes: noise reshapes the landscape, an effect the LNA (which is Gaussian and unimodal by construction) cannot capture. This shows deterministic stability analysis can give the wrong number of attractors.

connection

Stochastic processes in biology sit at the confluence of several mathematical fields. The CME is a Continuous-Time Markov Chain on $\mathbb Z_{\ge0}^N$ ; its diffusion limit is governed by the Fokker–Planck Equation and equivalently by Stochastic Differential Equations interpreted in the Itô sense. The system-size expansion is a Singular Perturbation Theory in $\Omega^{-1/2}$ . Population-genetics models connect to diffusion processes on bounded intervals (Wright–Fisher diffusion, Kimura's equation) and to Coalescent Theory when run backward in time. Neural spike trains and stochastic ion channels relate to point processes and Hidden Markov Models. Extinction and rare events tie into Large Deviation Theory (the WKB / eikonal method for the CME) and First-Passage Time analysis.

🔬 Where this appears

⏳ History

- 1900–1906 — Andrey Markov introduces Markov chains; Bachelier and later Einstein/Smoluchowski (1905–06) develop diffusion and Brownian motion, the analytic backbone.

- 1924 — G. Udny Yule models species evolution as a pure-birth (Yule) process — one of the first explicitly biological stochastic processes.

- 1930–1932 — Sewall Wright and R. A. Fisher formalize random genetic drift; the Wright–Fisher model of allele frequency change is born; Kolmogorov's forward/backward equations (1931) provide the analytic engine.

- 1940 — Max Delbrück, Statistical Fluctuations in Autocatalytic Reactions — molecular noise enters biochemistry.

- 1943 — Luria–Delbrück fluctuation test uses the variance of a stochastic process to demonstrate pre-selective mutation.

- 1955 — Motoo Kimura derives the diffusion (Kimura) equation for allele frequencies; later the neutral theory (1968).

- 1958 — P. A. P. Moran introduces the Moran model, a continuous-time alternative to Wright–Fisher with exact fixation results.

- 1961 — N. T. J. Bailey and others develop stochastic epidemic theory; the chemical master equation is formalized (McQuarrie 1967).

- 1976–1977 — Daniel Gillespie publishes the exact Stochastic Simulation Algorithm, making CME dynamics computable; van Kampen's system-size expansion (1961, textbook 1981) systematizes the noise expansion.

- 2002–2008 — Elowitz's dual-reporter experiment distinguishes intrinsic from extrinsic noise; explosion of single-cell stochastic gene-expression studies.

open-problem

Several deep challenges remain. (i) Rare events and extinction: the Gaussian LNA fails for tails; computing extinction times requires WKB/large-deviation methods whose pre-exponential factors and validity in multidimensional, oscillatory systems are still being worked out. (ii) Spatial stochasticity: the well-mixed CME ignores diffusion-limited reactions; reaction–diffusion master equations (RDME) suffer mesh-dependence pathologies for bimolecular reactions, and rigorous convergent particle-based methods remain computationally heavy. (iii) Moment closure: nonlinear propensities generate an infinite, unclosed hierarchy of moment equations; no universally accurate closure exists, and closures can produce spurious negative variances. (iv) Inference: estimating CME parameters from snapshot single-cell data (likelihood is generally intractable) drives ongoing work in finite-state-projection methods, neural likelihood approximations, and optimal-transport-based trajectory inference. (v) Multiscale coupling between fast molecular noise and slow phenotypic/evolutionary dynamics lacks a fully rigorous, general averaging framework.

Formal Statement

Why It Holds / Derivation

Building Intuition

Worked Examples

Common Misconceptions

Connections

Real-World Applications

Historical Development

Subtleties & Limits

See Also