Chemical reaction network

A chemical reaction network specifies chemical species, reactions between them, and their stoichiometric vectors. It can be modelled deterministically by reaction-rate equations or stochastically by molecule-copy-number jumps.

At a positive concentration vector $c$, the total incoming deterministic reaction flux equals the total outgoing flux at every reaction complex. Here $c^y={\prod }_i{c}_i^{y_i}$ is the law of mass action monomial. This is stronger than balancing net production of each chemical species, but does not require pairwise detailed balance between reverse chemical reactions. A directed cycle can have nonzero circulation while satisfying complex balance.

A formal combination of chemical species appearing on one side of a chemical reaction. The coefficients $y_i$ are nonnegative integers. A chemical reaction network can be viewed as a directed graph whose vertices are reaction complexes and whose edges are chemical reactions. The empty complex represents external production or removal.

Equal constant production, individual linear removal, and joint bimolecular removal give $d({\mu }_1-{\mu }_2)/dt=-\beta ({\mu }_1-{\mu }_2)$. For $\beta >0$ the unique stable equilibrium is symmetric and solves $\lambda =\beta \mu +C{\mu }^2$. For $\beta =0<C$, every point on ${\mu }_1{\mu }_2=\lambda /C$ is an equilibrium: the difference is conserved and trajectories approach the equilibrium on their own constant-difference line.

// Target: mathematical-biology.bigb

If individual removal is absent, simultaneous removal of one molecule from each species leaves $X_1-X_2$ unchanged. Independent constant-rate births change it by plus or minus one, each at rate $\lambda$. The difference is therefore a symmetric continuous-time random walk on the integers with variance growing as displayed. It has no normalizable stationary law. This explains why a deterministic attracting line of equilibria does not provide a stationary stochastic distribution.

// Target: probability-and-statistics.bigb

For $E=C\mu /(\beta +C\mu )<1$ and mean lifetime $\tau =(\beta +C\mu {)}^{-1}$, the normalized restoring matrix is $M={\tau }^{-1}\left(\begin{array}{cc}1& E\\ E& 1\end{array}\right)$ and the diffusion matrix is $D=(\tau \mu {)}^{-1}\left(\begin{array}{cc}2& E\\ E& 2\end{array}\right)$. Solving the normalized stationary fluctuation-dissipation relation for a reaction network gives the displayed covariance. The sum mode remains damped, but the difference mode's restoring rate is $\beta =(1-E)/\tau$. Its noise grows without bound relative to that restoring rate as $E\to 1$, giving strong negative correlations and invalidating small-fluctuation closure near exclusive joint removal.

// Target: mathematical-biology.bigb

A chemical species is a class of chemically identical entities treated as one component of a reaction system.

$\mathrm{R}{\mathrm{O}}^{-}$, the conjugate base of an alcohol.

The conjugate base of a propanol; the straight-chain form is CH3CH2CH2O−. A specified isomer matters in product structures.

A lithium salt of ethoxide. Its reactive oxygen donor and lithium-centred acceptor can participate in complementary contacts with an ester; solvation and aggregation affect a detailed orbital description.

$(\mathrm{C}{\mathrm{H}}_3{)}_3\mathrm{C}{\mathrm{O}}^{-}$, a bulky alkoxide base, often used in elimination reactions.

$\mathrm{C}{\mathrm{H}}_3\mathrm{C}{\mathrm{H}}_2{\mathrm{O}}^{-}$, the conjugate base of ethanol.

$\mathrm{C}{\mathrm{H}}_3{\mathrm{O}}^{-}$, the conjugate base of methanol.

The mole fraction of species $i$ is $x_i=N_i/{\sum }_j{N}_j$, where $N_i$ is its amount in moles. It is dimensionless and equals the volume mixing ratio for an ideal gas mixture.

For an ideal gas mixture, the volume mixing ratio $x_i$ equals the mole fraction, number fraction, and ratio of partial to total pressure: $x_i=N_i/N=n_i/n=P_i/P$, with ${\sum }_i{x}_i=1$. This equivalence fails for general nonideal mixtures.

The stoichiometric vector of reaction $r$ records its net change in each chemical species. If the copy-number state is $x$, firing the reaction changes it to $x+{\nu }_r$.

A reaction-rate equation is an ordinary differential equation for continuously varying species concentrations. Under the law of mass action, its reaction rates are monomials in those concentrations.

The displayed logarithmic derivative compares the loss and production flux sensitivities of species $i$ to the mean amount ${\mu }_j$. At a steady state $v_i^{+}=v_i^{-}=v_i$, it converts the normalized negative drift Jacobian into $M_{ij}=H_{ij}/{\tau }_i$, where ${\tau }_i={\mu }_i/v_i$ is the mean molecular lifetime. For constant production and loss proportional to ${\mu }^{\alpha }$, the scalar elasticity is $H=\alpha$.

// Target: mathematical-biology.bigb

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