From reaction list to stochastic evidence

A working tutorial from chemical master equations and exact simulation to burst distributions, linear noise, feedback, switching lifetimes, and growth-coupled memory.

CCBS · EST 5053 · Fall 2026 · lecture 5 · lecturer-written exposition · work through after the core lecture · see also the lineage

The core asks why a cell can have reproducible behavior and variable molecular histories. This exposition develops the tools needed to answer that question for a new network. A single addition–removal process provides our first complete calculation. Translation and promoter switching then reveal what this calibration leaves out. Finally, reaction-level covariance connects noise to Lecture 4's stability analysis, and growth connects a fixed-volume model to a dividing cell.

Parts 1–4 develop the local state, CME, event sampler and estimators. They include the conditions under which a long history represents a stationary ensemble. Part 5 derives the burst and distribution claims that the core introduces as extensions. Part 6 is a worked route through the linear-noise and cellular-memory problems. Part 7 asks which measurements could distinguish the mechanisms. These are study routes, not a second single-lecture timetable. Each checkpoint has an answer; try it before opening the solution.

We retain Lecture 3's reaction bookkeeping. Species are $X_1,\ldots,X_s$. Their random counts form $N_X(t)=(N_{X_1}(t),\ldots,N_{X_s}(t))^T$, and $n$ denotes one realized integer state. A species symbol is not a random variable. Named examples use shorter count symbols, such as $N$ for one species or $M,P$ for RNA and protein counts.

For channel $r$, reactant vector $\alpha_r$ and product vector $\beta_r$ list the numbers consumed and formed. The jump is $\gamma_r=\beta_r-\alpha_r$. Putting those columns together gives $\Gamma$. Propensities determine when those jumps occur.

1Begin with the question

A stochastic model is unnecessary until the observable asks for something a deterministic mean cannot answer.

Write the biological question before the reactions. “What is the mean expression level?” might be answered by a rate equation. “What fraction of cells have zero transcript?”, “How long until the first resistant cell appears?”, and “Does a threshold see isolated bursts?” cannot. They depend on the distribution, a boundary, or an individual history.

Now name the sampling operation. A snapshot across 1,000 cells is an ensemble. A 12-hour movie of one lineage is a time series. A population average after lysis is neither a single-cell distribution nor a path. The same mean can arise from all three while hiding different biology. Raj and colleagues made this distinction tangible by counting individual mRNAs and relating punctuated transcription to promoter states.⁷

Variation also has sources. Elowitz and colleagues used two fluorescent reporters in the same cell to separate reporter-specific fluctuations from shared cell-to-cell changes.⁵ Swain, Elowitz, and Siggia formalized that intrinsic-versus-extrinsic decomposition.⁶ A single broad histogram cannot perform this separation because measurement error, hidden cell state, and reaction-event variation are all superposed.

Choose a path, distribution or mean observable and state how it will be measured before selecting a model.

2Convert concentration to count

Concentration alone does not reveal discreteness. Volume completes the conversion.

For concentration $c$ in moles per liter and volume $V$ in liters, multiply by Avogadro's constant $N_A$ to get the molecule number $N$. Use the rough conversion $N_A\approx6\times10^{23}\,\mathrm{mol}^{-1}$:

A micromolar species in a femtoliter cell therefore has hundreds of copies. Invert the same relation: one copy corresponds to about $1/600\,\mu\mathrm M\approx1.7\,\mathrm{nM}$. The comparison tells us which concentrations lie near an integer boundary.

Events supply a second counting scale. For independent events with constant rate over a window, split the window into many intervals with a small firing probability. The binomial count approaches a Poisson law with mean and variance both equal to the expected event count $K$, so its relative fluctuation is $\sqrt K/K=K^{-1/2}$. A large molecule pool can still inherit a rare regulator's event history. Feedback, correlated packets and a changing rate require a different calculation.

Low count identifies a scale that needs further interpretation. Rare promoter states matter even when their protein output is abundant. Extinction matters when zero is absorbing. A system near a switching threshold may amplify small fluctuations. Conversely, a high-count species can inherit variation from a low-count regulator. Use the conversion to locate the possible source, then return to the observable.

Convert a concentration and compartment volume into a copy-number scale, and use that scale to test whether integer boundaries matter.

3Choose a local compartment

