Stochastic Processes in Biology

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–).

A stochastic process in biology is a family ${X (t)}_{t \in T}$ of random variables indexed by time, representing a fluctuating biological quantity — copy numbers of mRNA, allele frequencies, neuron membrane potential, or population sizes. These processes matter because biology operates at scales where discreteness and randomness dominate: a bacterium may contain only a handful of transcription-factor molecules, so the law of large numbers does not apply and deterministic ODEs misrepresent the dynamics. The central machinery is the theory of Markov processes — continuous-time jump processes governed by master equations, and their diffusion limits governed by Fokker–Planck and stochastic differential equations. A surprising fact: noise is not merely a nuisance but a functional resource — stochastic gene expression generates phenotypic diversity in clonal populations (bet-hedging), and molecular fluctuations can drive cells across decision boundaries that deterministic systems would never cross.

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/ 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.

Formal Statement

Markov Process and Transition Semigroup

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

P (X (t) \in A ∣ X (t_{n}), \dots, X (t_{1})) = P (X (t) \in A ∣ X (t_{n})) .

⟳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} \geq 0$ ( $i \neq = j$ ) and $q_{ii} = - \sum_{j \neq = i} q_{ij}$ . The transition probabilities $P_{ij} (t) = P (X (t) = j ∣ X (0) = i)$ satisfy

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

⟳Combine the Chapman–Kolmogorov equation

P (t + s) = P (t) P (s)

with

P (h) \approx I + h Q

for small

h

, then take

h \to 0

to get the forward/backward Kolmogorov equations.

Chemical Master Equation

Consider $N$ chemical species with copy-number vector $n \in Z_{\geq 0}^{N}$ undergoing $M$ reactions. Reaction $r$ changes the state by stoichiometric vector $ν_{r}$ and fires with propensity $a_{r} (n)$ , where $a_{r} (n) d t$ is the probability of one firing in $[t, t + d t)$ . The probability $P (n, t)$ obeys the Chemical Master Equation (CME):

\frac{\partial P ( n , t )}{\partial t} = r = 1 \sum M [a_{r} (n - ν_{r}) P (n - ν_{r}, t) - a_{r} (n) P (n, t)] .

⟳Balance probability flux into state

n

(from

n - ν_{r}

via reaction

r

) against flux out of

n

. This is a gain–loss bookkeeping over an infinitesimal interval.

For mass-action kinetics in volume $Ω$ , a bimolecular reaction $S_{i} + S_{j} \to \dots$ has propensity $a_{r} = \frac{c _{r}}{Ω} n_{i} n_{j}$ , while a unimolecular reaction has $a_{r} = c_{r} n_{i}$ , where $c_{r}$ is the stochastic rate constant.

System-Size Expansion and the Linear Noise Approximation

Write $n = Ω ϕ (t) + Ω ξ$ , where $ϕ$ is the macroscopic concentration. Van Kampen's expansion of the CME in powers of $Ω^{- 1/2}$ yields, at leading order, the deterministic rate equation

\frac{d ϕ}{d t} = r \sum ν_{r} \tilde{a}_{r} (ϕ),

⟳Substitute the ansatz into the CME, Taylor-expand the step operators

E^{\pm ν} \approx 1 \pm Ω^{- 1/2} \partial_{ξ} + \dots

, and collect the

Ω^{1/2}

terms.

and at next order a linear Fokker–Planck (Ornstein–Uhlenbeck) equation for $ξ$ with drift matrix (Jacobian) $A_{ij} = \partial_{ϕ_{j}} \sum_{r} ν_{r i} \tilde{a}_{r}$ and diffusion matrix $B = \sum_{r} ν_{r} ν_{r}^{⊤} \tilde{a}_{r} (ϕ)$ . The stationary covariance $Σ$ solves the Lyapunov equation

A Σ + Σ A^{⊤} + B = 0.

⟳Collect the

Ω^{0}

terms to get a Fokker–Planck equation with linear drift

A ξ

and constant diffusion

B

; its Gaussian stationary solution has covariance solving the continuous Lyapunov equation.

Why It Holds / Derivation

From CME to the Chemical Langevin Equation

We sketch Gillespie's derivation of the diffusion approximation.

