Stochastic dynamics: the master equation and Gillespie

A molecular reaction is an event. Give every event an integer jump and a state-dependent clock, then derive the probability equation, exact sample paths, and the deterministic limit from the same object. CCBS · EST 5053 · Fall 2026 · lecture 5 · core page (full) · 2026-09-15 (Tue) 10:40–12:15 · students present, lecturer-backed Take a result we can already prove with Lecture 4. Molecules are added at a constant rate

b

and removed independently at per-molecule rate

d>0

. If

m(t)

is a smooth count, its equation is

\dot m=b-dm

. The only fixed point is

m_*=b/d

, and the slope

-d

is negative, so deviations relax. We have established stability. Does that tell us how much one cell fluctuates, how long a deviation lasts, or whether it ever contains zero molecules?

Lecture 4 also gave us stable switches and oscillations as possible behaviors. Can we build those behaviors into a cell? The toggle switch and repressilator, published together in 2000, made that question a landmark of the emerging synthetic-biology program: construct a small genetic circuit, predict its dynamics, and measure whether a living cell realizes them.²⁰¹ The comparison will open this lecture. We will then develop the stochastic tools on the simpler addition-removal system, where each new prediction can be checked.

The claim

At molecular scale, a reaction network is an integer state moved by random events. State-change vectors and event clocks define the process. The chemical master equation (CME) evolves probabilities over those states. There is no general closed-form solution, but the Gillespie direct method samples histories without enumerating the state space. Its exactness is relative to a well-mixed, nonexplosive Markov jump model with constant hazards between events. The biological observable decides whether a path, a distribution, or the deterministic limit is the useful description.

0–8 min: Gap. Build an oscillator, then ask what a single-cell history adds.
8–14 min: Kinds. Separate hidden state, reaction variation and readout noise.
14–22 min: Scale. Estimate copies, derive Poisson, then count packets.
22–30 min: Map. Define the local volume, exchange criterion and turnover question.
30–39 min: Object. Write the integer state, jump, propensity and units.
39–48 min: Probability. Derive the chemical master equation.
48–52 min: Boundary. Make zero legal and conserve probability.
52–64 min: Calibration. Poisson, the observable generator and exact moments.
64–75 min: Sampling. Two uniforms sample the next event.
75–84 min: Estimation. Independent paths, time weights, uncertainty and extinction.
84–89 min: Views. Keep paths, ensembles and the deterministic limit separate.
89–92 min: Contract. State the world for which the direct method is exact.
92–95 min: Claim. Specify an observable, estimator and falsifier.

Selection and drop order

The page retains thirteen points. Sections 2 and 4 are branches around an eleven-point dependency spine. The 95-minute route reads three of the seven map examples and keeps the spine definition and exchange criterion. If time runs short, leave the umbrella matrix, detailed geometric estimate, identical-pair derivation and plotted system-size comparisons for reading. Keep the CME, generator setup, absorbing zero, calibration, hand-traced event and the two estimation procedures. Movies, playgrounds and addenda stay outside the clock.

Figure 1. A stable level still admits fluctuating histories. This is the addition-removal model, not repressilator data. All views use

b=4\,\mathrm{min}^{-1}

d=0.2\,\mathrm{min}^{-1}

, and initial count zero. The smooth curve is the exact expected count, each staircase a simulated history, and the bars the analytic count distribution at 30 minutes. The equations below derive their common mechanism.

Part 1

Say what you are looking at

The word noise names four different things. Separate them, then count the events behind the observable and decide whether that number is small.

1A smooth trajectory is not a single cell

A deterministic model can be correct for a concentration and unable to answer the path-level question.

The repressilator asks a cell to keep time. LacI represses TetR, TetR represses λ cI, and λ cI represses LacI. If one repressor rises, the next falls, releasing the third, which eventually suppresses the first. These are regulatory links, not chemical conversion arrows. Synthesis and removal of the three proteins are composite reaction channels whose effective laws include repression. A separate GFP reporter reads out the circuit.

Opening experiment. Designed oscillations and an actual cell. Excerpts of Elowitz and Leibler's original Figures 1c and 2a–c: deterministic and stochastic simulations at left; microcolony images and the GFP trace of a tracked lineage at right. Simulation axes count repressor proteins; the measured axis is reporter fluorescence in arbitrary units. Septation marks are shown below the measured trace. The comparison is qualitative, not a parameter-matched fit of the same observable. Insets on the simulations show autocorrelation.¹

The experimental oscillator has a history, not just a period. Peaks vary in timing and height, even along a lineage. The original paper already compared deterministic and stochastic simulations. It did not establish that every irregularity was intrinsic reaction noise. This is the question we inherit: which part of the observed variability comes from reaction events, shared cell conditions, or the measurement?

Watch the original microcolony (optional, 20 seconds)

Original repressilator microcolony movie, provided by the authors through Biological Circuit Design. Brightness reports GFP, not a direct count of all three repressors. Compare a tracked lineage with the whole colony; loss of population synchrony is not the same observable as an irregular period in one cell. A later redesign is discussed in the optional follow-up.

The circuit-building work helped make such single-cell questions experimentally concrete. Stochastic gene-expression and cell-fate models already existed, including the phage decision studied by Arkin and colleagues in 1998. The dual-reporter experiments that followed supplied ways to separate sources of variation.¹⁵⁷ Our first task is to choose what a model should predict.

Suppose half the cells contain no transcript and half contain twenty. Their mean is ten, although no cell is near ten. Or suppose every cell eventually reaches the same level, but the time of the first transcript determines whether a downstream switch fires. In both cases a mean trajectory suppresses the observation that carries the biology.

There are three questions to keep separate from the first minute:

object	question it answers	example observable
one count history $N_X(t)$ (defined in Section 5)	what happened to this cell, and when?	first-passage time, pulse duration, event order
a distribution $p(n,t)$	how likely is each possible count state $n$ ?	probability of zero, tail probability, fraction above threshold
a smooth concentration $c(t)$	what macroscopic trajectory emerges?	concentration level or deterministic attractor

A stochastic model is warranted when the observable depends on an integer boundary, a distribution, an event time, or an individual history. The word cell is not enough. A millimolar metabolite in a bacterium may be well described by concentration, while one promoter copy remains discrete in a much larger cell.

Room check, 60 seconds. A population-average assay reports ten transcripts per cell. List two single-cell distributions with that mean and different probabilities of zero. Which downstream behavior could distinguish them?

Name the observable before choosing the model. A path, a distribution and a smooth concentration answer three different questions, and a mean trajectory is silent about the first two whatever its parameters are.

2Not every spread is noise

A wide histogram is an observation. Four different things produce one, and only one of them is what a stochastic reaction model describes.

Before writing a probability model, ask two questions in order. Is the spread random at all? If it is, where does the randomness enter?

A hidden state can make a deterministic rule look random. Here is my umbrella rule: when I leave home or the office, I carry the umbrella if it is raining and the umbrella is where I am. Otherwise I leave it where it is. Given the rain and both locations, the decision is completely determined. Someone recording only “carried it / did not carry it” sees an irregular sequence because that record omits the state that controls the decision.

What is the umbrella's state?

Let $L_t$ be my departure location and $U_t$ the umbrella's location, each home or office. Let $R_t\in\{0,1\}$ say whether it rains. My action is $B_t=\mathbf1_{\{L_t=U_t\}}R_t$ . My next location is the other place. The umbrella moves there exactly when $B_t=1$ . The pair $(L_t,U_t)$ stores the memory of previous decisions.

Umbrella at departure?	Weather	Carry it?	At the next departure?
yes	rain	yes	available
yes	dry	no	unavailable
no	either	no	available

If we additionally model successive rain observations as independent with probability $p_{\mathrm{rain}}$ , availability $A_t=\mathbf1_{\{L_t=U_t\}}$ is a two-state discrete-time Markov chain. With rows for current state 0,1 and columns for next state 0,1,

P=\begin{bmatrix}0&1\\1-p_{\mathrm{rain}}&p_{\mathrm{rain}}\end{bmatrix}.

The decision is deterministic conditional on the input; the weather model supplies the randomness. Correlated weather needs more state, and the carry/not-carry sequence alone need not be Markov. This transition matrix contains probabilities per trip, unlike the continuous-time rate matrix introduced in Section 6.

A biological dataset can likewise mix cells with different plasmid counts, cell-cycle positions, or environments. Conditioning on those variables can explain some of the spread. It need not remove all spread or split a histogram into sharp groups. That is why we should name the state before assigning every unexplained fluctuation to a reaction.

A count model also leaves molecular encounters unresolved. Two cells with the same current molecule counts need not have the same next reaction time. Positions, solvent motion and encounter histories are absent from that state. Stochastic reaction kinetics represents their effect through conditional event probabilities. Whether counts alone give memoryless clocks is a physical modelling assumption, not a consequence of merely failing to measure something.

Randomness is relative to the description

A classical molecular-dynamics model can track trajectories while a coarser diffusion model assigns probabilities. The descriptions retain different information. Coarse-graining may produce effective randomness and memory, so a Markov approximation must be justified by mixing and timescale separation. Enlarging the state can expose hidden structure; it does not promise to eliminate every kind of uncertainty.

Genuine randomness still arrives by two different routes, and one experiment separates them. Put two copies of the same reporter in one cell and give them different colours. Whatever is shared by the cell, its size, its ribosome pool, its stage in the cycle, moves both readings together. Whatever belongs to one copy's own reaction history moves them apart. Writing $c$ and $y$ for the two normalised readings across a population,

\eta_{\mathrm{int}}^{2}=\frac{\langle(c-y)^{2}\rangle}{2\langle c\rangle\langle y\rangle},\qquad \eta_{\mathrm{ext}}^{2}=\frac{\langle cy\rangle-\langle c\rangle\langle y\rangle}{\langle c\rangle\langle y\rangle}.

(1)

The first is intrinsic noise, the events of this copy of the machinery, and it is what the rest of this lecture models. The second is extrinsic noise, shared cell state, and it is outside the model unless that state is put into it.⁷⁸ The decomposition is a definition relative to the chosen reporters, not a partition of nature.

The instrument contributes an additional spread. Molecule counting gives a discrete observable, but discreteness does not imply Poisson statistics. Fluorescence is an observation model involving abundance, maturation, background and measurement error. Some variability in these factors is biological, not instrument noise. Counting assays also require calibration and have errors.

A product of positive factors can be approximately log-normal when the sum of their logarithms has an approximately Gaussian distribution. Furusawa and colleagues studied biological mechanisms for such abundance statistics, not a universal detector artifact.⁹ Gamma laws also fit many bacterial protein distributions.⁶ A fitted family is useful evidence only together with the count, cell state and measurement process it describes. The optional distribution map makes these distinctions explicit.

Figure 2. Distinguish state variation, event noise, and the readout. The left mixture is drawn as the sum of its two components. The scatter panels simulate 2,000 reporter pairs, once with a shared factor and once independently. The right panel compares an illustrative Poisson mass function and log-normal density with the same mean. These are example distributions, not experimental data or universal assay laws.

Say which of the four you are looking at before you model it. The chemical master equation describes the third, conditional on the state you declared, and the other three are still there when it is right.

3Count the events, not the molecules

The relative spread of a molecular pool is set by the number of independent events that made it. That number can be far smaller than the pool.

Two counts, in order. How many molecules, and how many events.

Molecules first. For molar concentration $c$ and volume $V$ in liters,

N=cVN_A.

(2)