A geometric volume specifies where to count. A dynamical compartment also needs an exchange timescale and a named observable.

A dendritic spine gives a concrete volume

A dendritic spine is a small protrusion on a neuron's dendrite. The dendrite is an input-receiving branch. A spine head receives a synaptic input and connects to the branch through a narrow neck. The head's calcium signal can help regulate that synapse. Spine size and shape vary, so the following geometry is a representative model choice.²⁶

A submicrometer head has a subfemtoliter volume. Approximate the head by a sphere of radius $r_h\approx0.3\,\mu\mathrm m$. Its volume is

Thus the example volume is about $0.1\,\mathrm{fL}$. At 100 nM free calcium it contains about six free ions on average. At 10 nM it contains about 0.6. These means concern unbound cytoplasmic calcium. Calcium bound to a protein or stored in an organelle belongs to a different state variable. An instantaneous count remains integer even when its mean is below one.

Derive exchange through the neck

The head exchanges its contents through a restricted cross section. Let $c_h$ and $c_d$ be the head and dendrite concentrations of an ideal diffusing solute. Approximate the neck by a cylinder of length $L$ and area $A=\pi r_{\mathrm{neck}}^2$. Let $D$ be its diffusion coefficient. Assume a well-mixed head, a large dendritic reservoir and a nearly linear concentration profile along the neck.

Fick's law turns that concentration difference into an exchange rate. The gradient is approximately $(c_h-c_d)/L$. Multiply the diffusive flux density by the neck area, then divide by the head volume to obtain its concentration change:

The units check is useful: $DA/L$ is volume per time, so dividing $V$ by it gives time. A smaller neck radius lengthens equilibration as $r_{\mathrm{neck}}^{-2}$. The small head can mix internally while exchange with the dendrite remains slower.

Compare exchange with the process being observed

Four timescales answer different questions. A head can be a useful state variable while remaining coupled to its dendrite. Treating that state as approximately independent is a stronger simplification.

A calcium transient can remain local because clearance competes with spread. Let $x=c_h-c_{h,0}$ and $y=c_d-c_{d,0}$ denote deviations from a common resting concentration. Let $j(t)$ be an effective source in concentration per time. After buffering has been represented consistently in the coefficients, a linear compartment model is

The pulse divides between two removal routes. Suppose the dendrite stays at baseline, $y=0$, and an initial excess $x_0$ receives no further source. Define the relaxation rate $k_{\mathrm{relax}}=1/\tau_{\mathrm{clear}}+1/\tau_{\mathrm{exchange}}$. The solution is $x(t)=x_0e^{-k_{\mathrm{relax}}t}$. Integrating the neck loss gives the fraction that leaves by that route:

Local clearance dominates when $\tau_{\mathrm{clear}}\ll\tau_{\mathrm{exchange}}$. A sustained source can maintain a local signal through continuing influx and removal. The anatomical connection remains open.

A low-count pool can turn over rapidly

Mean count and renewal rate are independent pieces of information. A pool can retain the same mean while its molecules arrive and leave faster. In the constant-addition, linear-removal model developed below, multiplying both rates by the same factor leaves the stationary mean and count distribution unchanged. It shortens the correlation time by that factor. A time-integrating observer can then obtain more independent information from the same small mean pool.

The count conversion survives even when compartment independence fails. The relation $\bar N=cVN_A$ applies to any specified observation volume with the corresponding volume-averaged concentration. Predicting the fluctuations seen by a sensor additionally requires transport, binding, removal and the sensor's response time. A low count flags a scale worth modeling. It does not specify a probability law or its temporal correlations.

Define the local volume, molecule and readout. Derive or measure exchange, compare it with local reaction and mixing times, and keep the turnover rate separate from the instantaneous molecule count.

4Make hidden memory explicit

A Markov state contains enough present information to determine every instantaneous event rate.

The chemical master equation assumes that, conditional on state $N_X(t)=n$, the next-event rates do not depend on the earlier path. If a promoter that just closed has a different reopening hazard from one that has been closed for an hour, the on/off label alone is not Markov. Add an age, a sequence of conformational substates, or a non-Markov waiting-time model.

The same principle handles delays. Replacing transcription initiation by a reaction that instantly produces mature mRNA assumes elongation is either negligible or absorbed into an exponential effective step. If a fixed elongation delay shapes the autocorrelation, add an explicit chain or a delayed event. The simulator cannot repair missing state. It only samples the state model you wrote.