Counting reactions. Let $K_{r} (t, t + τ)$ be the number of firings of reaction $r$ in $[t, t + τ)$ . Then exactly $n (t + τ) = n (t) + \sum_{r} ν_{r} K_{r}$ .

Poisson structure. If $τ$ is small enough that propensities are approximately constant, $K_{r} \sim Poisson (a_{r} (n) τ)$ , independent across $r$ . This is the first timescale condition.

Macroscopic infinitesimal regime. If additionally $a_{r} (n) τ ≫ 1$ (many firings, the second condition — large copy numbers), the Poisson is well-approximated by a Gaussian: $K_{r} \approx a_{r} τ + a_{r} τ N_{r} (0, 1)$ .

Langevin form. Substituting and writing $τ N_{r} \to d W_{r}$ ,

d n = r \sum ν_{r} a_{r} (n) d t + r \sum ν_{r} a_{r} (n) d W_{r} .

⟳Replace the exact integer Poisson increments by Gaussian increments with matched mean and variance; the

a_{r}

multiplicative noise is the hallmark of intrinsic fluctuations.

Consistency. Equation [6] is the Chemical Langevin Equation; its Fokker–Planck equation has drift $\sum_{r} ν_{r} a_{r}$ and diffusion $\sum_{r} ν_{r} ν_{r}^{⊤} a_{r}$ , matching the LNA when linearized. The two conditions of steps 2–3 carve out the regime where both deterministic ( $Ω \to \infty$ ) and discrete-stochastic descriptions reconcile. $■$

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.

Building 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 $\geq 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 = Var (n) / ⟨ n ⟩$ : $F = 1$ flags a pure Poisson (constitutive) source, $F > 1$ flags bursting, and $F < 1$ flags negative feedback or sub-Poissonian control.

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

<svg viewBox="0 0 640 360" xmlns="http://www.w3.org/2000/svg">
<rect width="640" height="360" fill="#0a0e17"/>
<line x1="60" y1="320" x2="600" y2="320" stroke="#94a3b8" stroke-width="1.5"/>
<line x1="60" y1="40" x2="60" y2="320" stroke="#94a3b8" stroke-width="1.5"/>
<text x="330" y="350" fill="#e2e8f0" font-size="14" text-anchor="middle">time t</text>
<text x="20" y="180" fill="#e2e8f0" font-size="14" transform="rotate(-90 20 180)" text-anchor="middle">copy number n(t)</text>
<line x1="60" y1="150" x2="600" y2="150" stroke="#d4a853" stroke-width="2" stroke-dasharray="6 4"/>
<text x="470" y="142" fill="#d4a853" font-size="12">deterministic steady state n* = k/γ</text>
<polyline fill="none" stroke="#e2e8f0" stroke-width="1.6" points="60,300 75,300 75,280 95,280 95,260 110,260 110,240 130,240 130,200 150,200 150,180 170,180 170,160 195,160 195,170 215,170 215,150 240,150 240,130 265,130 265,160 290,160 290,140 315,140 315,170 340,170 340,150 365,150 365,135 390,135 390,165 415,165 415,145 440,145 440,170 465,170 465,150 490,150 490,135 515,135 515,160 540,160 540,145 565,145 565,158 590,158"/>
<text x="120" y="225" fill="#e2e8f0" font-size="11">transient relaxation</text>
<text x="400" y="115" fill="#e2e8f0" font-size="11">stationary fluctuations ~ √n*</text>
<rect x="430" y="130" width="160" height="45" fill="none" stroke="#94a3b8" stroke-width="0.8" stroke-dasharray="3 3"/>
</svg>

Worked Examples

Birth–Death Process: mRNA Statistics

Given: mRNA is produced at constant rate

k = 20 min^{- 1}

and degraded at per-molecule rate

γ = 0.5 min^{- 1}

. Reactions:

\emptyset k m

m γ \emptyset

, propensities

a_{1} = k

a_{2} = γ n

. Step 1 — Mean equation.

\frac{d ⟨ n ⟩}{d t} = k - γ ⟨ n ⟩ \Rightarrow ⟨ n ⟩_{ss} = k / γ = 20/0.5 = 40