Do it once for a bacterium and never look it up again. E. coli is about $1\,\mathrm{fL}=10^{-15}\,\mathrm{L}$ . One molecule is $1/N_A\simeq1.7\times10^{-24}$ mol, and dividing by the volume gives $1.7\times10^{-9}$ M. Round it, because the point is to carry it in your head:

\pf{\text{one copy per bacterium}\;\simeq\;1\,\mathrm{nM}}.

(3)

Lecture 2 uses the same rounded rule, and everything follows from it by multiplication. A regulator at 10 nM in a bacterium is about ten copies. A signalling protein at 1 µM is about a thousand. A metabolite at 1 mM is about a million. A mammalian cell is roughly two thousand times bigger, so the same 10 nM is about twenty thousand copies there.

Now the events. Suppose the observable is refreshed by events that occur at some total rate over a window $T$ . Cut the window into $m$ slots so short that no slot can hold two events. Each slot either fires or does not, independently, with the same probability $p$ . That is $m$ coin flips, so

\Pr\{K=k\}=\binom{m}{k}p^{k}(1-p)^{m-k}.

(4)

Hold the mean $\mu=mp$ fixed and let the slots shrink. Two lines finish it. The empty case is $(1-\mu/m)^{m}\to e^{-\mu}$ , and the ratio of neighbours is

\frac{\Pr\{K=k+1\}}{\Pr\{K=k\}}=\frac{m-k}{k+1}\cdot\frac{\mu}{m-\mu}\;\longrightarrow\;\frac{\mu}{k+1}.

Climbing that ladder from the empty case gives the whole distribution:

\Pr\{K=k\}=e^{-\mu}\frac{\mu^{k}}{k!}.

(5)

This Poisson count requires independent rare opportunities with a constant event rate over the window. Lecture 2's rapid mixing makes that a useful starting point for a constant source. It is not true for every reaction count over a fixed window: substrate depletion, promoter switching, or changes in the number of eligible molecules make the clock state dependent. Section 5 keeps that dependence rather than assuming a Poisson count for the whole history.

The moments come from the same binomial, before the limit. The mean is $mp=\mu$ . The variance is $mp(1-p)\to\mu$ . So the variance equals the mean, and

\mathrm{CV}=\frac{\sqrt{\operatorname{Var}(K)}}{\langle K\rangle}=\frac{\sqrt{\mu}}{\mu}=\pf{\mu^{-1/2}}.

(6)

A hundred events give ten per cent. Ten thousand give one per cent. One event gives order one.

Figure 3. The Poisson law, built rather than asserted. The middle panel evaluates

\mathrm{Binomial}(m,\mu/m)

\mu=3

for

m=6

and

m=30

against its limit. The right panel illustrates delivery timing conditional on exactly a thousand molecules, singly or in ten packets. The CV calculations below it concern the unconditioned Poisson-arrival ensemble with that mean, not variation of these fixed endpoints.

The count of events is not the count of molecules. Count arrivals over a window with negligible removal. Let each event deliver a fixed packet of $b_{\mathrm{pkt}}$ molecules, and let $K$ be Poisson with mean $\mu$ . Then $N=b_{\mathrm{pkt}}K$ , mean $\bar N=b_{\mathrm{pkt}}\mu$ , and variance $b_{\mathrm{pkt}}^2\mu=b_{\mathrm{pkt}}\bar N$ . The Fano factor is $b_{\mathrm{pkt}}$ and

\mathrm{CV}(N)=\frac{\sqrt{b_{\mathrm{pkt}}\bar N}}{\bar N}=\pf{\sqrt{b_{\mathrm{pkt}}/\bar N}}.

(7)

Compare pools with mean a thousand molecules. Single arrivals give $\mathrm{CV}=1/\sqrt{1000}\simeq0.03$ ; packets of a hundred give $\mathrm{CV}=\sqrt{100/1000}\simeq0.32$ . The second pool has only ten arrival events on average. Translation from short-lived transcripts supplies a biological version of packets, but the number translated per transcript is random. With removal, the stationary variance also depends on how packets are filtered. Those distinctions are worked out in the burst addendum.

Comparison with the exact conversion and with measurement

The exact constant is $N_A=6.02214076\times10^{23}\,\mathrm{mol}^{-1}$ . One molecule per femtoliter is 1.66 nM, and 1 µM in 1 fL is 602 molecules. The rounded mental rule is accurate within a factor of about 1.7, sufficient for these scale comparisons.

The bacterial proteome survey found approximately inverse-mean scaling of relative variance at low expression. That scaling is consistent with counting and bursting models, not proof that every protein is Poisson. A lower envelope observed in that dataset is not a theorem forbidding feedback from producing sub-Poisson fluctuations.⁶

For fixed packets arriving as a Poisson process over a window without appreciable removal, estimate the spread as $\sqrt{b_{\mathrm{pkt}}/\bar N}$ . Here $\bar N$ counts molecules and $\bar N/b_{\mathrm{pkt}}$ counts events. Identify the event size and observation window before applying the rule.

4Where noise matters, and where it does not

First specify the volume that contains the molecules. Then compare the reaction, exchange and observation times that determine how their fluctuations are seen.

Figure 4 separates a count conversion from a noise reference. Each horizontal position is the mean count $\bar N=cVN_A$ , with $c$ the mean concentration. The vertical position uses Section 3's independent-arrival reference. The dashed curve asks what fixed packets of thirty would do under that section's window assumptions. The seven named examples locate count scales. Their actual biological noise need not follow either curve.

A dendritic spine is a tiny synaptic protrusion on a neuron's input-receiving branch, or dendrite. A small head connects to the dendrite through a narrow neck. For the example, take the head volume to be $0.1\,\mu\mathrm m^3=0.1\,\mathrm{fL}=10^{-16}\,\mathrm L$ . This is the protrusion's volume. It is much smaller than the whole neuron.

The calcium label counts free ions in that head. Using $N_A\approx6\times10^{23}\,\mathrm{mol}^{-1}$ , 100 nM in this volume gives $\bar N\approx(10^{-7})(10^{-16})(6\times10^{23})=6$ ions. The figure's right-hand comparison instead fixes every concentration at 10 nM, giving about 0.6 ions in the spine. A fractional mean describes an average over repeated observations. Every instantaneous count is an integer. Buffer-bound calcium and calcium stored inside organelles are separate pools.

Figure 4. Count scales under a declared noise reference. Mean count is computed from concentration and volume. The solid and dashed curves are the single-arrival and fixed-packet reference laws from Section 3. They are not measured calcium noise. The calcium point uses 100 nM free calcium in a dendritic spine. The right panel uses 10 nM throughout and reports the corresponding mean counts and reference CVs. Exchange, buffering, turnover and the readout's averaging time require additional information.

The volume comparison holds the concentration fixed. At 10 nM the chosen volumes contain about 0.6 molecules in a spine, six in a bacterium, two hundred in a yeast cell and twelve thousand in a mammalian cell. Under the Poisson reference, larger counts give smaller relative fluctuations. This conclusion depends on the counting model as well as the volume. Cell type alone does not determine the observed noise.

When can the spine act as its own compartment?

A separate calcium signal requires exchange to be slow relative to local calcium handling. Let $\tau_{\mathrm{exchange}}$ describe equilibration through the neck and $\tau_{\mathrm{clear}}$ describe removal from the head's cytoplasm. Pumps can clear a calcium increase before much spreads into the dendrite when $\tau_{\mathrm{clear}}\ll\tau_{\mathrm{exchange}}$ . The head remains connected. Its internal mixing must also be fast enough for a single well-mixed state to describe the observable.

Calcium measurements support this separation under specific conditions. In rat hippocampal CA1 neurons, Sabatini and colleagues estimated clearance at about 12–15 ms. They inferred native equilibration times above a second after accounting for the calcium-binding indicator's effect on transport. That slower exchange is an inference, not a direct dye-free measurement or a property of every molecule in every spine. The exposition derives the geometric estimate and explains the evidence.³⁰

Turnover changes how many independent fluctuations a measurement can average. A small pool can replace its contents many times during a long observation. For constant addition and first-order removal, write $b$ for the addition rate and $d$ for each molecule's removal rate. Increasing both at fixed $b/d$ preserves the stationary count distribution but shortens its memory. The instantaneous spread stays the same, while a time average becomes more precise. Section 10 will turn that distinction into estimators and error formulas.

Estimate a protein burst before assigning a packet size. Suppose one mRNA is translated once every ten seconds and survives for a mean three minutes. It then makes about $180/10\simeq20$ proteins. More generally, if translation events occur at rate $k_{\mathrm{tl}}$ per transcript and mRNA removal at rate $d_m$ ,

\bar B_{\mathrm{tl}}=k_{\mathrm{tl}}\,\tau_m=\frac{k_{\mathrm{tl}}}{d_m},\qquad \tau_m=\frac{1}{d_m}=\frac{t_{1/2,m}}{\ln2}.

The mean burst size $\bar B_{\mathrm{tl}}$ is proteins per transcript. A faster translation initiation rate or longer transcript lifetime makes it larger. A transcriptional burst is different: it counts RNAs made in one promoter-on episode, with mean $k_{\mathrm{tx}}/\beta$ if that episode ends at constant rate $\beta$ . Neither burst size has one value shared by all genes and conditions.

Estimate meets measurement

Under weak lac expression in E. coli, Cai, Friedman and Xie inferred mean bursts of $5\pm2$ active β-galactosidase tetramers, or $20\pm8$ protein monomers. They used catalytic amplification and calibration, not direct movies of every translation event.²² The illustrative estimate reaches the right scale for this reporter. The packet of thirty in Figure 4 is a comparison line, not a universal measured burst size.

Three features require additional counting information. A promoter locus has only a few copies even in a large cell. Bursts can make many molecules share one initiating event. A small local volume can have its own dynamics when exchange is sufficiently slow. Each feature can matter in either a bacterium or a eukaryotic cell.

Where measured proteins sit

Taniguchi and colleagues measured abundances from roughly $10^{-1}$ to $10^4$ copies per cell. At higher expression the relative variance approached a floor near $0.1$ , corresponding to CV of order 30%. Shared cellular variation limited how much increasing abundance reduced the observed spread. This is a result under the survey's conditions, not a law for every abundant protein.⁶ Newman and colleagues found the same shape in budding yeast, and also found that noise sorts by what a protein does: proteins that answer to the environment are noisy, and the protein-synthesis machinery is quiet.²⁸

Room check, 60 seconds. A signaling protein has mean concentration 100 nM. Estimate its mean count and Poisson-reference CV in the three chosen volumes: a bacterium, a mammalian cell and a dendritic spine. Which transport and turnover times would you need before interpreting those CVs as fluctuations in a real measurement? If molecules arrive in packets, also specify the packet law and observation window.

Specify the molecule, volume and averaging window before estimating noise. Use the count conversion to locate discreteness, an exchange-versus-reaction comparison to justify a local compartment, and a turnover model to decide how much independent information the observation contains.

Part 2

Build the object, and meet the infinity

Start with Lecture 3's reaction list. Its jumps and clocks now define a probability-conserving dynamical system: the chemical master equation.

5One reaction list becomes integer jumps and clocks

The stochastic model is specified before the random numbers appear.

Keep Lecture 3's reaction accounting. For an elementary reaction among species $X_1,\ldots,X_s$ , the non-negative integer coefficients $\alpha_{ir}$ and $\beta_{ir}$ count reactants and products. Their vectors are $\alpha_r$ and $\beta_r$ , and one event changes the count by $\gamma_r=\beta_r-\alpha_r$ :