A deterministic decision can look random

Suppose I alternate trips between home and office, and carry my umbrella only if it is raining and the umbrella is at my departure location. Given location, umbrella location, and rain, my action is deterministic. A spectator who records only “carried / did not carry” has omitted the mechanism. That record need not itself be a Markov chain.

For a tractable model, assume rain on each departure is an independent Bernoulli event with probability $p$. Let the state before a trip be 0 if the umbrella is unavailable, 1 if available. From 0 I necessarily arrive where it is, so the next state is 1. From 1, rain makes me take it, leaving it available for my next departure; no rain leaves it behind. With current-state rows and next-state columns,

Under stationary sampling, the carrying probability is $p/(2-p)$, not $p$. The random input is our weather model; the decision rule remains deterministic. Biological intrinsic noise is likewise defined relative to a chosen system boundary, but molecular reaction events are not established to be deterministic merely because an everyday decision can be.

Add the hidden state needed to make event hazards depend on the present rather than an unrecorded history.

5Write jumps and propensities

Each reaction channel is an integer displacement paired with an event rate.

Keep the constant from Lecture 3 attached to its concentration law. For elementary channel $r$, write

The reactant order is $q_r=\sum_i\alpha_{ir}$. Since $v_r$ is an event flux in molarity per time, $k_r$ has units $\mathrm M^{1-q_r}/\mathrm{time}$. Stoichiometry acts afterward: the channel contributes $\gamma_r v_r$ to $\dot c$. In particular, $2A\xrightarrow{k}C$ with $v=kc_A^2$ contributes $-2kc_A^2$ to $\dot c_A$, not $-kc_A^2$.

Fix a well-mixed volume $V$ in liters and let $N_A$ be Avogadro's constant. Define $\Omega=N_AV$, so a count state $n$ has concentration $c=n/\Omega$. The factor $\Omega$ has units inverse molar. A reaction event changes concentration by $\gamma_r/\Omega$.

For an integer state $n\in\mathrm{Z}_{\ge0}^s$ containing the counts of $s$ species, channel $r$ has jump $\gamma_r$ and propensity $a_r(n)$. The definition concerns the identity of an event:

Two channels can share the same jump. Basal transcription and an activated route may both add one RNA; their destination probability contains the sum of their propensities. A channel rate is not a destination probability unless that destination uniquely identifies the channel.

For the elementary step $A+B\xrightarrow{k}C$, the jump is $(-1,-1,+1)$. Each of the $n_A$ molecules can partner with any of the $n_B$ molecules, giving $n_An_B$ eligible pairs. Multiply by a per-pair hazard to obtain the propensity.

For elementary $2A\xrightarrow{k}C$, choosing a first molecule and a distinct partner gives $n_A(n_A-1)$ ordered pairs. Each unordered pair appears twice. With three molecules there are six ordered pairs but only three physical pairs. If $\kappa$ is one physical pair's firing hazard,

Absorbing the factorial is a change of convention, not an approximation. However, $k_{\rm count}$ is not our concentration constant $k$. Its units are inverse time, whereas $k$ for this bimolecular channel has units inverse molar per time. We must determine their relationship before inserting a measured concentration-scale constant into an SSA.

Derive the conversion while keeping $k_r$ fixed

There are two ingredients in this conversion. Molecular counting determines the state dependence. The macroscopic limit determines the normalization of the per-combination hazard. Matching a deterministic equation alone would not determine the finite-count model.

For $A+B\xrightarrow{k}C$, the count law is $a=\kappa n_An_B$. Divide by $\Omega$ to get concentration event flux. Substituting $n_A=\Omega c_A$ and $n_B=\Omega c_B$ gives $a/\Omega=\kappa\Omega c_Ac_B$. Therefore $\kappa=k/\Omega$ matches $v=kc_Ac_B$.

For $2A\xrightarrow{k}C$, retain the unordered-pair count. The same calculation gives

At macroscopic volume with fixed nonzero concentration, $1/\Omega$ is negligible compared with $c_A$. Matching $v=kc_A^2$ requires $\kappa\Omega/2=k$. Thus

