Stochastic chemical kinetics

Stochastic chemical kinetics models molecular copy numbers as a continuous-time Markov chain. Each reaction has a state-dependent reaction propensity function and changes the state by its stoichiometric vector when it fires.

This noise-relevant event size weights absolute stoichiometric jump sizes by their molecular turnover flux. At stationarity the denominator is twice the balanced production or removal flux $v_i$, so the normalized diffusion diagonal is $D_{ii}=2⟨s_i⟩/({\tau }_i{\mu }_i)$. Fixed birth bursts of size $b$ and unit deaths give $⟨s⟩=(b+1)/2$. This quantity is not the event-count-weighted average absolute jump size.

// Target: mathematical-biology.bigb

At stationarity, the mean number ${\mu }_i$ divided by the molecular removal flux $v_i$ is the mean time a molecule remains in the system. Birth or death event rates must be multiplied by the numbers of molecules changed in each event before computing this lifetime.

// Target: mathematical-biology.bigb

Compartment-based stochastic diffusion represents molecular motion by nearest-neighbour jumps on a spatial lattice. Jump rates of order $D/h^2\pm v/(2h)$ converge to an advection--diffusion equation as compartment width $h$ tends to zero.

The asymmetric simple exclusion process permits at most one particle per lattice site and biases jumps in one direction. A mean-field hydrodynamic limit has advective flux proportional to $c(1-c)$ because both an occupied departure site and vacant arrival site are required.

An immigration--death process has constant birth rate $b$ and linear death rate $dy$ in state $y$. Its stationary distribution is Poisson with mean $b/d$.

A population jumps from $n$ to $n+2$ at a constant immigration rate $\lambda$ and from $n$ to $n-1$ at rate $\beta n$. Applying its Markov jump-process generator to $n$ and $n^2$ gives the displayed mean and variance equations. For $\beta >0$ and finite initial second moment, the stationary mean and variance are $2\lambda /\beta$ and $3\lambda /\beta$. Interpreting the pair-arrival rate as proportional to population instead would produce a different branching process.

A stochastic quasi-steady-state approximation replaces fast species by their stationary conditional distribution given the slow species. Averaging slow propensities over that distribution yields a reduced Markov jump process.

A slow-scale stochastic simulation algorithm applies the Gillespie algorithm to propensities averaged over fast conditional equilibrium, avoiding explicit simulation of every fast reaction.

The molecular copy number of a chemical species is the nonnegative integer number of its molecules in a specified reactor or compartment.

Genomic copy number is the count of copies of a specified DNA sequence in a defined biological unit, such as a cell or a haploid-equivalent genome. The sequence and denominator matter: a two-copy autosomal locus in a diploid cell is one copy per haploid genome, while a deletion or duplication can alter that locus without changing the ploidy of the whole cell. Such segment-level differences underlie copy number variation.^[ref]

DNA replication also changes the number of DNA molecules during the cell cycle. A replicated chromosome contains sister copies without thereby becoming two independently counted chromosomes under the usual centromere-based convention. Report the cell-cycle stage when comparing absolute DNA copies per cell.

Sequencing depth can estimate relative genomic copy number by comparing coverage with a reference region, but library preparation, sequence-dependent mapping and population mixtures affect that comparison. In a hypothetical mixture where half the cells carry two copies and half carry four, the mean is three copies per cell; that average does not show that any individual cell carries three. Genomic copy number also differs from the number of expressed RNA or protein molecules. This distinction matters when connecting gene dosage to a molecular abundance model.

The reaction propensity $a_r(x)$ is the instantaneous rate at which reaction $r$ fires when the copy-number state is $x$. Conditional on the present state, its reaction clock has an exponential distribution of rate $a_r(x)$.

A power-law propensity uses the same reactant monomial as a deterministic law of mass action, scaled by reactor volume. For example, a bimolecular channel may be assigned $a(x)=\alpha x^2/V$. This convention differs at small copy number from the combinatorial propensity $\alpha x(x-1)/V$ and must therefore be stated explicitly.

For reactions with stoichiometric vectors ${\nu }_r$ and reaction propensity functions $a_r$, the probability mass function $p(x,t)$ obeysIt is the forward operator of a Markov jump process equation for the copy-number continuous-time Markov chain.

The Gillespie algorithm samples statistically exact trajectories of a well-mixed chemical reaction network modeled by a continuous-time Markov chain. At copy-number state $x$, each reaction propensity function $a_r(x)$ is a nonnegative event rate and each stoichiometric vector ${\nu }_r$ specifies the update. Assume the rates have no explicit time dependence. Put $a_0(x)={\sum }_r{a}_r(x)$. If $a_0(x)=0$, no further chemical reaction can occur: the state is an absorbing state.

For $a_0(x)>0$, sample a waiting time $\tau$ from the exponential distribution with rate $a_0(x)$, independently choose channel $r$ with probability $a_r(x)/a_0(x)$, and update $(t,x)\leftarrow (t+\tau ,x+{\nu }_r)$. Recompute the reaction propensity functions and repeat. The no-event probability is $S(\tau )=e^{-a_0\tau }$ because $S^{\prime }=-a_{0}S$; multiplying by $a_r$ gives the joint waiting-time/channel probability density function $a_r{e}^{-a_0\tau }$. This factorization explains the two draws.^[ref]

For $A\to B$ with reaction propensity function $a(x)=k{x}_A$, the waiting time has mean $1/(k{x}_A)$, and $x_A=0$ is an absorbing state. Exactness concerns the specified chemical master equation; spatial mixing errors and explicitly time-varying rates require a different model or sampling rule.

Tau-leaping accelerates the Gillespie algorithm by approximating reaction rates as constant over an interval $\tau$, then drawing all reaction counts for that interval at once. At state $x$, with reaction propensity functions $a_r(x)$ and stoichiometric vectors ${\nu }_r$, the explicit step isThe leap condition requires the propensities to change little during the interval. Unlike event-by-event exact simulation, this introduces an approximation error controlled through the choice of $\tau$.^[ref]

For $A\to B$ with rate $k{x}_A$, the approximate number converted is Poisson with mean $k{x}_A\tau$. The exact number converted over the interval is instead binomial with parameters $x_A$ and $1-e^{-k\tau }$. The two agree to first order for small $k\tau$, but the Poisson draw can exceed $x_A$ and produce a negative count. Smaller or adaptive leaps and modified sampling rules address this limitation; freezing the rate is not automatically safe at low copy number.

A moment hierarchy is the coupled family of equations obtained by applying a Markov jump-process generator to powers of the state. Nonlinear reaction propensity functions make a moment of one order depend on higher-order moments, so the hierarchy generally does not close after finitely many equations.

For a Markov jump-process generator $L$, differentiating an integrable observable gives $d\mathbb{E}[f(X_t)]/dt=\mathbb{E}[Lf(X_t)]$, when the generator identity and interchange are justified. Choosing powers gives equations for moments. These need not form a closed finite system.

A moment closure replaces higher moments in a moment hierarchy by functions of retained lower moments. It converts an infinite hierarchy into a finite approximate system.

A central-moment closure sets selected high-order central moments to zero. For example, setting the third central moment of $X$ to zero gives