\sum_{i=1}^{s}\alpha_{ir}X_i\xrightarrow{k_r}\sum_{i=1}^{s}\beta_{ir}X_i, \qquad v_r(c)=k_r\prod_i c_i^{\alpha_{ir}}.

Keep $k_r$ as Lecture 3's concentration-scale constant. The displayed $v_r$ is reaction-event flux in molarity per time, before stoichiometry converts it to species changes. If the total reactant order is $q_r=\sum_i\alpha_{ir}$ , then $k_r$ has units $\mathrm M^{1-q_r}/\mathrm{time}$ . We will derive a count-scale law without changing what this constant means.

A composite channel uses the same count bookkeeping and a bare squiggle, with its effective law specified separately. A species symbol $X_i$ names the molecule, not its count. Write $N_{X_i}(t)$ for its random count and collect those counts into the random vector

N_X(t)=(N_{X_1}(t),\ldots,N_{X_s}(t))^T\in\mathrm{Z}_{\geq0}^{s}.

(8)

A particular integer state is $n=(n_1,\ldots,n_s)^T$ . Its component $n_i$ is a realized count of species $X_i$ . We use $s$ for the number of species so that $n$ remains available for the state. Concentrations are written $c_i$ . For a one-species example, $N(t)$ abbreviates that species' random count.

Fix the well-mixed volume $V$ in liters. Let $N_A$ be Avogadro's constant and define the count-to-concentration factor $\Omega=N_AV$ , with units inverse molar. At count state $n$ , the concentration is $c_i=n_i/\Omega$ . One event changes concentration by $\gamma_r/\Omega$ .

The stoichiometry has not changed. Stack the jump columns into $\Gamma=(\gamma_1,\ldots,\gamma_R)$ , with $R$ reaction channels. Lecture 3 called its concentration vector $x$ . Here we write $c$ , keeping species, counts and concentrations distinct. The same accounting gives

\underbrace{\dot c=\Gamma v(c)}_{\text{Lecture 3: concentration dynamics}},\qquad \underbrace{n\mapsto n+\gamma_r}_{\text{one event of channel }r}.

The new object is the propensity $a_r(n)$ . It is defined by the short-time probability

\Pr\{\text{one }r\text{ event in }[t,t+dt)\mid N_X(t)=n\}=a_r(n)dt+o(dt).

(9)

A propensity is the current rate of reaction events. Its units are events per time, not a probability. For an elementary channel it combines the number of eligible reactant combinations with their reaction hazard. For a composite channel it is an effective event-rate law, just as Lecture 3 supplies an effective flux separately from a bare squiggly arrow. It need not have a mass-action form.

Constant addition $\varnothing\rightsquigarrow X$ has jump $+1$ and clock $b$ . First-order removal $X\rightsquigarrow\varnothing$ has jump $-1$ and clock $dn_X$ . Here $n_X$ is the realized count of species $X$ . These are composite channels, so neither effective clock is written on its squiggle.

Count eligible pairs before assigning their hazard

For elementary $A+B\xrightarrow{k}C$ , each of the $n_A$ molecules can partner with any of the $n_B$ molecules. There are $n_A n_B$ eligible pairs. For elementary $2A\xrightarrow{k}C$ , choose one molecule in $n_A$ ways and a different one in $n_A-1$ ways. This counts each physical pair twice, once in each order. The unordered-pair count is therefore

\binom{n_A}{2}=\frac{n_A(n_A-1)}{2!}.

Let $\kappa$ be the firing hazard of one unordered pair at the stated volume. Then $a(n)=\kappa\binom{n_A}{2}$ . This is a count of available combinations times their individual hazard. It is zero when fewer than two molecules exist.

Recover the same concentration constant in the macroscopic limit

For different reactants, consistency with $v=kc_Ac_B$ gives $\kappa=k/\Omega$ . Indeed, $a(n)/\Omega=\kappa n_An_B/\Omega=kc_Ac_B$ . For identical reactants the pair count has an extra denominator. Using $c_A=n_A/\Omega$ ,

\frac{a(n)}{\Omega} =\frac{\kappa}{2\Omega}n_A(n_A-1) =\frac{\kappa\Omega}{2}c_A\left(c_A-\frac1\Omega\right).

At large copy number this must approach the declared event flux $kc_A^2$ . Therefore $\kappa\Omega/2=k$ , not $\kappa\Omega=k$ . The conversion and the resulting propensity are

\kappa=\frac{2k}{\Omega},\qquad a(n)=\frac{2k}{\Omega}\frac{n_A(n_A-1)}{2!} =\frac{k}{\Omega}n_A(n_A-1).

The factorial has canceled against the conversion to the per-pair hazard. We have not redefined $k$ . A count-only text may call $k_{\rm count}=\kappa/2!=k/\Omega$ simply “ $k$ .” Here that would conceal a change of units and volume dependence, so we retain the distinction.²

One convention across reaction orders

A propensity has units events/time, dimensionally inverse time. A concentration event flux has units molar/time. Rate constants acquire these units only after multiplication by the appropriate concentrations or counts. Every $k$ below retains its concentration-law meaning, with $\Omega=N_AV$ .

Channel / order	Concentration event flux	Count propensity	Constant units
constant source	$v=j$	$a=\Omega j=b$	$j:\mathrm{M}/\mathrm{time}$ , $b:1/\mathrm{time}$ (one molecule/event)
first-order removal	$v=dc_A$	$a=dn_A$	$d:1/\mathrm{time}$ in both
elementary $A+B$	$v=kc_Ac_B$	$a=kn_An_B/\Omega$	$k:1/(\mathrm{M}\,\mathrm{time})$ , $k/\Omega:1/\mathrm{time}$
elementary $2A$	$v=kc_A^2$	$a=kn_A(n_A-1)/\Omega$	same units as the preceding row, with per-unordered-pair hazard $2k/\Omega$

For the general elementary mass-action model, define the falling factorial $(n)_m=n(n-1)\cdots(n-m+1)$ , with $(n)_0=1$ and value zero when $m>n$ . The same conversion reads

a_r(n)=k_r\Omega^{1-q_r}\prod_i(n_i)_{\alpha_{ir}} =\underbrace{\left(k_r\Omega^{1-q_r}\prod_i\alpha_{ir}!\right)}_{\kappa_r:\ \text{hazard per combination}} \prod_i\binom{n_i}{\alpha_{ir}}.

Each factorial belongs to one reactant species. It is not a factorial of the total reaction order. This mass-action construction assumes exchangeable reactant combinations in a well-mixed volume. A composite effective flux alone does not determine a unique stochastic mechanism. The exposition derives the general conversion and its limits.

Stoichiometry supplies a different factor of two. Each $2A\xrightarrow{k}C$ event removes two $A$ molecules. Its deterministic contribution is $\dot c_A=-2kc_A^2$ , and its exact conditional count drift is $-2a(n)$ . At finite count, $a(n)/\Omega=kc_A(c_A-1/\Omega)$ , not $kc_A^2$ . The rate-constant conversion is exact within the stated jump model. Replacing the falling factorial by a power is a large-copy-number approximation.

Figure 5. The stochastic reaction contract has three columns. State says what presently exists. Jump says what one firing changes. Propensity says how quickly that firing occurs now. The master equation and simulator must both be generated from this same table.

Turn any reaction list into a stochastic model by filling three columns: the integer state, one jump vector per channel, and one propensity per channel in events per time. Everything after this section is generated from that table and from nothing else.

6The chemical master equation

Probability reaches a state through predecessor events and leaves it through clocks running there.

Let $p(n,t)=\Pr\{N_X(t)=n\}$ . To occupy state $n$ at time $t+dt$ , one of two mutually exclusive things happens to first order in $dt$ :

The system was already at $n$ and no reaction fired.
It was at predecessor $n-\gamma_r$ and reaction $r$ fired once.

The probability of two or more events is $o(dt)$ . Therefore

p(n,t+dt)=p(n,t)\left[1-\sum_r a_r(n)dt\right]+\sum_r p(n-\gamma_r,t)a_r(n-\gamma_r)dt+o(dt).

(10)

Subtract $p(n,t)$ , divide by $dt$ , and take the limit:

Chemical master equation (CME)

\frac{\partial p(n,t)}{\partial t}=\sum_r\left[ \underbrace{a_r(n-\gamma_r)p(n-\gamma_r,t)}_{\text{incoming probability flow}} -\underbrace{a_r(n)p(n,t)}_{\text{outgoing probability flow}}\right].

(11)

This is the dynamical equation for the entire count distribution. Write one equation for each admissible integer state $n$ . Probabilities are dimensionless; each term has units inverse time. An inaccessible predecessor contributes zero. From a reaction table and an initial distribution, the CME predicts zero-count probabilities, fractions above a threshold, and the distribution at any time. It is the probability-level counterpart of Lecture 3's $\dot c=\Gamma v(c)$ , not an equation for a noisy concentration trajectory.

This is the chemical master equation. The first term for channel $r$ is evaluated at the predecessor because its clock rings before the jump. The second term removes probability from the state being described. The equation is linear in $p$ , even if propensities are nonlinear in $n$ .

For constant addition and first-order removal, write $p_n(t)=\Pr\{N(t)=n\}$ . Reading every arrow that touches count $n$ gives

\dot p_n=b p_{n-1}+d(n+1)p_{n+1}-(b+dn)p_n.

(12)

Figure 6. The count ladder writes equation (12). Each arrow carries the probability at its starting state multiplied by the clock running there. Focus on one node and read incoming minus outgoing flow.

The same equation on a finite state graph

A promoter can be off or on. Let $\alpha$ be its off-to-on switching rate and $\beta$ the reverse rate. For the probability column $p=(p_0,p_1)^T$ , the same flow accounting is

\frac{dp}{dt}=Qp,\qquad Q=\begin{bmatrix}-\alpha&\beta\\\alpha&-\beta\end{bmatrix}.

Column $j$ lists departures from state $j$ : off-diagonal entries are rates into other states, and the diagonal is minus their sum. Each column sums to zero. This is a continuous-time Markov chain. $Q$ contains rates, not transition probabilities. The infinite count ladder is the same construction with infinitely many states. These arrows are state transitions, not chemical reaction steps.

How big is this system of equations?

Equation (12) is one equation for every non-negative integer, coupled to its neighbours. Its state space is infinite. Even allowing only 0 to 100 copies of each of ten species gives $101^{10}$ possible states before constraints. Conservation may make a state space finite and much smaller, but enumeration still grows rapidly.

Numerical state-space truncation needs a bound on the probability and observables discarded. Some special networks have analytic solutions, including Section 8's ladder. For larger networks, direct sampling avoids enumerating every possible count. A rare observable can still require many samples.

Write the master equation for any reaction list you can draw. Pick one state, add the flow arriving on each incoming arrow, subtract the flow leaving on each outgoing one, and evaluate each propensity at the state the arrow starts from. No index gymnastics is involved and the same reading works on a finite graph.

7Zero is a biological boundary

Negative counts are not states. Which clocks remain at zero determines what the biology can do next.

At $n=0$ , the removal propensity $dn$ is zero. There is no state $-1$ , so set $p_{-1}=0$ or omit that predecessor explicitly. Equation (12) becomes

\dot p_0=d p_1-bp_0.

(13)