The combinatorial denominator cancels because the same chemical constant has a different conversion to an unordered-pair hazard. It does not justify setting $\kappa=k/\Omega$ for identical reactants. That would halve the macroscopic event flux. Gillespie gives the distinct-reactant and identical-reactant conversions explicitly in Section 2, pages 37–38.¹

The general formula counts each reactant species separately. Define $h_r(n)=\prod_i\binom{n_i}{\alpha_{ir}}$ as the number of unordered reactant combinations. Under the elementary stochastic mass-action assumption, each has the same hazard $\kappa_r$, so $a_r(n)=\kappa_r h_r(n)$. At large counts,

To recover the specified $v_r(c)$, set $\kappa_r=k_r\Omega^{1-q_r}\prod_i\alpha_{ir}!$. Define the falling factorial $(n)_m=n(n-1)\cdots(n-m+1)$, with $(n)_0=1$ and zero value when $m>n$. Substitution gives the exact propensity of the stated count model:

The factorial product is not a factorial of the total reaction order. A formal reactant vector $\alpha=(2,1)$ gives $h=\binom{n_A}{2}n_B$, with denominator $2!1!$, not $3!$. Its mass-action normalization would be $\kappa=2k/\Omega^2$. This is a combinatorial example, not evidence that a biological overall reaction is an elementary three-body collision. Binding and catalysis should be resolved at the mechanism's appropriate scale.

The conversion is exact. Replacing falling factorials by powers is not

At an admissible count state, the normalized propensity can be written exactly as

When every required reactant count is large compared with its stoichiometric coefficient, replacing each factor by $c_i$ recovers $v_r(c)$. For $2A$, the ratio $a/(\Omega v)=1-1/n_A$ for $n_A>0$. At one molecule the exact propensity is zero, while substituting concentration directly into $\Omega kc_A^2$ would permit an impossible event.

1. Order the state vector once and keep that order in every jump.
2. Compute products minus reactants for each component of $\gamma_r$.
3. Count eligible reactant combinations at the current integer state.
4. Multiply by a rate constant whose units make $a_r$ events per time.
5. Check that $n+\gamma_r$ stays admissible whenever $a_r(n)>0$.

Construct a jump-and-propensity table with legal destinations, eligible reactant combinations and consistent units.

6One small interval

The CME follows by keeping all events of order $dt$ and discarding events of smaller order.

Let $p(n,t)=\Pr\{N_X(t)=n\}$. To occupy $n$ at $t+dt$, either the system was already there and no channel fired, or it was at a predecessor $n-\gamma_r$ and channel $r$ fired once. The probability of two events is $o(dt)$. Therefore

Subtract $p(n,t)$, divide by $dt$, and take the limit. No Gaussian assumption enters. No random term is added to an ODE. The probability flow comes directly from discrete reaction events.

Account for all first-order event probabilities in one short interval before taking the continuous-time limit.

7The general CME

Gain from every predecessor minus loss through every outgoing channel.

Read it in this order. The argument $n-\gamma_r$ identifies the state that reaches $n$ after jump $\gamma_r$. The first propensity is evaluated at that predecessor because the event clock runs before the jump. The second term removes probability from $n$ at the rate its own clock runs.

The CME is linear in $p$ even when propensities are nonlinear in $n$. Its difficulty is dimension: a modest number of species creates a huge or infinite count lattice. Gillespie's direct method avoids enumerating that lattice by sampling paths instead.¹

Write a chemical master equation by evaluating each incoming clock at its predecessor and each outgoing clock at the current state.

8Boundaries and conservation

The admissible state space determines which predecessor terms exist.

Let $b$ be a constant addition propensity in events/time, and $d$ an individual removal hazard in inverse time. For constant addition and first-order removal,

At $n=0$, define $p_{-1}=0$. The removal clock is already zero. Hence $\dot p_0=dp_1-bp_0$. Sum equation (5) from zero to infinity. After shifting indices, every addition gain cancels an addition loss and every removal gain cancels a removal loss. Thus $d\sum_n p_n/dt=0$.

Run three checks on any CME: every gain points from a legal predecessor, every outgoing clock appears once with a minus sign, and summing over the legal state space conserves probability. If probability disappears, either an index is wrong or an absorbing/outside state has not been named.

Material conservation survives every random event