. Step 2 — Stationary distribution. The CME stationary solution is Poisson with parameter

λ = k / γ

P_{ss} (n) = \frac{λ ^{n} e ^{- λ}}{n !}, λ = \frac{k}{γ} .

⟳Detailed balance for the birth–death chain gives

P (n + 1) / P (n) = a_{1} / a_{2} (n + 1) = k / [γ (n + 1)]

; iterating yields the Poisson form.

Step 3 — Variance and Fano factor. Poisson

\Rightarrow Var (n) = λ = 40

, so

F = Var / ⟨ n ⟩ = 1

. Step 4 — Coefficient of variation.

CV = σ / ⟨ n ⟩ = 1/ λ = 1/ 40 \approx 0.158

⟨ n ⟩_{ss} = 40 molecules, Var = 40, F = 1, CV \approx 15.8%

Wright–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;

E [X_{\infty}] = X_{0}

), fixation probability equals initial frequency:

π_{fix} (i) = \frac{i}{N} .

⟳

X_{t} = i / N

is a martingale under neutrality; optional stopping at the absorbing boundaries

{0, 1}

gives

E [X_{\infty}] = π_{fix} \cdot 1 + (1 - π_{fix}) \cdot 0 = X_{0}

For one mutant: $π_{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

π_{fix} (p) = \frac{1 - e ^{- 2 N se p}}{1 - e ^{- 2 N se}},

with

N_{e}

the effective size; for one beneficial mutant of small

s

π_{fix} \approx \frac{1 - e ^{- 2 s}}{1 - e ^{- 2 N s}} \approx 2 s

when

1 ≪ 2 N s

If $s = 0.01$ : $π_{fix} \approx 2 s = 0.02$ — twenty-fold above neutral.

Step 3 — Expected neutral fixation time (for a mutation that does fix) scales as

\approx 2 N

generations (Moran) or

\approx 4 N

(diploid Wright–Fisher).

π_{neutral} = 0.001; π_{sel} (s = 0.01) \approx 0.02; T_{fix} \sim O (N) generations

Common Misconceptions

"Adding Gaussian white noise to the ODE always gives the right stochastic model." False in general. Intrinsic biochemical noise is multiplicative and state-dependent (

a_{r} (n)

in [6]), not additive constant noise. Moreover the diffusion approximation requires large copy numbers; near extinction (

n \sim 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. 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.

Connections

Stochastic processes in biology sit at the confluence of several mathematical fields. The CME is a Continuous-Time Markov Chain on $Z_{\geq 0}^{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 $Ω^{- 1/2}$ . Population-genetics models connect to Diffusion Process|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 Process|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.

Real-World Applications

Single-cell gene expression / single-molecule mRNA counting (smFISH, scRNA-seq). The two-state (telegraph) model and its Fano factor decode burst size and frequency from measured distributions — foundational in systems biology (Elowitz, Raj, van Oudenaarden, ~2002–2008).
Ion channel gating (Hodgkin–Huxley with channel noise). Channels are Markov chains over conformational states; the Gillespie / Langevin treatment predicts spontaneous firing and variability in small neurons (single-channel patch-clamp, Neher–Sakmann).
Population genetics & molecular evolution. Wright–Fisher diffusion and the Kimura equation underlie estimates of effective population size, selection coefficients, and the neutral theory; the coalescent powers modern phylogenetic and demographic inference.
Epidemic dynamics. Stochastic SIR / SIS models give the probability of outbreak extinction, the critical community size for measles persistence, and final-size distributions — essential when case numbers are small (early outbreak, eradication endgame).
Bacterial persistence and bet-hedging. Stochastic switching between growth states explains antibiotic persisters; the rate of phenotype switching is tuned to environmental fluctuation statistics (Balaban et al.).
Enzyme kinetics at single-molecule resolution. Stochastic Michaelis–Menten and the "molecular memory" of turnover-time distributions (English, Xie et al.) reveal dynamic disorder invisible to bulk assays.
Synthetic biology / circuit design. Noise propagation analysis (using the Lyapunov equation [5]) guides design of robust toggle switches, oscillators (repressilator), and feedback controllers that suppress or exploit fluctuations.

Historical Development

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.

Subtleties & Limits

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.