Now sum equation (12) over every non-negative count. Shift the index in each incoming sum. Every internal addition flow cancels one outgoing addition flow, and every removal flow cancels in the same way. The result is

\frac{d}{dt}\sum_{n=0}^{\infty}p_n(t)=0.

(14)

This is the probability analogue of a conservation law. If the sum does not remain one, an index, boundary, or unlisted outside state is wrong.

The biological meaning can be stronger. For an autocatalyst whose addition clock is proportional to its own count, both addition and removal stop at zero. Zero is then absorbing. A deterministic solution may approach zero without ever reaching it, while a molecular path can hit zero exactly and never recover.

Figure 7. Two reaction mechanisms give two different zeros. Constant addition lets the process leave zero. Autocatalytic addition cannot restart without a molecule already present. The boundary behavior is encoded by the propensities, not patched into the simulation afterward.

Audit any master equation you write by summing it over all states: the total must be exactly zero, and a non-zero answer localises the error to a boundary, an index shift, or a state you forgot to list. Then read the boundary as biology, because whether a clock survives at zero is the difference between a pause and an extinction.

8One solvable system calibrates the model

Constant addition plus first-order removal makes equation, distribution, moments, and simulation meet at one answer.

Take $b>0$ and $d>0$ . Define the rightward probability current across the edge from $n$ to $n+1$ as $J_n=bp_n-d(n+1)p_{n+1}$ . Equation (12) reads $\dot p_n=J_{n-1}-J_n$ . At stationarity the current is the same along every edge. Equation (13) forces $J_0=0$ , so it vanishes everywhere on this ladder. Consequently,

b p_n=d(n+1)p_{n+1}\quad\Longrightarrow\quad p_{n+1}=\frac{b/d}{n+1}p_n.

(15)

Define $\mu=b/d$ . Iteration gives $p_n=p_0\mu^n/n!$ . Normalization finishes the derivation:

1=\sum_{n=0}^{\infty}p_n=p_0\sum_{n=0}^{\infty}\frac{\mu^n}{n!}=p_0e^{\mu}\quad\Longrightarrow\quad p_n=e^{-\mu}\frac{\mu^n}{n!}.

(16)

The stationary count is Poisson with

\langle N\rangle=\operatorname{Var}(N)=\mu=\frac{b}{d},\qquad F=\frac{\operatorname{Var}(N)}{\langle N\rangle}=1.

(17)

Can we get a mean or variance without solving every probability? Choose an observable $g(n)$ , a numerical quantity computed from the current count. For example, $g(n)=n$ measures count, $g(n)=n^2$ supplies the second moment, and $g(n)=\mathbf1_{\{n=0\}}$ tests whether the cell is empty. We need the expected rate of change of that observable.

Condition on $N(t)=n$ and look ahead by a small time $h$ . Addition changes $g$ by $g(n+1)-g(n)$ with probability $bh+o(h)$ . Removal changes it by $g(n-1)-g(n)$ with probability $dnh+o(h)$ . No event contributes zero change. Thus

\begin{aligned} \mathrm E[g(N(t+h))-g(n)\mid N(t)=n] ={}&bh[g(n+1)-g(n)]\\ &+dnh[g(n-1)-g(n)]+o(h). \end{aligned}

The infinitesimal generator acts on observables

The generator $\mathcal L$ maps a function $g$ to its conditional expected instantaneous change:

(\mathcal Lg)(n)=\lim_{h\downarrow0}\frac{\mathrm E[g(N(t+h))-g(n)\mid N(t)=n]}{h}.

The parentheses mean “apply the operator to $g$ , then evaluate at state $n$ .” They do not denote a new species or multiplication by the count. Its units are observable units per time. Dividing the preceding event calculation by $h$ gives

(\mathcal{L}g)(n)=b[g(n+1)-g(n)]+dn[g(n-1)-g(n)].

(18)

At zero, omit the removal term because its propensity is zero. For a general network, the same rule is $(\mathcal Lg)(n)=\sum_r a_r(n)[g(n+\gamma_r)-g(n)]$ . It is a rate-weighted finite difference, the jump-process counterpart of a directional derivative.

Average over the current state to obtain observable dynamics. Provided the required expectations are finite, $d\langle g(N)\rangle/dt=\langle(\mathcal Lg)(N)\rangle$ . This is the same probability bookkeeping as the CME, with the sums collected by observable instead of destination state. In Section 6's finite-state convention $\dot p=Qp$ , the observable column is acted on by $Q^T$ , since $d(g^Tp)/dt=(Q^Tg)^Tp$ .

Choose the observable	Change at addition / removal	Generator
$g(n)=n$	$1,\ -1$	$b-dn$
$g(n)=n^2$	$2n+1,\ -2n+1$	$2bn+b-2dn^2+dn$

Write $m=\langle N\rangle$ and $s_2=\langle N^2\rangle$ . The second row gives $\dot s_2=2bm+b-2ds_2+dm$ . Variance is $\sigma^2=s_2-m^2$ , so $\dot\sigma^2=\dot s_2-2m\dot m$ . Substitution gives

\frac{d\langle N\rangle}{dt}=b-d\langle N\rangle,\qquad \frac{d\operatorname{Var}(N)}{dt}=b+d\langle N\rangle-2d\operatorname{Var}(N).

(19)

Choose $b=4\,\mathrm{min}^{-1}$ and $d=0.2\,\mathrm{min}^{-1}$ . The predicted stationary mean and variance are $b/d=4/0.2=20$ , with relative standard deviation $1/\sqrt{20}\simeq0.22$ . Stability has given us a level. Event statistics have added a width.

Numerical comparison

Five thousand independent paths, each initialized at zero and sampled at 60 minutes, give sample mean $19.975$ and unbiased variance $19.567$ . At twelve relaxation times the exact transient mean is $20(1-e^{-12})\simeq19.99988$ . The standard error of the sampled mean is approximately $\sqrt{20/5000}\simeq0.063$ . This is a finite-ensemble comparison, not an exact equality of a histogram and a curve.¹⁹

The zero-current argument used the ladder's topology and boundary. A network with cycles can instead have a stationary distribution carrying circulating currents. Stationarity means that each node's incoming and outgoing currents balance, not that every edge separately balances. This distinction will matter for driven chemistry.

Figure 8. The solvable process is a unit test, not a universal noise law. The analytic curve and simulated histogram agree within finite-sample variation. A new implementation should pass this check before a nonlinear biological result is interpreted.

Why this closes

The propensities are affine in count, so the mean and variance equations close. If a clock contains $N^2$ , then the mean equation contains $\langle N^2\rangle$ , whose equation can contain a third moment. Agreement between the deterministic equation and the exact stochastic mean is special here.

Back to Lecture 3: the same Gamma, different questions

Choose each count coordinate as the observable in the general generator. The exact first-moment equation and its concentration version at fixed volume are

\frac{d\,\mathrm E[N_X]}{dt}=\Gamma\,\mathrm E[a(N_X)],\qquad \frac{d\,\mathrm E[N_X/(N_AV)]}{dt}=\frac{\Gamma}{N_AV}\,\mathrm E[a(N_X)].

The right-hand sides have units molecules/time and molar/time, respectively. For affine clocks, expectation passes through the rate law and the mean closes exactly. For a bimolecular channel between distinct species $X_i$ and $X_j$ , however,

\mathrm E[N_{X_i}N_{X_j}]=\mathrm E[N_{X_i}]\mathrm E[N_{X_j}]+\operatorname{Cov}(N_{X_i},N_{X_j}).

A covariance enters the exact mean reaction rate. The symbols inside the expectation are random counts, not species names.

Lecture 3's $\dot c=\Gamma v(c)$ is the macroscopic law obtained when scaled count fluctuations become negligible and the count propensities approach the corresponding concentration fluxes. Replacing $\mathrm E[a(N_X)]$ by $a(\mathrm E[N_X])$ at small copy number is an additional approximation, not an identity of the CME. Section 11 tests the system-size limit.

Calibrate a new implementation against the stationary targets: mean and variance $b/d$ , and Fano factor one. Allow for transients and finite-sample uncertainty. A persistent discrepancy beyond those effects calls for an implementation or sampling audit before a biological interpretation.

Part 3

Infinitely many states, one trajectory at a time

Here is the payoff. A cell model can have unboundedly many count states and no closed-form distribution. We can still sample its histories exactly, one event at a time, using two uniform random numbers per event. The loop fits on one page and there is no timestep to choose.

9Sample the event, not the timestep

Even when the distribution is too large to compute directly, its next event can be sampled without enumerating all states.

Section 6 gave us a potentially infinite system of coupled equations. Section 8 solved a useful special case, but most networks lack such a closed form. We can still sample their histories and estimate observables with controlled sampling uncertainty. The direct method visits only the current state and its enabled reactions. It does not first compute $p(n,t)$ on the whole state space.

The reason it works is one property of the model. Between events the state does not change, so no propensity changes either, so the waiting time to the next event is drawn from a fixed distribution rather than an evolving one. That single fact turns an infinite-dimensional problem into two draws from a uniform.

A small-timestep simulator repeatedly asks whether an event occurred during $\Delta t$ . It becomes slow when most steps are empty and biased when $\Delta t$ is too large. The direct method asks a different question: how long until the next event?

At state $n$ , let $a_0(n)=\sum_r a_r(n)$ . Let $S(\tau)$ be the probability that no event occurs during the next $\tau$ . Conditional on no event, the state and every propensity are unchanged. Therefore

S(\tau+d\tau)=S(\tau)[1-a_0d\tau]\quad\Longrightarrow\quad \frac{dS}{d\tau}=-a_0S,\qquad S(0)=1.

(20)

Solving gives $S(\tau)=e^{-a_0\tau}$ . If $u_1$ is uniform on $(0,1)$ , inverse-transform sampling gives the next waiting time

\tau=-\frac{\ln u_1}{a_0(n)}.

(21)

Which channel fires? Independent exponential clocks race. Conditional on an event, channel $r$ wins with probability

\Pr\{r\mid\text{an event}\}=\frac{a_r(n)}{a_0(n)}.

(22)

Now perform one complete step. At $n=3$ , with $b=4$ and $d=0.2$ per minute, the two propensities are $a_+=4$ and $a_-=0.6$ , so $a_0=4.6$ . Take $u_1=0.25$ . Equation (21) gives $\tau=0.3014\,\mathrm{min}$ . Take $u_2=0.93$ . Its threshold is $u_2a_0=4.278$ . Addition occupies $[0,4)$ and removal occupies $[4,4.6)$ , so removal wins. Advance the time, update $n:3\mapsto2$ , and rebuild both clocks.

Figure 9. Two uniforms produce one direct-method event. The first becomes an exponential waiting time. The second lands on a propensity interval. Recomputing every state-dependent propensity after the jump is part of the algorithm.

Compute every $a_r(n)$ and their sum $a_0$ .
If $a_0=0$ , the state is absorbing. Stop.
Draw $u_1,u_2$ independently and uniformly from $(0,1)$ .
Set $\tau=-\ln(u_1)/a_0$ .
Choose the first channel whose cumulative propensity exceeds $u_2a_0$ .
Update $t\leftarrow t+\tau$ and $n\leftarrow n+\gamma_r$ .
Recompute propensities in the new state and repeat.

That is the whole algorithm, and it is the one Gillespie published in 1977 in exactly this form.³ Notice what is absent. No timestep, so nothing to converge. No truncation of the state space, so no upper count to justify. No solution of anything. The cost per event does not depend on how many states exist, only on how many reaction channels there are, which is why a network with more states than atoms in the universe still runs.