Probability conservation and molecular conservation are different statements. Resolve elementary binding as $A+B\xrightleftharpoons[k_-]{k_+}C_{AB}$. The two jumps are $(-1,-1,+1)^T$ and its negative. Each preserves $n_A+n_{C_{AB}}$ and $n_B+n_{C_{AB}}$, not merely their averages.

More generally, if a row $\ell^T$ satisfies $\ell^T\Gamma=0$, then $\ell^T N_X(t)$ is constant on every path of the closed reaction model. With fixed binding totals of 3 A-units and 5 B-units, choose complex count $z\in\{0,1,2,3\}$ as the only independent coordinate. The clocks become $a_+(z)=k_+(3-z)(5-z)/(N_AV)$ and $a_-(z)=k_-z$. At the upper boundary association vanishes. At zero dissociation vanishes. This is Lecture 3's conservation law used to reduce a stochastic state space.

Locate a probability-bookkeeping error by summing the CME, and distinguish an absorbing zero from a zero with a surviving source.

9The generator shortcut

The generator obtains observable equations without writing every state probability.

Choose a numerical observable $f(n)$ of the current count state, such as a count, its square, or an indicator that the cell is empty. We want its expected instantaneous change without solving all state probabilities. Condition on $N_X(t)=n$. In a short interval $h$, reaction $r$ changes the observable by $f(n+\gamma_r)-f(n)$ with probability $a_r(n)h+o(h)$. No event contributes zero increment. Hence

Divide by $h$ and let it vanish. The resulting linear operator is the infinitesimal generator $\mathcal L$. It takes an observable to another function of state, with units observable/time:

Then $d\langle f(N_X)\rangle/dt=\langle(\mathcal{L}f)(N_X)\rangle$. Choosing $f(n)=n_i$ gives the mean of $N_{X_i}$. Choosing $f(n)=n_in_j$ gives the second moment $\mathrm E[N_{X_i}N_{X_j}]$. Each term asks only how one event changes the selected observable.

For nonlinear propensities, the hierarchy usually opens. A mean equation contains a second moment, which contains a third. Moment closure then becomes an approximation. For affine propensities, polynomial moments of a given order close at that order. The constant-addition, linear-removal process and the two-stage gene-expression example exploit this special structure.

Identical reactants reveal two separate finite-count corrections

Consider only the channel $2A\xrightarrow{k}C$ from Section 5. Let $N(t)$ denote its random $A$ count, $\bar c_A=\mathrm E[N]/\Omega$ its mean concentration, and $\sigma_c^2=\operatorname{Var}(N)/\Omega^2$ its concentration variance. The jump is $-2$, so the exact mean equation is

Use $\mathrm E[N(N-1)]=(\mathrm E[N])^2+\operatorname{Var}(N)-\mathrm E[N]$, then divide by $\Omega$:

The variance term comes from averaging a nonlinear propensity. The negative $\bar c_A/\Omega$ term comes from excluding a molecule as its own partner. Both must be accounted for when comparing with $\dot c_A=-2kc_A^2$. Neither changes the meaning of $k$.

A Poisson count has $\operatorname{Var}(N)=\mathrm E[N]$, so these two corrections cancel at that instant. The dimerization process does not in general preserve a Poisson distribution. An initial agreement therefore does not establish exact deterministic mean dynamics for all time.

A matrix example fixes the transpose convention

A promoter opens with hazard $\alpha$ and closes with hazard $\beta$. In off/on order, using a column of probabilities,

The columns of $Q$ sum to zero. For the on-state indicator $f=(0,1)^T$, $Q^Tf=(\alpha,-\beta)^T$: off gains indicator at rate $\alpha$, on loses it at rate $\beta$. The on-probability therefore obeys $\dot p_1=\alpha-(\alpha+\beta)p_1$. Define its stationary value $p_{\mathrm{on}}=\alpha/(\alpha+\beta)$; integration gives

The umbrella matrix used discrete-time row probabilities; this example uses continuous-time probability columns. State the convention before comparing matrices.

Derive an observable equation from conditional event increments, and identify the moment closure needed before using a nonlinear deterministic mean.

10Constant addition and linear removal

This standard immigration–death process is valuable because every representation is tractable.

The state is $n\ge0$. The jump and propensity table is:

