Reaction Kinetics in Biology

Full Derivation from Mass Action

An enzyme $E$ binds substrate $S$ to form a complex $ES$ which then converts to product $P$:

Mass action kinetics for each species:

The total enzyme is conserved: $[E]_0 = [E] + [ES]$, so $[E] = [E]_0 - [ES]$. The steady-state approximation assumes $d[ES]/dt = 0$(valid when $[S] \gg [E]_0$):

Solving for $[ES]$:

where we define the Michaelis constant:

The rate of product formation is $v = k_2[ES]$, giving the Michaelis-Menten equation:

Key properties: $v \to V_{\max}$ as $[S] \to \infty$ (saturation),$v = V_{\max}/2$ when $[S] = K_M$, and$v \approx (V_{\max}/K_M)[S]$ for $[S] \ll K_M$ (first-order regime). The ratio $k_{\text{cat}}/K_M$ is the catalytic efficiency, bounded above by the diffusion limit $\sim 10^9$ M$^{-1}$s$^{-1}$.

The Hill Equation

The phenomenological Hill equation describes cooperative binding:

where $n$ is the Hill coefficient:$n > 1$ indicates positive cooperativity (sigmoidal curve),$n = 1$ is Michaelis-Menten, and $n < 1$ is negative cooperativity. For hemoglobin binding oxygen: $n \approx 2.8$ (4 binding sites, not fully cooperative).

MWC (Monod-Wyman-Changeux) Concerted Model

The MWC model provides a mechanistic basis for cooperativity. The protein exists in two conformational states: R (relaxed, high affinity) and T (tense, low affinity). All subunits switch together (concerted transition). Key parameters:

• • $L = [T_0]/[R_0]$: allosteric constant (equilibrium between T and R with no ligand)
• • $K_R$: dissociation constant for ligand binding to R state
• • $K_T$: dissociation constant for ligand binding to T state
• • $c = K_R/K_T$: ratio of affinities ($c < 1$ means R binds more tightly)

For a protein with $n$ identical subunits, the partition functions for the R and T states with $k$ ligands bound are:

Defining $\alpha = [S]/K_R$, the fractional saturation is:

The effective Hill coefficient at the midpoint can be derived from the MWC parameters:

For large $L$ and small $c$ (strong allosteric effect), $n_H \to n$. The MWC model explains how hemoglobin achieves $n_H \approx 2.8$ with 4 subunits:$L \approx 10^4{-}10^6$, $c \approx 0.01$.

Mutual Repression and Bistability

Two genes that mutually repress each other can create a bistable switch (Gardner et al. 2000). The dynamics are:

where $\alpha$ is the maximum expression rate and $n$ is the Hill coefficient of repression. The nullclines are curves where each derivative vanishes:

Fixed points occur at nullcline intersections. By symmetry, there is always one fixed point at $[A] = [B]$. For $n > 1$ and sufficiently large $\alpha/\gamma$, two additional stable fixed points appear (one with A high/B low, one with A low/B high), making the symmetric point unstable.

Bifurcation Analysis

To find the bifurcation point, we analyze the stability of the symmetric fixed point$[A]^* = [B]^* = A^*$ where $A^* = \alpha/(\gamma(1 + A^{*n}))$. Linearizing around this point with perturbations $\delta A$, $\delta B$:

Let $f = \frac{n\alpha A^{*n-1}}{(1+A^{*n})^2}$. The eigenvalues are$\lambda_{\pm} = -\gamma \pm f$. The symmetric point becomes unstable when$\lambda_+ > 0$, i.e., when $f > \gamma$:

This is a pitchfork bifurcation. The critical Hill coefficient for bistability depends on $\alpha/\gamma$: for the symmetric case, $n > 1$ is necessary but not sufficient. As $\alpha/\gamma$ increases, bistability emerges for smaller$n$ values.

The Chemical Master Equation

For a system with discrete molecule counts $\mathbf{n} = (n_1, n_2, \ldots)$, the probability $P(\mathbf{n}, t)$ evolves according to the chemical master equation:

where $a_\mu(\mathbf{n})$ is the propensity (rate) of reaction $\mu$ and$\boldsymbol{\nu}_\mu$ is the stoichiometry vector. For a simple birth-death process ($\emptyset \xrightarrow{k} X$, $X \xrightarrow{\gamma} \emptyset$):

The steady-state solution is a Poisson distribution:

For a Poisson distribution, $\text{Var}(n) = \langle n \rangle$, so the Fano factor $F = \text{Var}/\text{Mean} = 1$. Deviations from $F = 1$ reveal additional sources of noise (e.g., transcriptional bursting gives $F = 1 + b$ where $b$ is the burst size).

The Gillespie Algorithm

The Gillespie (stochastic simulation) algorithm generates exact sample paths of the master equation. The derivation proceeds as follows. Given the current state $\mathbf{n}$ at time $t$:

Step 1: Compute all propensities $a_\mu(\mathbf{n})$and the total propensity $a_0 = \sum_\mu a_\mu$. The probability that no reaction occurs in $[t, t+\tau)$ and then reaction $\mu$ fires in $[t+\tau, t+\tau+d\tau)$ is:

Step 2: Sample the waiting time from the exponential distribution:$\tau = -\ln(r_1)/a_0$ where $r_1 \sim U(0,1)$.

Step 3: Choose which reaction fires: select $\mu$ such that$\sum_{j=1}^{\mu-1} a_j < r_2 a_0 \leq \sum_{j=1}^{\mu} a_j$ where $r_2 \sim U(0,1)$.

Step 4: Update state: $\mathbf{n} \to \mathbf{n} + \boldsymbol{\nu}_\mu$,$t \to t + \tau$.

Stochastic effects matter when molecule numbers are small. The coefficient of variation scales as $\text{CV} \sim 1/\sqrt{N}$, so for $N < 100$ molecules, fluctuations exceed 10% of the mean and deterministic ODEs become unreliable.

The Repressilator Model

The repressilator (Elowitz & Leibler, 2000) consists of three genes forming a cyclic negative feedback loop: gene 1 represses gene 2, gene 2 represses gene 3, and gene 3 represses gene 1. The dynamics are:

where $j = ((i-2) \bmod 3) + 1$ gives the cyclic repression, $\alpha$ is the maximum transcription rate, $\alpha_0$ is the leakiness, $\beta$ is the translation rate, and $\delta_m, \delta_p$ are the degradation rates.

Conditions for Oscillation: Hopf Bifurcation

The system has a symmetric fixed point where all proteins and mRNAs are equal: $m^* = m_i$,$p^* = p_i$ for all $i$. To determine stability, we linearize around this fixed point. The 6D Jacobian has a circulant block structure that can be diagonalized using the discrete Fourier transform with roots of unity $\omega_k = e^{2\pi i k/3}$,$k = 0, 1, 2$.

For each mode $k$, the eigenvalues satisfy the 2D characteristic equation:

where $f'(p^*) = \frac{n\alpha p^{*n-1}}{(1+p^{*n})^2}$ is the repression sensitivity. The $k=0$ mode (uniform perturbation) is always stable. The $k=1$ and $k=2$modes (complex conjugate pair) give oscillatory instability. Setting $\lambda = i\omega$and separating real and imaginary parts, the Hopf bifurcation criterion is:

The oscillation frequency at the bifurcation point is:

This condition requires sufficiently strong repression (large $\alpha$ and $n$) relative to the degradation rates. For the parameters used by Elowitz & Leibler ($n \approx 2$, strong promoters), the repressilator robustly oscillates with a period of several hours, consistent with the cell division cycle.