Live event history

Every path starts at zero. The solid smooth curve is the exact transient mean $m(t)=(b/d)(1-e^{-dt})$ ; the dashed horizontal line is its stationary limit $b/d$ . Lower the removal rate to see why they must not be confused. The jagged path is sampled by the direct method.

addition

b

4.0 /min removal

d

0.20 /min

Theory and path readout appear here.

Sample a history of the stated Markov reaction network from its jump-and-propensity table. Use the two-uniform event loop without choosing a timestep or enumerating the probability state space. Check the physical assumptions in Section 12 before treating exact sampling of the model as an adequate account of the cell.

10Turn trajectories into estimates

The event loop produces histories. An estimator specifies how those histories answer a mean, probability or time-occupancy question.

Define the observable before averaging. Write $Y(t)=g(N_X(t))$ , where $g$ is a scalar function of the count state. Choosing one species count gives its mean. Choosing an indicator, equal to one when that count exceeds a threshold and zero otherwise, gives an exceedance probability. There are two sampling operations to distinguish: independent paths at one time and intervals along one path.

Independent paths estimate a fixed-time ensemble

Sample every path at the same physical time. Generate $M$ independent Gillespie histories with the same parameters and the same initial distribution. For $M\geq2$ , let $Y_\ell(t)=g(N_X^{(\ell)}(t))$ be the value from path $\ell$ . Then

\widehat m_g(t)=\frac1M\sum_{\ell=1}^{M}Y_\ell(t),\qquad \widehat{\mathrm{SE}}[\widehat m_g(t)]=\frac{s_g(t)}{\sqrt M},\qquad s_g^2(t)=\frac1{M-1}\sum_{\ell=1}^{M}[Y_\ell(t)-\widehat m_g(t)]^2.

The target is $\mathrm E[g(N_X(t))]$ . For an indicator, the same average estimates a probability. For a mean curve, repeat the calculation at each desired time. No stationarity or ergodicity assumption is required. Independent paths and finite variance give the $M^{-1/2}$ standard-error scaling at a fixed time. Errors at different times on the resulting mean curve are correlated.

Extending the paths is different from adding paths. Simulating farther into the future does not supply more independent observations at the original time $t$ . It changes the time that can be studied. Increasing $M$ reduces sampling error for the fixed-time target. At later times the target variance can itself change, so its error need not decrease.

A single path gives a residence-time average

A Gillespie state persists for its entire holding interval. Suppose the state is $n_j$ on $[t_j,t_{j+1})$ , and the recorded intervals cover $[0,T]$ . Integrating the piecewise-constant path gives

\overline g(T)=\frac1T\int_0^T g(N_X(t))\,dt =\frac1T\sum_j g(n_j)(t_{j+1}-t_j).

The weights are elapsed times. In a worked four-minute record, the count is 2 for 1 minute, 3 for 0.2 minutes and 2 for 2.8 minutes. Its time mean and fraction of time at count at least 3 are

\overline N_4=\frac{2(1)+3(0.2)+2(2.8)}4=\pf{2.05},\qquad \overline{\mathbf1_{N\geq3}}=\frac{0.2}{4}=\pf{0.05}.

Giving the three recorded states equal weight would instead give mean $7/3$ . That samples event records. It overrepresents states whose clocks ring quickly. Predetermined equally spaced observation times approximate the time integral as the grid is refined, but nearby observations remain correlated. Include the final interval up to the observation horizon, even if the process has become absorbing.

Figure 10. Match the estimator to the sampling operation. Left: the worked event record above. Interval widths supply the time weights. Right: independent-path standard error divided by the observable's fixed-time standard deviation,

1/\sqrt M

. The curve assumes finite variance. It does not describe the precision gained by adding frames to one path.

The live simulation uses the same weights over minutes 20–30. Only the part of each holding interval inside that window contributes. Its displayed path mean and its exact expected window mean answer a finite-window question. Neither is automatically the stationary value $b/d$ .

Ergodicity licenses a stationary interpretation

The finite-window time average is defined without ergodicity. An additional question is whether a longer history estimates a stationary ensemble mean. Let $\pi(n)$ be the stationary law. In an appropriate ergodic regime, for an integrable observable,

\overline g(T)\longrightarrow\sum_n g(n)\pi(n)\qquad\text{as }T\to\infty.

Here ergodicity means that one typical long history samples the stationary law relevant to its accessible states. Constant addition with $b>0$ and removal with $d>0$ provides our working example. A nonstationary initial condition adds finite-window bias. Discarding an initial interval can reduce that bias, but does not prove that a slowly switching process has explored all relevant states.

Correlation sets the precision of a stationary time average. Let $\sigma_g^2$ be stationary variance and $\rho_g(\tau)$ the normalized autocorrelation. Define $\tau_{\mathrm{int}}=\int_0^\infty\rho_g(\tau)\,d\tau$ . When correlation decays sufficiently fast and $T$ is long compared with that memory,

\mathrm{SE}(\overline g(T))\approx\sqrt{\frac{2\sigma_g^2\tau_{\mathrm{int}}}{T}},\qquad M_{\mathrm{eff}}\approx\frac{T}{2\tau_{\mathrm{int}}}.

For the stationary addition-removal count, $\rho_N(\tau)=e^{-d\tau}$ , so $\tau_{\mathrm{int}}=1/d$ . Faster turnover at fixed $b/d$ preserves snapshot variance while increasing the information in a fixed-duration movie. The exposition derives the finite-window variance and the limits of the error approximation.

An autocatalyst for which the long-path substitution fails

Return to the absorbing boundary of Section 7. Take composite autocatalytic addition $X\rightsquigarrow2X$ and removal $X\rightsquigarrow\varnothing$ , with propensities $a_+(n)=\beta n$ and $a_-(n)=\beta n$ . Here $\beta>0$ has units inverse time and $N(0)=1$ . Equal clocks give $d\mathrm E[N]/dt=0$ , hence $\mathrm E[N(t)]=1$ at every finite time. Yet every path eventually reaches zero and remains there, with probability one.

Rare survivors carry the mean. In this critical case, $\Pr[N(t)>0]=1/(1+\beta t)$ and $\mathrm E[N(t)\mid N(t)>0]=1+\beta t$ . At $\beta t=99$ , only 1% survive and their mean count is 100. The full ensemble mean is still 1. Along almost every single path, however, $T^{-1}\int_0^T N(t)dt\to0$ . The positive population has no stationary ergodic regime that could justify replacing this ensemble mean by one long-path average.

The exposition derives extinction and explains the noncommuting average and long-time limit. Absorption alone is not a universal test for disagreement. When removal exceeds autocatalytic addition, both the long-time ensemble mean and the path mean approach zero.

Choose independent paths for a fixed-time ensemble and residence-time weights for a path integral. State the target before invoking ergodicity. Report uncertainty using the number of independent paths or the correlation time of a stationary record, rather than the number of stored events.

11One model has three honest views

A path, a probability distribution, and a deterministic trajectory are related outputs, not interchangeable evidence.

The same Gillespie implementation supplies the samples needed for all three views. Retain one event history to study a path. Repeat the run independently and take snapshots at the same time to estimate a distribution or mean curve. Increase volume at fixed concentration and compare scaled paths to examine a deterministic limit. The program is the same. The sampling and scaling operations differ.

View	What to obtain from the simulation	What it answers
One path	Keep event times and states.	Individual histories, residence times and first passages.
Probability law and ensemble mean	Repeat independent paths. Apply Section 10's fixed-time estimators.	Fractions of realizations, uncertainty and expected trajectories.
Deterministic concentration limit	Repeat at increasing volumes, scaling propensities and counts consistently.	The smooth law approached as relative count fluctuations vanish.

The calibration model has an additional exact relation. Define $m(t)=\mathrm E[N(t)]$ . Its affine propensities make equation (19) close exactly, so this finite-system ensemble mean follows the deterministic equation. A perturbation $\delta m$ of that mean relaxes according to

m(t)=m(0)e^{-dt}+\frac bd(1-e^{-dt}),\qquad \frac{d}{dt}\delta m=-d\,\delta m.

The slope $-d$ is Lecture 4's stability test: perturbations of the mean decay. Stability does not stop reactions. Even at stationarity, addition and removal each have mean event rate $b$ . Their difference vanishes while reactions continue.

Every comparison needs three labels:

Initial condition. A stationary sample and a population initialized at zero answer different questions.
Observation time. At $t\ll1/d$ , the ensemble has not reached its stationary Poisson law.
Averaging operation. A fixed-time path ensemble, a residence-time average and a system-size limit use different operations. Section 10 specifies when a long time average can estimate a stationary expectation.

The deterministic limit is also a declared comparison. Scale a reference volume $V_{\mathrm{ref}}$ by the dimensionless factor $R_V=V/V_{\mathrm{ref}}$ while holding concentration and macroscopic rate constants fixed. This is distinct from the dimensional conversion factor $\Omega=N_AV$ in Section 5. In our model use $a_+(n)=R_Vb$ and $a_-(n)=dn$ , where $b$ is the reference-volume addition propensity. Define $z=N/R_V$ , the count per reference volume. A jump changes it by $1/R_V$ , and at stationarity

\langle z\rangle=\frac bd,\qquad \operatorname{Var}(z)=\frac{b}{dR_V},\qquad \operatorname{CV}(z)=\sqrt{\frac{d}{R_Vb}}.

This is how a smooth concentration law emerges from faster, smaller relative jumps. More generally, density-dependent jump processes converge to deterministic chemical kinetics under the scaling made precise by Kurtz.⁵

Figure 11. The smooth trajectory is earned by system-size scaling. Each colored line is a computed normalized count path

z=N/R_V

. Actual molarity is

z/(N_AV_{\mathrm{ref}})

. The dashed curve is the deterministic solution. Larger systems fluctuate more in absolute molecule number and less in relative concentration.

The boundary can survive the limit question

A large mean may justify an ODE near its central trajectory while leaving a rare extinction or threshold-crossing probability biologically decisive. Model adequacy belongs to the observable, not to copy number alone.

Attach three labels to every comparison you make between a model and data: the initial condition, the observation time, and the averaging operation. Two of the three are usually left implicit, and a disagreement that traces to an implicit label is not evidence about the mechanism.

Part 4

Make a biological claim

An exact algorithm cannot rescue an incomplete state. Finish by naming the world the clocks describe and the observation that could prove it inadequate.

12Exact for a stated world

The direct method exactly samples the stated well-mixed Markov jump process.

The exactness claim has four hypotheses. Molecular amounts are represented as discrete counts. The reaction volume is well mixed. Every hazard is determined by the current modeled state. Propensities stay fixed between reaction events. Under those assumptions, the direct method has no timestep approximation and samples the same path law defined by the CME.²⁴

A finite simulation horizon also requires nonexplosion: only finitely many events occur in a finite time interval. Our calibration and autocatalytic examples satisfy this condition. A general reaction network must be checked before assuming its event loop can reach every requested time.

It is not exact molecular physics. Spatial gradients can make location part of the state. Crowding and rebinding can create memory. Transcriptional elongation can create a non-exponential delay. Hidden conformations can change a waiting-time law. Cell growth changes volume and therefore count-scale propensities. Each failure points to a repair: enlarge the state, add space or age, or change the model class.

Figure 12. Exactness is nested inside a model boundary. The direct method samples the inner process exactly. Processes in the outer biological world matter only if the observable is sensitive to them, but then they must be represented or declared absent.