The deterministic concentration analogue is addition minus first-order removal. The stochastic model adds an integer boundary and a distribution, but it does not change the mechanism. This makes disagreement among the representations especially diagnostic.

Use the addition and removal clocks to identify the stationary count scale and the correct zero-state equation.

11Poisson, step by step

Stationary neighboring fluxes generate a recursion, and normalization finishes it.

Assume $b,d>0$. Define current across edge $n\leftrightarrow n+1$ of the state graph by $J_n=bp_n-d(n+1)p_{n+1}$. The CME says $\dot p_n=J_{n-1}-J_n$. At stationarity the current is constant along the ladder. The zero boundary gives $0=\dot p_0=-J_0$, so that constant is zero. Only now conclude $bp_n=d(n+1)p_{n+1}$. Rearrange and define $\mu=b/d$:

Iterate: $p_1=\mu p_0$, $p_2=\mu^2p_0/2!$, and in general $p_n=p_0\mu^n/n!$. Normalize:

Thus $N\sim\operatorname{Poisson}(\mu)$. Its mean and variance both equal $\mu$. This is the reference process against which additional gene-expression clocks can be diagnosed.

Zero edge current is stronger than stationary probability. A cyclic state network can sustain circulating current with every state probability constant. Our boundary argument licenses detailed balance for this nearest-neighbor ladder; it does not say cells at steady state are at thermodynamic equilibrium.

Obtain the Poisson law from boundary-forced zero current and normalization, without assuming pairwise balance on an arbitrary network.

12Moments, step by step

Apply the generator to $n$ and $n^2$, then convert the second moment to variance.

For $f(n)=n$, addition changes $f$ by 1 and removal changes it by −1. Equation (6) gives

For $f(n)=n^2$, an addition changes the square by $(n+1)^2-n^2=2n+1$; a removal changes it by $(n-1)^2-n^2=-2n+1$. Therefore

Use $\sigma^2=\langle N^2\rangle-\langle N\rangle^2$, reserving $V$ for volume. Differentiation subtracts $2\langle N\rangle(b-d\langle N\rangle)$ from equation (10); the terms $2b\langle N\rangle$ cancel, and the second moments combine into variance:

At stationarity, equation (9) gives $\langle N\rangle=b/d$. Equation (11) then gives $\sigma^2=b/d$. The probability recursion and moment method agree, which is precisely why this network is a good implementation test.

Derive both mean and variance from event increments, including the subtraction of the evolving squared mean.

13Transients and system size

The calibration model lets us separate four effects: initialization, fluctuation memory, system size and observation duration.

Starting from mean $M_0$, equation (9) solves to

The relaxation time is $1/d$. Sampling at one tenth of that time does not produce a stationary ensemble. Initialization matters even when the eventual stationary law is known.

The entire transient law has a direct construction. Start with exactly $n_0$ molecules. Each independently survives with probability $e^{-dt}$, giving a binomial survivor count. New arrivals form a Poisson process; an arrival at time $u$ survives with probability $e^{-d(t-u)}$. Independent thinning gives a Poisson new-survivor count with mean $\int_0^t b e^{-d(t-u)}du$. The contributions are independent:

For $n_0=0$, the transient law is Poisson with a changing mean. At $b=4$, $d=0.2$ and $t=60$ minutes, its mean is $20(1-e^{-12})\simeq19.99988$. The sampled mean from 5,000 paths has standard error about $\sqrt{20/5000}=0.063$.

The stable state still fluctuates

Lecture 4's linearized matrix is the scalar $A=-d$, restoring mean perturbations as $e^{-dt}$. Events continually inject variance. At the stationary mean their squared-jump-weighted rate is $D=(+1)^2b+(-1)^2d(b/d)=2b$. If $C$ is stationary variance, equation (11) becomes

This algebraic Lyapunov equation balances covariance; it is not itself a nonlinear global-stability proof. The same structure appears in network linear-noise approximations; here the second moment is exact. Also, $\langle N(t+\tau)\mid N(t)\rangle=b/d+[N(t)-b/d]e^{-d\tau}$. Multiply by the centered count at time $t$ and average at stationarity: $\operatorname{Cov}[N(t),N(t+\tau)]=(b/d)e^{-d|\tau|}$. The restoring time is also the fluctuation memory time.

Change system size without changing the concentration kinetics. Fix a reference volume $V_{\mathrm{ref}}$ and define the dimensionless multiplier $R_V=V/V_{\mathrm{ref}}$. The physical conversion factor remains $\Omega=N_AV=R_VN_AV_{\mathrm{ref}}$. If the reference addition propensity is $b$, use $R_Vb$ in the enlarged volume and keep the per-molecule removal hazard $d$ unchanged. Stationary count mean grows as $R_V$, standard deviation as $\sqrt{R_V}$, and relative spread as $R_V^{-1/2}$. Kurtz gives the rigorous density-dependent convergence to deterministic kinetics.³

The plotted normalization is count per reference volume. Write $z=N/R_V$. Then $\mathrm E[z]=b/d$ and $\operatorname{Var}(z)=b/(dR_V)$ at stationarity. Actual molarity is $c=z/(N_AV_{\mathrm{ref}})$. This separates a dimensionless size comparison from the unit conversion used to derive propensities.

Match initial conditions and observation windows before comparing a transient ensemble with a stationary distribution or system-size limit.

14Derive the waiting time

While no event occurs, the state is fixed, so the total hazard is constant.

Let $a_0(n)=\sum_r a_r(n)$. If $S(\tau)$ is the probability that no event occurs for another $\tau$, then

Hence $S(\tau)=e^{-a_0\tau}$. If $U_1$ is uniform on $(0,1)$, solve $U_1=S(\tau)$:

This skips every empty interval. It is especially efficient when events are rare because the sampled waiting time can be arbitrarily long without repeated null steps.

Sample a next-event time from the survival probability and state why the hazards must remain fixed while waiting.

15Derive the winning channel

Conditional on an event, each channel owns a fraction of the total hazard.

The joint density that the next event occurs at time $\tau$ and is channel $j$ is $a_j e^{-a_0\tau}$. Integrating over time yields $\Pr(J=j)=a_j/a_0$. Draw a second uniform $U_2$ and choose the first $j$ satisfying

The order of channels changes no probability if the intervals have the correct widths. After the chosen jump, every propensity that depends on a changed species must be recomputed. The clocks were conditionally independent only while the old state remained fixed.

Select the firing channel from cumulative propensity and update the state only after identifying that channel.

16Implement the direct method

A correct simulator is a short loop surrounded by explicit assertions and event logs.

1. Initialize time $t$, integer state $n$, stop condition, and a reproducible random seed.
2. Evaluate every $a_r(n)$. Assert that each value is finite and nonnegative.
3. Set $a_0=\sum_ra_r$. If it is zero, stop at an absorbing state.
4. Draw $U_1,U_2\in(0,1)$. Set $\tau=-\ln U_1/a_0$.
5. If $t+\tau$ passes the observation time, record the unchanged state there and stop.
6. Choose the channel by equation (15), advance time, and apply $n\leftarrow n+\gamma_j$.
7. Assert integer nonnegative state, log time and channel, then return to step 2.

Store events rather than only regularly sampled output. A step path can be rendered later at any observation grid. Event times are also needed to test waiting-time models, while a regularly sampled trajectory may hide short episodes.

A complete reference implementation

The following Python uses only the standard library. Each returned pair records a state immediately after a jump; the last pair is an observation-time sentinel, not an event. The state between successive timestamps is constant. Parameters use minutes, with $b$ in additions/minute and $d$ in inverse minutes.

import math
import random

def addition_removal(b=4.0, d=0.2, end=60.0, n0=0, seed=123):
    assert all(math.isfinite(z) and z >= 0 for z in (b, d, end))
    assert isinstance(n0, int) and n0 >= 0
    rng = random.Random(seed)
    t, n = 0.0, n0
    history = [(t, n)]
    while t < end:
        total = b + d*n
        if total == 0:
            break
        u = rng.random()
        while u == 0:  # require an open lower endpoint for the logarithm
            u = rng.random()
        next_t = t - math.log(u)/total
        if next_t >= end:
            break
        n += 1 if rng.random()*total < b else -1
        assert n >= 0
        t = next_t
        history.append((t, n))
    history.append((end, n))
    return history

The horizon test must occur before applying the jump. Otherwise an event from the future contaminates the endpoint sample. At zero, removal has zero propensity; clipping a negative count would hide a broken event selector. A zero-total-hazard state still contributes its entire remaining residence time.