State the four hypotheses whenever you call a simulation exact: discrete counts, a well-mixed volume, hazards fixed by the present modelled state, and propensities constant between events. When one of them fails, name the repair it points to rather than the algorithm, because the algorithm was not the thing that broke.

13Make the observable choose the model

A stochastic model is complete only when it names what it predicts and what would reject it.

For a new biological network, produce this six-line contract:

State. Name each integer component and the legal state space.
Events. List every reaction, jump vector, and propensity with units.
Probability. Write one interior CME equation and every special boundary equation.
Path. Hand-trace one direct-method event with both uniform draws recorded.
Observable and estimator. Choose a mean, Fano factor, zero probability, waiting time or first-passage statistic. State whether it comes from independent paths or a time-weighted record, with its observation window and uncertainty.
Falsifier. Name one data pattern the model cannot produce and one missing state or process that could repair it.

For example, a stationary Fano factor reliably above one, after accounting for sampling uncertainty, rejects the constant-addition, first-order-removal Poisson baseline. It does not identify the replacement. Promoter switching, translation packets, feedback, shared cell-cycle state, and mixed subpopulations can all widen a count distribution. Models that agree on a histogram can disagree on inter-event waiting times.¹¹ Keep event times when the claim concerns a mechanism's internal steps.

Figure 13. The model ends at a falsifier. The observable comes first, because it decides which state and representation are adequate. A visually plausible jagged trace is not a validation criterion.

The core story

An engineered oscillator makes single-cell histories a biological question. The state and observable decide which sources of variation matter. The same reaction accounting as Lecture 3 supplies integer jumps and propensities. The CME evolves their probabilities, the generator evolves observables, and Gillespie samples their histories without enumerating the state space. The estimator specifies how those histories become a probability, mean or time-occupancy claim. The deterministic rate equation emerges under a stated scale limit and is the exact mean only in special cases. A falsifiable observable decides which description the cell's behavior requires.

Lecture 6 begins with a concrete hidden state. A promoter switches on at rate $\alpha$ and off at rate $\beta$ , makes RNA at rate $k$ while on, and RNA is removed at rate $\gamma n$ . Its stationary on-probability is $p_{\mathrm{on}}=\alpha/(\alpha+\beta)$ . Can we replace the promoter by constant addition $b=k p_{\mathrm{on}}$ ? That replacement preserves the stationary mean, but not generally the fluctuations. Compare the switching time $1/(\alpha+\beta)$ with the RNA lifetime $1/\gamma$ . Lecture 6 will make this separation of timescales an explicit reduction argument.

Turn a reaction list into the six-line modelling contract, including an observable and a falsifier. Use a rejected prediction to identify what must change, while keeping open the possibility that several mechanisms fit the same histogram.

Playground

Outside the lecture clock

Seven models put the event loop to work. Tune them, compare histories with ODEs, and try the optional follow-ups. These are not additional timed lecture points.

▶Seven systems, with tunable parameters

Use the same reactions to compare a stochastic history with a deterministic prediction. Change a parameter, keep the seed fixed, and ask what changed.

The stochastic engine uses Section 9's direct method. Combinatorial clocks turn unavailable mass-action channels off. The toggle instead has explicitly declared effective Hill propensities. At scheduled changes of season, the bet-hedging model advances to the boundary and starts the next race with the new hazards. The dashed comparison uses the corresponding continuous rate laws, not the average of the displayed stochastic path.

Reaction playground

Choose a system to reveal its parameters and units. Apply parameters reruns with the same seed; New seed changes the random history. Reset model restores that example's defaults. All examples use illustrative model time units, not a universal minute.

system

Random seed Show deterministic ODE

Pick a system.

1 · A switch that changes its own mind

One dynamic species, four directed channels, and a bistable concentration law. With chemostatted concentrations $a$ and $b$ , $\dot x=k_1ax^{2}-k_2x^{3}+k_3b-k_4x$ has three fixed points at the default parameters, two attracting and one between them. Deterministically the system stays in its initial basin while the parameters remain fixed. Changing the controls can remove bistability.

\mathrm{A}+2\mathrm{X}\xrightarrow{k_1}3\mathrm{X},\qquad 3\mathrm{X}\xrightarrow{k_2}\mathrm{A}+2\mathrm{X},\qquad \mathrm{B}\xrightarrow{k_3}\mathrm{X},\qquad \mathrm{X}\xrightarrow{k_4}\mathrm{B}.

A count trajectory can change basins without a change outside the cell. At the default volume, the macroscopic fixed points correspond to about 1 and 38 molecules. The high-to-low first-passage benchmark on the count ladder is 131 time units and rises to $2.3\times10^8$ at 3.75 times the volume with the same concentration law. The two escape directions need not have the same waiting time. A stochastic memory has a lifetime, which depends strongly on copy number. These benchmarks use the default kinetic parameters and stated thresholds, not every setting of the controls.¹² That steepness is measurable. Acar and colleagues sorted the two expression states of a bistable galactose circuit and watched them relax back, and the escape rates fell precipitously as the barrier grew, from a steady state reached in about ten hours near the critical point to cells effectively locked in one state.²⁹

Watch for: a run with no flips at all, and how tempting it is to report that one as the answer.

2 · A chemical vote, and what it does with a tie

Two opinions and an undecided state. Disagreement makes both parties undecided; an undecided molecule meeting an opinion adopts it.

\mathrm{X}+\mathrm{Y}\xrightarrow{r}\mathrm{X}+\mathrm{B},\qquad \mathrm{X}+\mathrm{Y}\xrightarrow{r}\mathrm{Y}+\mathrm{B},\qquad \mathrm{X}+\mathrm{B}\xrightarrow{r}2\mathrm{X},\qquad \mathrm{Y}+\mathrm{B}\xrightarrow{r}2\mathrm{Y},\qquad r=1/n.

This is approximate majority, a population protocol that computes which opinion was in the majority and commits every molecule to it.¹³ Consensus is absorbing: once one opinion reaches zero, no propensity is left. From 110 against 90 the majority wins 184 runs out of 200. From a dead tie the answer is a coin flip, 112 out of 200 in the same test, and the coin is the reaction noise. The same four reactions are the topology of the cell-cycle switch.¹⁴

Watch for: how long the system sits near the tie before committing, and that the commitment is irreversible.

3 · The cycle that does not come back

Rabbits breed, foxes eat rabbits and become foxes, foxes die. Every one of those is a lifetime compressed into an arrow, so all three are composite and their rate laws are declared rather than derived: mass action at rates $k_1$ , $k_2$ and $k_3$ .

\mathrm{R}\rightsquigarrow 2\mathrm{R},\qquad \mathrm{R}+\mathrm{F}\rightsquigarrow 2\mathrm{F},\qquad \mathrm{F}\rightsquigarrow\varnothing.

The deterministic Lotka-Volterra model has a centre surrounded by closed orbits. At the default initial state it stays exactly at that centre; change an initial count to see a deterministic orbit. There is no restoring force toward one preferred amplitude. Count fluctuations can change that amplitude and eventually reach a zero boundary, which is absorbing for the lost species. In the archived default-parameter benchmark, 38 of 40 seeded runs reached an extinction within 150 time units while the deterministic populations persisted.

Watch for: which one dies. Losing the prey and losing the predator are different endings with different biology.

4 · Bursts, from a promoter that will not sit still

A gene that switches on and off, transcribing only while on. Every arrow is composite: what makes a promoter accessible is not resolved here, and neither is transcription. Their declared rates are $\alpha$ off to on, $\beta$ back, $k$ for transcription while on, and $\gamma$ per transcript for removal.

\mathrm{G}_{\text{off}}\leftrightsquigarrow\mathrm{G}_{\text{on}},\qquad \mathrm{G}_{\text{on}}\rightsquigarrow\mathrm{G}_{\text{on}}+\mathrm{M},\qquad \mathrm{M}\rightsquigarrow\varnothing.

This is Section 3's packet argument with a mechanism attached. The stationary Fano factor is $1+k\beta/[(\alpha+\beta)(\alpha+\beta+\gamma)]$ , which is 8.5 at the shipped rates against 1 for Poisson, and a long run measures 8.9. The mean is 5 transcripts. Nothing about the mean says the promoter is switching. Everything about the width does. Lecture 6 is this system, worked.

Watch for: the flat stretches at zero. A cell in one of those is not a cell with a low mean.

5 · Losing on purpose, so the lineage survives

Two phenotypes, one environment that changes, and no sensor anywhere. $\mathrm{F}$ grows fast and dies under stress. $\mathrm{S}$ grows slowly and survives it. Cells convert between them at a small constant rate, and crowding is a second-order death so the population stays bounded on its own.

\mathrm{F}\rightsquigarrow 2\mathrm{F},\qquad \mathrm{S}\rightsquigarrow 2\mathrm{S},\qquad \mathrm{F}\leftrightsquigarrow\mathrm{S},\qquad \mathrm{F}+\mathrm{F}\rightsquigarrow\mathrm{F},\ \ldots\qquad \mathrm{F}\rightsquigarrow\varnothing\ \text{under stress}.

Switching creates a reservoir before stress arrives. With no switching, a pure $\mathrm{F}$ culture cannot produce protected $\mathrm{S}$ cells. With switching at $0.02$ per model time unit, a slow-growing minority can survive the stress channel that removes fast cells. This minority costs growth in the good seasons. It is insurance, not a guarantee against extinction.¹⁶¹⁷

Compare the same six-season schedule with and without switching. In fifteen seeded benchmark runs per setting, the median final population is about 150 with switching and zero without it. Changing the season duration or the slow phenotype's cost changes the tradeoff.

Watch for: how the switching rate that improves survival changes when you alter season duration, stress severity, or the slow phenotype's cost. Compare many seeds, not only the most fortunate lineage.

6 · One enzyme, counted

Michaelis and Menten with a single enzyme molecule, so the turnovers are events rather than a rate.

\mathrm{E}+\mathrm{S}\xrightleftharpoons[k_{-1}]{k_{1}}\mathrm{C}\xrightarrow{k_{\mathrm{cat}}}\mathrm{E}+\mathrm{P}.

A turnover requires binding and successful catalysis, with possible unbinding and rebinding in between. Its waiting time is therefore not generally one exponential, nor just the sum of two if retries are allowed. In the special case of irreversible binding at fixed substrate, the two waits give mean $1/(k_1S)+1/k_{\mathrm{cat}}$ . Only in the saturating-substrate limit does its reciprocal approach $k_{\mathrm{cat}}$ . With depletion, later intervals have different statistics. The live readout reports the pooled inter-product intervals and states that limitation.

Watch for: raise the enzyme count to ten and watch the staircase turn into a line.

7 · The genetic toggle: an attractor is not an eternal memory

Return to Lecture 4's mutual repression. Repressor X suppresses synthesis of Y, and Y suppresses synthesis of X. Gene expression and removal are composite channels:

\varnothing\rightsquigarrow X,\quad X\rightsquigarrow\varnothing,\qquad \varnothing\rightsquigarrow Y,\quad Y\rightsquigarrow\varnothing.

Let $q$ be the count per concentration unit and write $x=n_X/q,\ y=n_Y/q$ . This scale is proportional to volume, not a conserved total. Take the repression threshold as one concentration unit, maximum synthesis $a$ , synthesis asymmetry $\rho$ , Hill exponent $h$ , and removal constant $d$ . The paired descriptions are

\begin{aligned} a_{X+}(n)&=\frac{qa}{1+(n_Y/q)^h},&a_{X-}(n)&=dn_X,&\dot x&=\frac{a}{1+y^h}-dx,\\ a_{Y+}(n)&=\frac{q\rho a}{1+(n_X/q)^h},&a_{Y-}(n)&=dn_Y,&\dot y&=\frac{\rho a}{1+x^h}-dy. \end{aligned}

The Hill clocks are lumped effective laws, under fast binding/promoter assumptions, not elementary mass action. At $a=8,\rho=d=1,h=2$ , the deterministic attractors are approximately $(x,y)=(7.87,0.127)$ and its mirror. The intermediate symmetric fixed point is a saddle. These illustrative values are calculated from the displayed equations, not fitted to Gardner and colleagues' experiment.²⁰

Try three contrasts: keep the parameters and change only the seed; increase the count scale from 2 to 20 while keeping initial concentrations fixed; reduce the Hill exponent toward 1 or change the synthesis ratio. A stochastic crossing of an unchanged separatrix and the deterministic disappearance of an attractor are different events. The simulation makes no promise that a finite observation window will contain a switch. The growth addendum below asks what happens when the landscape itself moves.

Where these came from, and where to take them

Approximate majority and the fox-and-rabbit system are from Erik Winfree's molecular computation course, which shipped them as stochastic simulator examples.¹⁸ The bistable switch is Schlögl's, in the form Vellela and Qian analysed.¹² A landmark application is the stochastic model of the lambda phage decision between lysis and lysogeny. Its biochemical network linked variable fate outcomes to the timing of early molecular events.¹⁵ That is the shape of the good problems here. Find a decision a cell makes once, write its reaction list, and ask what decides it.

The six original examples have reproducible benchmarks in playground_models.py. The new toggle's roots and geometry are checked in verify_feedback.py, and browser_audit.mjs tests the controls and displayed calculations. These files live in research/lecture05/analysis/.¹⁹

✦Further experiments and noise tools

Optional continuations, outside the lecture clock. Each question extends a calculation you can already begin from the reaction table.

An oscillator's noise can be diagnosed and engineered

Better observation changes the question we can answer. Potvin-Trottier, Lord, Vinnicombe and Paulsson followed the repressilator in a mother machine, keeping single-cell lineages under sustained growth. Oscillations persisted in all tracked cells under their conditions. This showed that some apparent failures in earlier short observations need not be failures of the oscillator. It did not make the original circuit's single-cell timing noise disappear.²¹

The circuit was then improved by separating several mechanisms. Reporter redesign reduced interference with the tagged-protein degradation machinery. Removing active repressor degradation lengthened the clock, but alone improved timing noise only modestly. Additional TetR-binding sites acted as a titration sponge: more TetR had to be present to keep the next gene repressed, changing the timing threshold. The final design kept a period of roughly fourteen generations with much more regular phase, allowing population synchrony without communication between cells.²¹

A simple calculation explains why a low copy-number threshold is noisy. Suppose a repressor count falls from $N_{\mathrm{peak}}$ to threshold $s$ by independent removal at rate $dn$ . Successive waiting times are independent exponentials. Their sum $T_s$ has

\mathrm E[T_s]=\sum_{n=s+1}^{N_{\mathrm{peak}}}\frac1{dn},\qquad \operatorname{Var}(T_s)=\sum_{n=s+1}^{N_{\mathrm{peak}}}\frac1{d^2n^2}.

The last few molecules dominate the variance. Raising the threshold avoids waiting for those especially uncertain final events. This is an illustrative removal model for threshold timing, not a full reproduction of growth dilution, reporter dynamics, or the redesigned oscillator. Explore: change the threshold and compare relative timing error at the same peak count. Is making the clock slower the same as making it more precise?

A burst is a mechanism, a size distribution, and a timescale

Translation gives a random, not fixed, packet size. During one mRNA's life, translation at rate $k_{\mathrm{tl}}$ races against removal at rate $d_m$ . At each race, translation wins with probability $q_b=k_{\mathrm{tl}}/(k_{\mathrm{tl}}+d_m)$ , after which the same competition starts again. The number $B$ of proteins made before removal is therefore

\Pr(B=j)=(1-q_b)q_b^j,\quad j=0,1,\ldots,\qquad \mathrm E[B]=\bar B=\frac{k_{\mathrm{tl}}}{d_m},\quad \operatorname{Var}(B)=\bar B(1+\bar B).

To replace this episode by one instantaneous jump, the mRNA lifetime must be short relative to the protein timescale. Resolving individual bursts in a movie additionally requires that episodes not overlap too heavily. A large Fano factor alone does not establish visibly isolated packets.²⁵²³

Observation and removal change the packet formula. For independent packets arriving as a Poisson process at rate $r_b$ , a window of length $T$ with no removal gives mean count $r_bT\mathrm E[B]$ and variance $r_bT\mathrm E[B^2]$ . With first-order removal at rate $dn$ , apply the generator to count and count squared instead. At stationarity,

\bar N=\frac{r_b\mathrm E[B]}{d},\qquad F=\frac12\left(1+\frac{\mathrm E[B^2]}{\mathrm E[B]}\right).

For geometric packets, $\mathrm E[B^2]=\bar B+2\bar B^2$ , hence $F=1+\bar B$ . This differs from both the fixed-packet window law in Figure 4 and the corresponding geometric-packet window law $F=1+2\bar B$ . The event mechanism and filtering have to travel with the number.

Which distribution should gene expression have?

The answer is a model-and-observable statement, not “RNA is gamma, protein is log-normal.” RNA and protein molecule counts are integers. Gamma and log-normal laws are continuous densities, useful as approximations or for continuous readouts, not exact count laws.

Assumptions	Observable / prediction	What to check next
constant addition, independent first-order removal	Poisson stationary count, for RNA or protein	Fano one, correct zero fraction, exponential temporal correlation
two-state promoter, RNA made while on	telegraph stationary RNA law; expressible as a Poisson mixture with a beta-distributed intensity	switching times relative to RNA removal; it need not be bimodal
Poisson arrivals of geometric packets, independent removal	negative-binomial count law in the burst limit	packet statistics and separation of lifetimes
appropriate large-packet/continuous-abundance limit	gamma density; also an empirical protein-distribution fit	low-count discreteness and whether fitted parameters identify kinetic rates
approximately Gaussian sum of logarithmic factors	log-normal positive abundance or readout	which biological and measurement factors multiply; zero counts need separate treatment

Golding and colleagues directly observed transcriptional episodes in bacteria, with about two transcripts per episode in their reporter and conditions. Their MS2-labelled RNAs were unusually stable, so division mattered strongly to their removal and partitioning.²⁴ Taniguchi and colleagues found gamma fits useful for most of their measured E. coli protein distributions, with log-normal fits less successful at low expression.⁶ Neither result is an organism-wide distribution law.

Promoter states are hypotheses about memory. A bacterial on/off model can represent slow access or activity changes. A mammalian promoter may need chromatin opening, an accessible inactive state, and an active state. A three-state cycle with two sequential inactive steps has a refractory waiting interval, unlike a single memoryless off state. But two states are often useful in mammalian systems too, and bacteria may require more. Counting states is a model-selection problem. Mammalian transcript counting supports transcriptional episodes, but does not by itself identify the molecular mechanism of the hidden state.¹⁰ Waiting-time data can distinguish mechanisms that fit nearly the same count histogram.¹¹ Explore: compare two models with the same mean and variance but different off-time distributions.

The linear-noise approximation reuses the Jacobian

Lecture 4's restoring matrix also shapes fluctuations. Write a deterministic count-scale drift $F(n)=\Gamma a(n)$ and let $\bar n$ be a stable fixed point in a regime of small relative fluctuations. Define its Jacobian $A=\partial F/\partial n|_{\bar n}$ . Independent reaction clocks add count increments with local covariance per time

D=\sum_r a_r(\bar n)\gamma_r\gamma_r^T =\Gamma\operatorname{diag}(a(\bar n))\Gamma^T.

Linearize the restoring drift and evaluate event noise at the fixed point. The linear-noise approximation (LNA) gives a covariance matrix $C$ satisfying

\dot C=AC+CA^T+D,\qquad AC_*+C_*A^T+D=0.

This is the algebraic Lyapunov equation at stationarity. Here $A$ has units inverse time, $C$ has count-squared units, and $D$ has count-squared per time. Use a stable independent-coordinate subsystem if exact conservation creates zero modes. For addition/removal, $A=-d$ and $D=b+d\bar n=2b$ , so $C_*=b/d$ . That covariance is exact for this affine model; the Gaussian distribution implied by LNA is still not an exact Poisson count law.

Explore: use the toggle's Jacobian near one attractor to estimate local covariance, then compare with SSA restricted to that basin. A local Gaussian cannot describe a switch to the other basin, extinction, or a strongly skewed low-count distribution. Reaction-level correlations also matter: a single event that removes X and adds Y contributes a negative off-diagonal entry to $D$ . Independently adding noise to each ODE would miss it.²

Suppress noise by changing the right mechanism

Two-stage expression has protein Fano factor $F_P=1+k_{\mathrm{tl}}/(d_m+d_p)$ , where $d_p$ is protein removal. At a fixed mean protein level, more transcripts translated fewer times each can reduce intrinsic protein fluctuations. The mean alone would not reveal that design choice.²⁵

A second route is negative feedback. For one species with addition $f(n)$ , removal $dn$ , and stable mean-field fixed point $f(\bar n)=d\bar n$ , LNA gives

C_*=\frac{d\bar n}{d-f'(\bar n)},\qquad F_{\mathrm{LNA}}=\frac{d}{d-f'(\bar n)}.

A negative addition slope strengthens restoration and can produce sub-Poisson variance in this model. It is not a free universal noise bound: an explicit noisy regulator, delay, shared resources, or coupled reaction increments changes the calculation. Explore: hold the mean fixed while varying feedback, calculate both $A$ and $D$ , then test the predicted variance and response time. This is the central calculation behind the problem-set noise-suppression questions.

✦Growth, division, and the lifetime of circuit memory

A cell is both the vessel and part of the mechanism. Growth can move an attractor, push a state toward a boundary, and create correlated or discrete fluctuations.

Two experiments separate robustness from immunity. Zhang and colleagues in Xiao-Jun Tian's group compared self-activation and mutual repression. Transfer into fresh medium stimulated growth and could erase the self-activator's memory; conditioned medium preserved it. The tested toggle recovered its expression state after the fast-growth phase. Their models connected this difference to regulatory topology and relative response times. The result is conditional robustness, not a theorem that every toggle withstands every growth perturbation.²⁶

Zhu, Chu and Fu at SIAT examined a mutually repressive circuit across growth conditions. The two expression arms responded unequally to growth, shifting their nullclines and creating a growth-dependent transition between bistability and monostability. They found no significant growth-rate difference between the two phenotypes in that comparison, so the mechanism was not simply one phenotype outgrowing the other. Their result complements the topology comparison: even a resilient wiring pattern can lose a state when its parameters change asymmetrically.²⁷

First distinguish counts from concentrations, as in Lecture 3. Let $N_i$ count molecules, $V$ be volume in liters, $c_i=N_i/(N_AV)$ be molar concentration, and $\lambda=\dot V/V$ be the instantaneous growth rate. Suppose synthesis flux $f_i(c,\lambda)$ has units molar/time and degradation rate $\delta_i$ has units inverse time. Between divisions, smooth mean-field bookkeeping gives

\dot N_i=N_AVf_i(c,\lambda)-\delta_iN_i,\qquad \dot c_i=\frac{\dot N_i}{N_AV}-c_i\frac{\dot V}{V} =f_i(c,\lambda)-(\delta_i+\lambda)c_i.

The dilution term follows from the changing denominator. If growth changes deterministically with nutrients, this is a deterministic, possibly time-dependent disturbance. If circuit expression changes growth through burden, the equations acquire a feedback loop. If growth differs randomly between cells, the same growth variable can perturb several circuit components together. Those are distinct models.

Growth addendum. Geometry, volume, and partitioning are different operations. Left: a computed self-activator rate balance with green addition and red removal

f^-(x)=x

, where a concentration-reducing pulse can cross the unstable threshold. Centre: nullclines of the symmetric toggle, where a common proportional concentration reduction preserves the side of the diagonal separatrix. These are geometric illustrations, not fits or literal concentration halving at cytokinesis. Right: the exact binomial partition law for a cell with twenty molecules before symmetric division; the daughter concentration varies even when division time and volume ratio are fixed.

Then specify what happens at division. For a symmetric split, daughter volume is $V^+=V^-/2$ . If molecules are independently and fairly allocated to a chosen daughter,

N_i^+\mid N_i^-=n\sim\mathrm{Binomial}(n,1/2),\qquad c_i^+=\frac{2N_i^+}{N_AV^-},\qquad \mathrm E[c_i^+\mid n,V^-]=c_i^-.

Halving count and volume leaves the expected concentration unchanged. It is incorrect to halve concentration just because the cell divided. The random partition gives conditional concentration variance $n/(N_AV^-)^2$ and relative standard deviation $1/\sqrt n$ . At twenty parental molecules that is about 22%; at two hundred, about 7%. Approximately binomial RNA partitioning was measured in Golding and colleagues' reporter system.²⁴

Even exactly scheduled division times can therefore coexist with stochastic state resets. Perfect proportional partitioning would instead give no concentration jump. Unequal daughter volumes, clustered molecules, replication, or shared partition fractions require a richer reset model and may correlate different species' perturbations. A stationary-cell CME with fixed hazards between events does not automatically include any of these effects.

Finally ask whether the disturbance crosses a basin or changes the basin. In a one-dimensional self-activator, a sufficiently strong concentration reduction can cross the unstable threshold. In an ideal symmetric toggle, a common proportional reduction preserves the sign of $x-y$ , so it need not cross the separatrix. Unequal growth responses can instead move or remove that separatrix, as the SIAT work emphasizes. Random partitioning can push cells differently even under the same deterministic growth program.

An extension a student team can actually test

Start with the toggle and self-activator at matched mean concentrations and declared stability margins. Compare four interventions: a common deterministic growth pulse; unequal growth-dependent synthesis; scheduled divisions with proportional partition; and the same schedule with binomial partition. Track state retention, first switching time, and nullclines. Keep volume explicit so growth dilution is not counted twice. A change of attractor is a Lecture 4 bifurcation; a random crossing of a fixed basin is a Lecture 5 path event; the growth term producing either belongs to Lecture 3. Do not claim that the fixed-volume playground already simulates division.

†References

M. B. Elowitz and S. Leibler, “A synthetic oscillatory network of transcriptional regulators,” Nature 403, 335–338 (2000). DOIThe designed three-repressor oscillator and its variable single-lineage behavior.
D. T. Gillespie, “Stochastic Simulation of Chemical Kinetics,” Annual Review of Physical Chemistry 58, 35–55 (2007). DOIPropensities, the CME, direct SSA, and the physical assumptions of stochastic chemical kinetics.
D. T. Gillespie, “Exact stochastic simulation of coupled chemical reactions,” Journal of Physical Chemistry 81, 2340–2361 (1977). DOIThe direct method, with the survival and competition argument in its original form.
D. T. Gillespie, A. Hellander, and L. R. Petzold, “Perspective: Stochastic algorithms for chemical kinetics,” Journal of Chemical Physics 138, 170901 (2013). full textModel-relative exactness and the hierarchy of stochastic simulation methods.
T. G. Kurtz, “The Relationship between Stochastic and Deterministic Models for Chemical Reactions,” Journal of Chemical Physics 57, 2976–2978 (1972). DOIThe density-dependent limit connecting jump processes to deterministic kinetics.
Y. Taniguchi, P. J. Choi, GW. Li, H. Chen, M. Babu, J. Hearn, A. Emili, and X. S. Xie, “Quantifying E. coli proteome and transcriptome with single-molecule sensitivity in single cells,” Science 329, 533–538 (2010). full textAbundances over five decades, the $1/\langle N\rangle$ noise scaling below ten copies, and the extrinsic floor above it.
M. B. Elowitz, A. J. Levine, E. D. Siggia, and P. S. Swain, “Stochastic Gene Expression in a Single Cell,” Science 297, 1183–1186 (2002). DOIThe paired-reporter experiment separating reporter-specific and shared variation.
P. S. Swain, M. B. Elowitz, and E. D. Siggia, “Intrinsic and extrinsic contributions to stochasticity in gene expression,” PNAS 99, 12795–12800 (2002). full textWhere the two definitions in equation (1) come from, and what they do and do not partition.
C. Furusawa, T. Suzuki, A. Kashiwagi, T. Yomo, and K. Kaneko, “Ubiquity of log-normal distributions in intra-cellular reaction dynamics,” Biophysics 1, 25–31 (2005). full textLog-normal abundances with standard deviation proportional to the mean, which is a constant relative spread rather than an event count.
A. Raj, C. S. Peskin, D. Tranchina, D. Y. Vargas, and S. Tyagi, “Stochastic mRNA Synthesis in Mammalian Cells,” PLoS Biology 4, e309 (2006). DOISingle-molecule transcript counts and promoter-state interpretation.
S. Braichenko, J. Holehouse, and R. Grima, “Distinguishing between models of mammalian gene expression: telegraph-like models versus mechanistic models,” Journal of the Royal Society Interface 18, 20210510 (2021). full textCount distributions can agree while inter-event timing distinguishes hidden mechanisms.
M. Vellela and H. Qian, “Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: the Schlögl model revisited,” Journal of the Royal Society Interface 6, 925–940 (2009). full textThe rate constants used in playground system 1, and the divergence between the deterministic and stochastic accounts of the same bistable network. The model is Schlögl’s, from 1972.
D. Angluin, J. Aspnes, and D. Eisenstat, “A simple population protocol for fast robust approximate majority,” Distributed Computing 21, 87–102 (2008). full textThe four-reaction consensus protocol, with its convergence time and error bounds.
L. Cardelli and A. Csikász-Nagy, “The cell cycle switch computes approximate majority,” Scientific Reports 2, 656 (2012). full textThe same four reactions found inside the mitotic entry switch.
A. Arkin, J. Ross, and H. H. McAdams, “Stochastic kinetic analysis of developmental pathway bifurcation in phage $\lambda$ -infected Escherichia coli cells,” Genetics 149, 1633–1648 (1998). full textA mechanistic stochastic model of a cell-fate decision that predates the engineered toggle and repressilator.
E. Kussell, R. Kishony, N. Q. Balaban, and S. Leibler, “Bacterial persistence: a model of survival in changing environments,” Genetics 169, 1807–1814 (2005). full textSpontaneous phenotype switching as insurance, and the optimal switching rate set by how often the environment changes.
M. Thattai and A. van Oudenaarden, “Stochastic gene expression in fluctuating environments,” Genetics 167, 523–530 (2004). full textWhy a population that hedges can outgrow one that optimises for the average environment.
E. Winfree, BE/CS/CNS 191, Molecular Programming, California Institute of Technology, winter 2017. The CRNSimulator stochastic examples, from which playground systems 2 and 3 are taken.
Lecture 5 numerical and figure ledger, seeds 20260902, 20260907, 20260908 and 20260915. The checked generators produce the plotted calculations and benchmarks. The playground exposes its actual seed and parameter settings; its ODE is integrated separately from the SSA.
T. S. Gardner, C. R. Cantor, and J. J. Collins, “Construction of a genetic toggle switch in Escherichia coli,” Nature 403, 339–342 (2000). DOIEngineered mutual repression and cellular memory. The playground uses an illustrative symmetric Hill model, not a fit to these data.
L. Potvin-Trottier, N. D. Lord, G. Vinnicombe, and J. Paulsson, “Synchronous long-term oscillations in a synthetic gene circuit,” Nature 538, 514–517 (2016). DOILong-term single-cell measurements, reporter/degradation interactions, and titration-based improvement.
L. Cai, N. Friedman, and X. S. Xie, “Stochastic protein expression in individual cells at the single molecule level,” Nature 440, 358–362 (2006). DOICatalytic amplification and calibrated protein burst sizes under weak bacterial expression.
V. Shahrezaei and P. S. Swain, “Analytical distributions for stochastic gene expression,” PNAS 105, 17256–17261 (2008). DOIGeometric bursts and negative-binomial protein counts in a declared fast-mRNA limit.
I. Golding, J. Paulsson, S. M. Zawilski, and E. C. Cox, “Real-Time Kinetics of Gene Activity in Individual Bacteria,” Cell 123, 1025–1036 (2005). DOITranscriptional episodes and approximately binomial division partitioning in an MS2-labelled RNA reporter.
M. Thattai and A. van Oudenaarden, “Intrinsic noise in gene regulatory networks,” PNAS 98, 8614–8619 (2001). DOITwo-stage expression moments, translational burst interpretation, and feedback.
R. Zhang, J. Li, J. Melendez-Alvarez, X. Chen, P. Sochor, H. Goetz, Q. Zhang, T. Ding, X. Wang, and X.-J. Tian, “Topology-dependent interference of synthetic gene circuit function by growth feedback,” Nature Chemical Biology 16, 695–701 (2020). DOISelf-activation versus toggle memory under growth-dependent dilution; full author manuscript inspected.
J. Zhu, P. Chu, and X. Fu, “Unbalanced response to growth variations reshapes the cell fate decision landscape,” Nature Chemical Biology 19, 1097–1104 (2023). DOI; Unequal growth dependence shifts the toggle's nullclines. The published text carries the control this section relies on: no significant difference in growth rates between the two phenotypes, so the fate variation is not nonlinear dilution from metabolic burden.
J. R. S. Newman, S. Ghaemmaghami, J. Ihmels, D. K. Breslow, M. Noble, J. L. DeRisi, and J. S. Weissman, “Single-cell proteomic analysis of S. cerevisiae reveals the architecture of biological noise,” Nature 441, 840–846 (2006). DOIThe yeast counterpart of the previous entry, by flow cytometry over 2,500 proteins. Noise tracks stochastic mRNA production and destruction, extrinsic variation dominates at high abundance, and noise sorts by protein function.
M. Acar, A. Becskei, and A. van Oudenaarden, “Enhancement of cellular memory by reducing stochastic transitions,” Nature 435, 228–232 (2005). DOIMeasured escape rates in a bistable galactose circuit, falling steeply with the calculated barrier height. The experimental counterpart of the first-passage benchmark computed above.
B. L. Sabatini, T. G. Oertner and K. Svoboda, “The Life Cycle of Ca2+ Ions in Dendritic Spines,” Neuron 33, 439–452 (2002). DOIFigures 2–4 separate calcium clearance, indicator-loaded equilibration and the inferred native exchange time.