Sources and References

Primary textbook — B. Ingalls, Mathematical Modeling in Systems Biology (MIT Press, 2013; freely available at math.uwaterloo.ca/~bingalls/)

Supplementary texts — U. Alon, An Introduction to Systems Biology: Design Principles of Biological Circuits (2nd ed., CRC Press, 2019); J.D. Murray, Mathematical Biology I: An Introduction (3rd ed., Springer, 2002); J. Keener & J. Sneyd, Mathematical Physiology (Springer, 2009)

Online resources — XPPAUT documentation (math.pitt.edu/~bard/xpp); BioModels Database (ebi.ac.uk/biomodels); Virtual Cell (vcell.org)

---

Chapter 1: Biochemical Kinetics and Reaction Networks

1.1 Mass Action Kinetics

Living cells are dynamic chemical systems — a human cell contains thousands of distinct molecular species participating in a web of interacting reactions. To reason quantitatively about cell behaviour, we model these reactions using the framework of chemical kinetics. The fundamental postulate of chemical kinetics is the law of mass action: the rate of an elementary reaction is proportional to the product of the concentrations of its reactants.

For a bimolecular reaction

\[ A + B \xrightarrow{k} C, \]

the reaction rate is \(v = k[A][B]\), where \(k\) is the rate constant (units depend on the order of the reaction). For a first-order, unimolecular reaction \(A \xrightarrow{k} B\), the rate is simply \(v = k[A]\). The intuition is elementary: two molecules react only when they collide, and the collision frequency is proportional to the product of their concentrations.

1.2 ODEs from Reaction Networks

Given a set of reactions and their rates, we construct ordinary differential equations (ODEs) for each molecular species by tracking how every reaction contributes to production or consumption. A key application is computing the steady state — the concentration vector at which all rates of change vanish. For a linear pathway with mass action kinetics, this can be done by hand.

Consider the reversible reaction

\[ A \underset{k_{-1}}{\stackrel{k_1}{\rightleftharpoons}} B. \]

The forward reaction \(A \to B\) occurs at rate \(k_1[A]\) and the reverse reaction \(B \to A\) at rate \(k_{-1}[B]\). Therefore:

\[ \frac{d[A]}{dt} = -k_1[A] + k_{-1}[B], \qquad \frac{d[B]}{dt} = k_1[A] - k_{-1}[B]. \]

At thermodynamic equilibrium, the net rates vanish: \(k_1[A]^* = k_{-1}[B]^*\), giving the equilibrium constant

\[ K_{eq} = \frac{[B]^*}{[A]^*} = \frac{k_1}{k_{-1}}. \]

This ratio of rate constants encodes the thermodynamic favourability of the reaction. The principle of detailed balance states that at equilibrium each individual reaction must be balanced by its reverse — this connects kinetic parameters to thermodynamics and constrains the parameters of any physically consistent model.

1.3 The Stoichiometric Matrix

For networks with many species and reactions, tracking each ODE separately becomes unwieldy. The stoichiometric matrix \(S\) provides a compact, systematic framework. We define:

• Each column of \(S\) corresponds to a reaction.
• Each row of \(S\) corresponds to a molecular species.
• The entry \(S_{ij}\) equals the net number of molecules of species \(i\) produced per occurrence of reaction \(j\) (negative for consumption, positive for production).

If \(\mathbf{X}\) is the vector of species concentrations and \(\mathbf{v}(\mathbf{X})\) is the vector of reaction rates, then the entire ODE system is compactly expressed as:

\[ \frac{d\mathbf{X}}{dt} = S \cdot \mathbf{v}(\mathbf{X}). \]

1.4 A More Complex Network: Michaelis-Menten via Stoichiometric Matrix

1.5 Conservation Laws and Moiety Conservation

Many biochemical reaction networks conserve the total amount of certain chemical moieties (structural subgroups), even as individual species interconvert. The systematic enumeration of conservation laws via the left null space of \(S\) is not merely a mathematical convenience — it directly guides model reduction and experimental design. A conserved quantity is both a constraint that reduces the dimensionality of the ODE system and a biological invariant that can be measured to verify model consistency. These conservation laws reduce the effective dimension of the system. Mathematically, a conservation law corresponds to a vector \(\mathbf{c}\) such that \(\mathbf{c}^T S = \mathbf{0}\) — that is, \(\mathbf{c}\) lies in the left null space of \(S\). Then

\[ \frac{d}{dt}\left(\mathbf{c}^T \mathbf{X}\right) = \mathbf{c}^T S \mathbf{v} = 0, \]

so the quantity \(\mathbf{c}^T \mathbf{X}\) is constant in time.

A ubiquitous example is enzyme conservation: a free enzyme \(E\) and its substrate-bound complex \(ES\) together conserve the total enzyme:

\[ [E] + [ES] = E_{tot} = \text{const}. \]

This allows us to eliminate one variable and reduce the system’s dimension. Conservation laws are automatically found by computing the left null space of \(S\), a linear algebra operation readily performed by software.

This theorem is the foundation of flux balance analysis: rather than solving nonlinear kinetic equations, one characterises all feasible steady-state fluxes as a cone and optimises within it.

---

Chapter 2: Enzyme Kinetics

2.1 Motivation and the Michaelis–Menten Scheme

Nearly all biochemical reactions in cells are catalysed by enzymes — protein catalysts that dramatically accelerate reaction rates without being consumed. Understanding enzyme kinetics is central to modelling cellular metabolism, signalling, and gene regulation.

The canonical model is the Michaelis–Menten scheme, proposed by Leonor Michaelis and Maud Menten in 1913. A substrate \(S\) binds reversibly to the enzyme \(E\) to form a complex \(ES\), which then releases product \(P\) and regenerates the free enzyme:

\[ E + S \underset{k_{-1}}{\stackrel{k_1}{\rightleftharpoons}} ES \xrightarrow{k_2} E + P. \]

By mass action kinetics, the ODEs for this scheme are:

\[ \frac{d[S]}{dt} = -k_1[E][S] + k_{-1}[ES], \]\[ \frac{d[ES]}{dt} = k_1[E][S] - k_{-1}[ES] - k_2[ES], \]\[ \frac{d[E]}{dt} = -k_1[E][S] + (k_{-1} + k_2)[ES], \]\[ \frac{d[P]}{dt} = k_2[ES]. \]

The enzyme conservation \([E] + [ES] = E_{tot}\) allows us to eliminate \([E]\). This still leaves a nonlinear system, which motivates the key approximation.

2.2 The Quasi-Steady-State Approximation

If the enzyme concentration is much smaller than the substrate concentration (\(E_{tot} \ll [S]\)), the complex \([ES]\) reaches a quasi-steady state rapidly and thereafter tracks the slow variable \([S]\). Setting \(d[ES]/dt \approx 0\):

\[ k_1(E_{tot} - [ES])[S] - (k_{-1} + k_2)[ES] = 0. \]

Solving for \([ES]\):

\[ [ES] = \frac{E_{tot}[S]}{K_M + [S]}, \]

where \(K_M = (k_{-1} + k_2)/k_1\) is the Michaelis constant. The net rate of product formation is then:

\[ v = k_2[ES] = \frac{V_{max}[S]}{K_M + [S]}, \]

with \(V_{max} = k_2 E_{tot}\). This is the celebrated Michaelis–Menten equation.

The parameter \(K_M\) has an elegant biological interpretation: it measures the enzyme’s affinity for the substrate (lower \(K_M\) means higher affinity). When \([S] \ll K_M\), the rate is approximately linear in \([S]\): \(v \approx (V_{max}/K_M)[S]\), resembling first-order kinetics. When \([S] \gg K_M\), the enzyme is saturated and \(v \approx V_{max}\), giving zeroth-order kinetics.

2.3 Validity of the Quasi-Steady-State Approximation

The QSSA is valid when the complex \([ES]\) equilibrates fast relative to the substrate. A more precise criterion (Segel and Slemrod, 1989) is:

\[ \frac{E_{tot}}{K_M + [S]_0} \ll 1, \]

where \([S]_0\) is the initial substrate concentration. This corrects the naive criterion \(E_{tot} \ll [S]_0\) by including the case when \([S]_0\) is itself small but \(K_M\) is large.

2.4 Enzyme Inhibition

Inhibitors are molecules that reduce enzyme activity. Their mechanisms fall into two main categories.

The distinction has profound practical implications. Competitive inhibitors (such as many drugs targeting enzyme active sites) can be titrated out by high substrate concentrations — the cell can in principle compensate. Non-competitive inhibitors are more potent regulatory tools because their effect is independent of substrate level.

2.4 Cooperativity and the Hill Equation

Many enzymes and receptors display cooperativity: binding of one ligand molecule alters the affinity of subsequent ligand molecules. A classic example is haemoglobin, where binding of one oxygen molecule increases affinity for subsequent molecules (positive cooperativity).

The Hill equation phenomenologically captures sigmoidal binding curves:

\[ v = \frac{V_{max}[S]^n}{K_{half}^n + [S]^n}, \]

where \(n\) is the Hill coefficient and \(K_{half}\) is the concentration giving half-maximal activity. For \(n = 1\), the Hill equation reduces to Michaelis–Menten kinetics. For \(n > 1\), the response is sigmoidal — the curve rises steeply around \(K_{half}\), creating a switch-like response that is highly biologically relevant. Hill coefficients \(n > 1\) indicate positive cooperativity; \(n < 1\) indicates negative cooperativity.

The sigmoidal response from cooperativity is a core design principle in cell biology: it allows biochemical switches that respond sharply to small changes in ligand concentration, filtering noise while amplifying biologically meaningful signals.

---

Chapter 3: Dynamics and Stability Analysis

3.1 Phase Plane Analysis

For a two-dimensional ODE system

\[ \frac{dx}{dt} = f(x, y), \qquad \frac{dy}{dt} = g(x, y), \]

the phase plane is the \((x, y)\) space in which the trajectory of a solution is traced. A central tool is the nullcline: the \(x\)-nullcline is the set of points where \(dx/dt = 0\) (i.e., \(f(x,y) = 0\)), and the \(y\)-nullcline where \(dy/dt = 0\). Fixed points (steady states) occur at the intersections of nullclines.

By drawing nullclines and noting the direction of flow in each region, one can often determine qualitative behaviour without solving equations explicitly. This geometric approach is particularly powerful for biological systems, where precise parameter values may be unknown.

3.2 Linearisation and the Jacobian

Near a fixed point \(\mathbf{x}^*\), the dynamics of small perturbations \(\delta\mathbf{x} = \mathbf{x} - \mathbf{x}^*\) are governed by the linearised system:

\[ \frac{d\,\delta\mathbf{x}}{dt} \approx J\,\delta\mathbf{x}, \]

where \(J\) is the Jacobian matrix evaluated at \(\mathbf{x}^*\):

\[ J_{ij} = \frac{\partial f_i}{\partial x_j}\bigg|_{\mathbf{x}^*}. \]

The stability of the fixed point is determined by the eigenvalues \(\lambda\) of \(J\). If all eigenvalues satisfy \(\text{Re}(\lambda) < 0\), the fixed point is locally asymptotically stable: perturbations decay and trajectories return to \(\mathbf{x}^*\). If any eigenvalue has \(\text{Re}(\lambda) > 0\), the fixed point is unstable.

3.3 Worked Linearisation: The Predator–Prey System

Before discussing bifurcations, we work through a complete linearisation example that illustrates the general method.

3.4 Bifurcation Theory

A bifurcation occurs when a qualitative change in system behaviour (number or stability of fixed points, appearance of a limit cycle) happens as a parameter is continuously varied. Bifurcation analysis reveals how biological systems switch between functional states.

Saddle-node bifurcation: As a parameter \(\mu\) increases through a critical value \(\mu_c\), two fixed points (one stable, one unstable) collide and annihilate. For \(\mu < \mu_c\) there are two fixed points; for \(\mu > \mu_c\) there are none. This mechanism underlies the sudden appearance or disappearance of stable states in biological switches.

The Hopf bifurcation is the primary mechanism generating sustained oscillations in biological models — circadian rhythms, cell cycle oscillations, and neural firing all arise through Hopf bifurcations.

3.5 Limit Cycles and the Goodwin Oscillator

A limit cycle is an isolated closed trajectory in the phase plane: nearby trajectories spiral toward it (stable limit cycle) or away from it (unstable limit cycle). Unlike the closed orbits of the Lotka–Volterra system (which form a continuous family parameterised by the conserved \(H\)), a limit cycle is structurally stable — it persists under small perturbations of the system.

The Goodwin oscillator is the canonical minimal model for a biological negative-feedback oscillator — a single gene producing a protein that ultimately represses its own synthesis through a chain of intermediate steps. With \(m\) intermediate species \(X_1, \ldots, X_m\) and the last species repressing production of the first:

\[ \frac{dX_1}{dt} = \frac{v_1}{1 + (X_m/K)^n} - k_1 X_1, \]\[ \frac{dX_i}{dt} = k_{i-1} X_{i-1} - k_i X_i, \quad i = 2, \ldots, m. \]

Each step degrades the previous species and produces the next, with Hill-function repression closing the loop.

3.7 Bistability

A system is bistable if it has two coexisting stable steady states. In the phase plane, this typically appears as an S-shaped (sigmoidal) nullcline intersecting another nullcline at three points: two stable nodes separated by an unstable saddle. The two stable states correspond to distinct biological phenotypes, and the saddle defines the threshold (separatrix) between the two basins of attraction.

Bistability requires some form of positive feedback or cooperative nonlinearity. Once a bistable system commits to one state, it remains there even if the stimulus that triggered the switch is removed — the system exhibits hysteresis as a parameter sweeps back and forth. This irreversibility is fundamental to biological decision-making (cell differentiation, the mitotic switch).

3.8 Separation of Time Scales

When a system contains processes occurring at very different rates, one can exploit separation of time scales to reduce model complexity. Fast variables relax quickly to a manifold determined by the slow variables; one can then apply a quasi-steady-state approximation for the fast variables and obtain a lower-dimensional system for the slow dynamics. This is the mathematical justification for the Michaelis–Menten QSSA and underlies many model reductions in cell biology.

---

Chapter 4: Metabolic Networks

4.1 Metabolic Pathways and Feedback Inhibition

Metabolic networks consist of chains and cycles of enzyme-catalysed reactions that convert nutrient substrates into energy and biosynthetic precursors. A canonical motif is the linear biosynthetic pathway:

\[ S \xrightarrow{E_1} M_1 \xrightarrow{E_2} M_2 \xrightarrow{E_3} P, \]

where \(S\) is the substrate, \(M_i\) are intermediates, and \(P\) is the end product. A common regulatory design is end-product (feedback) inhibition: the end product \(P\) inhibits the first enzyme \(E_1\), preventing wasteful overproduction of \(P\) when its concentration is already high. This is a classic example of negative feedback securing homeostasis — maintaining product levels near a set point despite fluctuating inputs.

The inhibition is often allosteric: \(P\) binds a regulatory site on \(E_1\), inducing a conformational change that reduces activity. The Hill equation captures the sigmoidal dependence of this inhibition on \([P]\).

4.2 Metabolic Control Analysis

A central question in metabolic network modelling is: how is control over the steady-state metabolic flux distributed among the enzymes? Metabolic Control Analysis (MCA) provides a quantitative framework through dimensionless sensitivity coefficients.

The flux control coefficient measures the relative change in steady-state flux \(J\) in response to a relative change in the activity (or level) of enzyme \(E_i\):

\[ C^J_{E_i} = \frac{\partial J}{\partial E_i} \cdot \frac{E_i}{J} = \frac{d\ln J}{d\ln E_i}. \]

A coefficient \(C^J_{E_i} = 0.8\) means that a 1% increase in enzyme \(i\) increases steady-state flux by 0.8% — enzyme \(i\) exerts substantial control.

This theorem, independently discovered by Burns, Kacser, and Heinrich and Rapoport in the early 1970s, is more surprising than it looks. The naive biochemist’s intuition is that one enzyme — the “rate-limiting step” — controls the flux, with all others playing a secondary role. If that were true, the sum would be dominated by a single coefficient near 1, with the rest near 0. What Kacser and Burns showed experimentally is that control is typically distributed: in many real pathways, no single enzyme has \( C^J_{E_i} > 0.5 \), and the sum constraint forces the coefficients to compete. The theorem thus quantitatively refutes the “rate-limiting step” dogma and suggests that interventions (drug targets, metabolic engineering) should account for systemic redistribution of control, not just the single bottleneck enzyme. It is a conservation law for regulatory influence.

Elasticity coefficients describe how an individual enzyme’s rate responds to changes in its substrates or products:

\[ \varepsilon^{v_i}_{S} = \frac{\partial \ln v_i}{\partial \ln [S]}. \]

These are local properties of individual enzymes (measured in isolation), while flux control coefficients are global properties of the network. The connectivity theorem relates the two:

\[ \sum_i C^J_{E_i} \varepsilon^{v_i}_{S} = 0, \]

for any intermediate metabolite \(S\). MCA thus provides a powerful framework for understanding how enzymatic parameters determine the behaviour of an entire metabolic system.

4.3 Flux Control Coefficients: A Worked Example

The summation theorem is most illuminating through a concrete calculation.

4.4 Flux Balance Analysis

For large-scale metabolic networks (thousands of reactions), the full kinetic details are unknown. Flux balance analysis (FBA) bypasses kinetics by imposing the steady-state constraint \(S \cdot \mathbf{v} = \mathbf{0}\) and then using linear programming to find the flux distribution \(\mathbf{v}\) that maximises a biological objective (typically growth or ATP production), subject to thermodynamic and capacity constraints. FBA has been remarkably successful at predicting growth phenotypes and the effects of gene knockouts in bacteria and yeast.

---

Chapter 5: Signal Transduction

5.1 Receptor–Ligand Binding

Cells communicate through signal transduction: extracellular signals (hormones, growth factors, neurotransmitters) bind to cell-surface receptors, triggering intracellular cascades. The basic binding reaction is

\[ R + L \underset{k_{-1}}{\stackrel{k_1}{\rightleftharpoons}} RL, \]

where \(R\) is the free receptor and \(L\) is the ligand. At steady state, the fraction of occupied receptors follows a hyperbolic (Michaelis–Menten-like) curve:

\[ f = \frac{[L]}{K_D + [L]}, \]

where \(K_D = k_{-1}/k_1\) is the dissociation constant (lower \(K_D\) means higher affinity). This is the dose-response relationship: it describes how the cellular response varies with ligand concentration. The steepness of the dose-response curve — its sensitivity — determines how precisely a cell can detect changes in ligand concentration.

5.2 Dose-Response: A Numerical Example

5.3 Signalling Cascades and the MAPK Cascade

Many signalling pathways employ multi-tier kinase cascades. The MAPK (MAP kinase) cascade is a three-tier phosphorylation cascade that amplifies upstream signals and is involved in cell proliferation, differentiation, and stress responses. At each tier, a kinase is activated (by phosphorylation) and in turn activates the next kinase, while phosphatases oppose each activation.

A remarkable property discovered by Goldbeter and Koshland is zero-order ultrasensitivity. When the kinase and phosphatase operating on a substrate both work near saturation (zero-order kinetics), the steady-state phosphorylation level is an extremely steep sigmoid function of the kinase-to-phosphatase ratio. The effective Hill coefficient of this response can be much greater than 1 even for non-cooperative enzymes.

The steady-state fraction of activated substrate \(W^*\) at each tier satisfies the Goldbeter–Koshland equation. If we denote the kinase activity \(u\), phosphatase activity \(v\), and Michaelis constants \(J = K_M^{kinase}/[S_{total}]\) and \(K = K_M^{phosphatase}/[S_{total}]\), then

\[ W^* = GK(u, v, J, K) = \frac{2u J}{v - u + vJ + uK + \sqrt{(v - u + vJ + uK)^2 - 4(v-u)uJ}}. \]

This function transitions sharply from 0 to 1 as the ratio \(u/v\) crosses 1, with the sharpness increasing as \(J\) and \(K\) decrease (i.e., as both enzymes approach saturation). When three such zero-order ultrasensitive stages are cascaded in the MAPK pathway, the overall input-output response can achieve an apparent Hill coefficient exceeding 30 — far more switch-like than anything achievable by allosteric cooperativity alone.

5.4 Adaptation

Many signalling systems display perfect adaptation: after a sustained step change in stimulus, the response returns exactly to its pre-stimulus level, encoding only changes in stimulus rather than absolute level. This is analogous to derivative sensing.

Perfect adaptation requires integral feedback control — a special network topology in which the output feeds back through an integrating node. Mathematically, perfect adaptation requires specific network motifs: a regulatory circuit that effectively computes the time integral of the output error. The analysis of adaptation in bacterial chemotaxis (Barkai and Leibler, 1997) showed that near-perfect adaptation can be achieved in a robust, parameter-insensitive manner through a methylation-based integral feedback.

5.5 Signalling Network Motifs

Biological signalling networks are not random wiring diagrams — they are built from recurring network motifs, small subgraphs that appear far more frequently than expected by chance. Uri Alon’s systematic analysis of the transcription network of E. coli identified the feedforward loop (FFL) as one of the most common motifs. In a coherent type-1 FFL, a transcription factor X activates both Y and Z, and Y also activates Z. The temporal logic of this circuit functions as a sign-sensitive delay: Z is activated only after a sustained (not transient) stimulus to X, while Z turns off rapidly when X is inactivated. This filtering of transient signals is biologically important — it prevents cellular responses to noise while still responding to genuine sustained signals.

The mathematical analysis of FFLs uses first-order linear ODEs: \(Z' = \beta_Z - \gamma_Z Z\) (when X and Y are both active), and one shows that the time to half-maximal Z expression is \(t_{1/2} = \ln(2)/\gamma_Z\), independent of the exact parameter values. The delay is thus robust.

5.6 Bistability and Oscillations in Signalling

Positive feedback in signalling can generate bistability, creating all-or-nothing switches. The activation of MPF (Maturation Promoting Factor) during mitotic entry is bistable: a cell either commits to division or stays in interphase, with hysteresis ensuring the commitment is irreversible.

Negative feedback with time delay generates sustained oscillations. The NF-κB signalling system oscillates due to negative feedback through the inhibitor IκBα, which is transcriptionally induced by NF-κB itself but requires time for synthesis and translocation. Similarly, calcium oscillations arise from the interplay of IP3-mediated calcium release and calcium-dependent inactivation — a system amenable to Hopf bifurcation analysis.

5.6 Gene Regulatory Networks: The Toggle Switch

Gene expression is regulated by transcription factors — proteins that bind specific DNA sequences and either activate or repress transcription. Networks of transcription factors implement computation at the cellular level.

\[ u\text{-nullcline: } u = \frac{\alpha_1}{1 + v^n}, \qquad v\text{-nullcline: } v = \frac{\alpha_2}{1 + u^n}. \]

For small \(n\) (weak cooperativity) and moderate \(\alpha_1 = \alpha_2\), the two nullclines intersect only once, at the symmetric fixed point \(u^* = v^*\) — the system is monostable. For large \(n\) (strong cooperativity) or large \(\alpha\), the nullclines develop an S-shaped character and intersect at three points: two stable nodes and one unstable saddle. The system is then bistable: it can reside in either “gene \(u\) on / gene \(v\) off” or “gene \(u\) off / gene \(v\) on”, corresponding to distinct cellular phenotypes.

---

Chapter 6: Gene Regulatory Networks

6.1 Transcription, Translation, and the Basic Gene Expression Model

Gene expression involves two main steps: transcription (synthesis of mRNA from a DNA template) and translation (synthesis of protein from the mRNA template). A minimal model treats mRNA concentration \(m\) and protein concentration \(p\) as dynamic variables:

\[ \frac{dm}{dt} = \beta_m - \gamma_m m, \qquad \frac{dp}{dt} = \beta_p m - \gamma_p p, \]

where \(\beta_m\) is the mRNA transcription rate, \(\gamma_m\) the mRNA degradation rate constant, \(\beta_p\) the translation rate constant (protein produced per mRNA per time), and \(\gamma_p\) the protein degradation rate constant. The steady-state protein level is:

\[ p^* = \frac{\beta_m \beta_p}{\gamma_m \gamma_p}. \]

The timescales are set by the half-lives: \(\tau_m = \ln 2 / \gamma_m\) and \(\tau_p = \ln 2 / \gamma_p\). In many bacteria, \(\tau_m \approx\) 2–5 minutes while \(\tau_p\) can be hours (protein is diluted by cell division as well as degraded). The separation of these timescales justifies a further quasi-steady-state approximation: setting \(m \approx \beta_m/\gamma_m\) gives an effective single-equation model for protein.

6.2 The Toggle Switch

The toggle switch (Gardner et al., Nature 2000) is a synthetic gene circuit comprising two mutually repressing genes. Gene A encodes a repressor of gene B, and gene B encodes a repressor of gene A. With Hill-function repression:

\[ \frac{dm_A}{dt} = \frac{\alpha_A}{1 + \left([P_B]/K_B\right)^n} - \gamma_{m_A} m_A, \]\[ \frac{dm_B}{dt} = \frac{\alpha_B}{1 + \left([P_A]/K_A\right)^n} - \gamma_{m_B} m_B. \]

Phase plane analysis reveals that for sufficiently strong repression (large Hill coefficient \(n\) or high repressor production), this network is bistable: either gene A is highly expressed and gene B suppressed, or vice versa. The two stable states correspond to distinct cellular phenotypes. Transient stimuli can flip the circuit from one state to the other. The toggle switch illustrates how bistability emerges from positive feedback (mutual repression creates a double-negative, equivalent to positive, feedback loop).

6.3 The Repressilator

The repressilator (Elowitz & Leibler, Nature 2000) is a three-gene synthetic circuit where gene A represses gene B, gene B represses gene C, and gene C represses gene A — forming a ring of three repressors. This negative feedback loop with an odd number of elements (here, three) generates sustained oscillations, as confirmed experimentally via GFP fluorescence.

The minimal protein-only model:

\[ \frac{dp_i}{dt} = \frac{\alpha}{1 + p_j^n} - \gamma p_i, \quad (i, j) \in \{(A,C), (B,A), (C,B)\}. \]

For large \(n\) (strong cooperativity), the system undergoes a Hopf bifurcation as the repression strength increases, giving rise to limit cycle oscillations. The repressilator demonstrates that oscillations do not require explicit positive feedback — they can arise from delayed negative feedback alone.

6.4 Noise in Gene Expression

In single cells, many proteins are present at low copy numbers (tens to hundreds of molecules), and biochemical reactions are inherently stochastic. Intrinsic noise arises from random fluctuations in the timing of individual transcription, translation, and degradation events — even genetically identical cells in identical environments will show cell-to-cell variability. Extrinsic noise arises from fluctuations in upstream components (transcription factors, ribosomes, polymerases) shared by many genes.

The Gillespie algorithm provides an exact method for simulating the stochastic dynamics of a well-mixed chemical system. At each step, two random numbers determine (i) the time to the next reaction event and (ii) which reaction fires, with probabilities proportional to the reaction propensities. For large molecule numbers, the stochastic simulations converge to the deterministic ODE predictions, but for low copy numbers, stochastic effects dominate.

For a simple gene expression model (birth–death process for protein), the stationary distribution of protein number is Poisson with mean \(\langle p \rangle = \beta/\gamma\) — intrinsic noise satisfies \(CV^2 = \sigma^2/\langle p\rangle^2 = 1/\langle p\rangle\), so noise is inversely proportional to mean protein number.

6.5 Circadian Oscillators

Circadian oscillators are biochemical clocks that generate ~24-hour oscillations in cellular activity, enabling organisms to anticipate daily environmental cycles. The simplest mathematical model capturing the core mechanism is the Goodwin oscillator: a negative feedback loop with a single gene producing a protein that represses its own transcription, with multiple intermediate steps introducing sufficient delay.

The Drosophila circadian oscillator involves the Period (PER) and Timeless (TIM) proteins, which dimerize and translocate to the nucleus to repress their own transcription. A key kinetic feature is the multi-step phosphorylation of PER (catalysed by the kinase DBT), which introduces the delay needed for oscillations. The mathematical analysis involves finding the Hopf bifurcation boundary in parameter space as the number of phosphorylation steps and the repression strength are varied.

---

Chapter 7: Electrophysiology

7.1 Membrane Potential and Ion Channels

Neurons and muscle cells are electrically excitable: they can generate rapid, transient changes in their transmembrane voltage — action potentials — that propagate along cell membranes to transmit information. The transmembrane voltage \(V\) arises from differences in ion concentrations across the lipid bilayer, maintained by active transporters (ion pumps).

The electrochemical equilibrium potential for a single ion species is given by the Nernst equation:

\[ E_{ion} = \frac{RT}{zF} \ln\frac{[ion]_{out}}{[ion]_{in}}, \]

where \(R\) is the gas constant, \(T\) temperature (Kelvin), \(F\) the Faraday constant, and \(z\) the ion valence. For \(K^+\) at body temperature with typical physiological concentrations, \(E_K \approx -90\) mV; for \(Na^+\), \(E_{Na} \approx +60\) mV. The resting membrane potential of a typical neuron (~−70 mV) lies between these extremes because multiple ion species contribute, weighted by their respective membrane permeabilities. The Goldman equation accounts for this:

\[ V_{rest} = \frac{RT}{F} \ln\frac{P_K[K^+]_{out} + P_{Na}[Na^+]_{out} + P_{Cl}[Cl^-]_{in}}{P_K[K^+]_{in} + P_{Na}[Na^+]_{in} + P_{Cl}[Cl^-]_{out}}, \]

where \(P_X\) denotes the membrane permeability to ion species \(X\). Ion flow across the membrane is mediated by ion channels — protein pores that open and close in response to voltage, ligands, or mechanical stimuli.

7.2 The Hodgkin–Huxley Model

The Hodgkin–Huxley (HH) model (1952, Nobel Prize 1963) is the quantitative foundation of modern computational neuroscience. Alan Hodgkin and Andrew Huxley developed the model from voltage-clamp experiments on the squid giant axon. The membrane is modelled as a capacitor in parallel with conductance elements for each ion species:

\[ C_m \frac{dV}{dt} = I_{app} - I_{Na} - I_K - I_L, \]

where \(C_m\) is the membrane capacitance, \(I_{app}\) an applied current, and the ionic currents are:

\[ I_{Na} = g_{Na} m^3 h (V - E_{Na}), \]\[ I_K = g_K n^4 (V - E_K), \]\[ I_L = g_L (V - E_L). \]

Here \(g_{Na}\), \(g_K\), \(g_L\) are maximal conductances, and \(m\), \(h\), \(n\) are gating variables between 0 and 1 representing the probabilities that individual channel gates are open. Each gating variable obeys a first-order kinetic equation:

\[ \frac{dw}{dt} = \frac{w_\infty(V) - w}{\tau_w(V)}, \quad w \in \{m, h, n\}, \]

where \(w_\infty(V)\) is the voltage-dependent steady-state gate value and \(\tau_w(V)\) is the voltage-dependent time constant, both determined from experimental data.

The action potential mechanism in the HH model: A sufficiently large depolarisation (membrane potential increasing from rest) rapidly opens Na\(^+\) channels (fast \(m\) activation), causing an inward Na\(^+\) current that further depolarises the membrane — a positive feedback loop driving the rapid upstroke. The rising voltage simultaneously begins to inactivate Na\(^+\) channels (slow \(h\) inactivation) and activate K\(^+\) channels (slow \(n\) activation), both of which drive repolarisation back toward rest. The refractory period following an action potential arises because the \(h\) gates remain inactivated for a time, preventing immediate re-firing.

7.3 The FitzHugh–Nagumo Model

The full 4-dimensional HH model is difficult to analyse analytically. The FitzHugh–Nagumo (FHN) model is a 2-dimensional caricature that retains the essential excitability:

\[ \frac{dv}{dt} = v - \frac{v^3}{3} - w + I_{app}, \]\[ \frac{dw}{dt} = \varepsilon(v + a - bw), \]

where \(v\) is a fast voltage-like variable, \(w\) is a slow recovery variable, \(I_{app}\) is the applied current, and \(\varepsilon \ll 1\) enforces the time-scale separation.

The \(v\)-nullcline is a cubic \(w = v - v^3/3 + I_{app}\); the \(w\)-nullcline is the line \(w = (v + a)/b\). Their intersection defines the fixed point. For small \(I_{app}\), the system is at a stable resting state (excitable: a sufficiently large perturbation initiates a large excursion before returning to rest). As \(I_{app}\) increases, a subcritical Hopf bifurcation occurs and the fixed point becomes unstable, giving repetitive firing (limit cycle). The phase plane of FHN is analytically tractable and reveals the threshold separatrix, refractory period, and oscillatory regime in a transparent geometric way.

7.4 FitzHugh-Nagumo: A Worked Bifurcation Analysis

7.5 Bursting Neurons

Many neurons do not fire regular spike trains but produce bursts — clusters of action potentials separated by quiescent periods. Bursting can be understood through slow-fast analysis: the fast subsystem (voltage and fast gates) is HH-like, while a slow variable (e.g., slow K\(^+\) current or intracellular Ca\(^{2+}\)) modulates the fast subsystem’s bifurcation parameter on a slower time scale.

As the slow variable evolves, it drives the fast subsystem through a bifurcation: when the fast subsystem is in its oscillatory regime, spikes fire (the burst); when the slow variable has accumulated enough to push the fast subsystem back to quiescence (via a fold bifurcation or homoclinic bifurcation), firing stops until the slow variable recovers and the next burst begins. This geometrical decomposition of slow-fast dynamics is central to the computational analysis of neuronal diversity.

7.6 Synaptic Transmission

Neurons communicate at synapses. An action potential in the presynaptic neuron causes Ca\(^{2+}\)-mediated release of neurotransmitter into the synaptic cleft. Neurotransmitter binds postsynaptic receptors, opening ion channels and generating a postsynaptic potential: an EPSP (excitatory postsynaptic potential) if the reversal potential drives depolarisation, or an IPSP (inhibitory postsynaptic potential) if it drives hyperpolarisation.

The simplest model of a postsynaptic neuron is the integrate-and-fire model:

\[ C_m \frac{dV}{dt} = -g_L(V - E_L) + I_{syn}(t), \]

with the rule that whenever \(V\) reaches a threshold \(V_{th}\), a spike is recorded and \(V\) is reset to \(V_{reset}\). Despite its simplicity, the integrate-and-fire model captures many features of neural coding and is widely used in network-scale simulations.

---

Chapter 8: Bifurcation Theory — Normal Forms and Biological Switching

8.1 Motivation: How Biological Systems Change Behaviour

Many of the most striking phenomena in cell biology — cell differentiation, the onset of oscillations in a circadian clock, the transition from interphase to mitosis — are not gradual changes but abrupt qualitative transitions. A system that was monostable suddenly becomes bistable; a steady state that was attracting suddenly becomes repelling and a limit cycle appears. Bifurcation theory is the mathematical study of these qualitative transitions: it catalogues the distinct ways that the phase portrait of a dynamical system can change as a parameter varies, and provides normal forms — minimal canonical models — that capture the universal behaviour near each transition.

The central insight is that near a bifurcation point, all systems of the same codimension and type look identical after an appropriate change of coordinates. This means that a few universal normal forms capture an enormous diversity of biological phenomena, and that qualitative predictions (the system will oscillate, or will be bistable) can be made on topological grounds without knowing precise parameter values.

8.2 Saddle-Node Bifurcation

The saddle-node bifurcation is the simplest and most ubiquitous bifurcation. It describes the creation or annihilation of a pair of fixed points — one stable, one unstable — as a parameter varies.

The bifurcation diagram for the normal form \(\dot{x} = \mu - x^2\) has a characteristic parabolic shape in the \((\mu, x)\) plane: the curve of fixed points is \(x = \pm\sqrt{\mu}\), existing only for \(\mu \geq 0\). The upper branch (\(x_+ > 0\)) is unstable (filled circles in a bifurcation diagram) and the lower branch (\(x_- < 0\)) is stable (open circles). At \(\mu = 0\), the two branches collide and annihilate — a fold in the bifurcation curve. Moving the parameter from \(\mu > 0\) to \(\mu < 0\) causes the stable steady state to disappear suddenly, potentially triggering a large-amplitude excursion to a distant attractor. This catastrophic jump is the mathematical signature of a tipping point.

8.3 Pitchfork Bifurcation

When a system has a symmetry (such as \(x \to -x\)), saddle-node bifurcations occur in symmetric pairs and the resulting bifurcation is the pitchfork. This is the archetypal bifurcation for systems with a symmetry that can be spontaneously broken.

The supercritical pitchfork models spontaneous symmetry breaking: for \(\mu < 0\), the system rests at the symmetric state \(x = 0\); as \(\mu\) increases through 0, the symmetric state becomes unstable and the system settles into one of the two symmetry-broken states \(\pm\sqrt{\mu}\). Which one it chooses is determined by noise or initial conditions.

8.4 Hopf Bifurcation — Normal Form and Worked Example

The Hopf bifurcation governs the birth of limit cycle oscillations from a fixed point. The normal form in polar coordinates \((r, \theta)\) for the supercritical case is:

\[ \dot{r} = \mu r - r^3, \qquad \dot{\theta} = \omega, \]

where \(r \geq 0\) is the amplitude of oscillation and \(\omega > 0\) is the angular frequency. For \(\mu < 0\), the only attractor is the fixed point \(r = 0\) (the origin, a stable spiral). For \(\mu > 0\), the fixed point \(r = 0\) is unstable and a stable limit cycle appears at \(r^* = \sqrt{\mu}\) with angular frequency \(\omega\). At \(\mu = 0\), the eigenvalues of the linearisation are exactly \(\pm i\omega\).

---

Chapter 9: Delay Differential Equations

9.1 Motivation: Time Delays in Biology

Many biological processes involve inherent time delays that cannot be adequately captured by ordinary differential equations. Gene regulation is a canonical example: after a transcription factor binds a promoter, a cascade of events must occur — RNA polymerase recruitment, transcription, mRNA processing, nuclear export, ribosome assembly, and translation — before the encoded protein appears in the cytoplasm. This total delay \(\tau\) can be tens of minutes. If the protein product then feeds back to repress its own gene, the time-delayed negative feedback can generate oscillations that a corresponding ODE model (without delay) would not predict.

Other examples include: immune response delays (cytokine production after pathogen recognition takes hours), population dynamics (gestation periods in predator-prey models), and pharmacokinetics (drug absorption and distribution times). The mathematical framework for these systems is the delay differential equation (DDE).

9.2 The Simplest DDE: Linear Stability Analysis

The canonical scalar DDE with negative feedback is:

\[ \dot{x}(t) = -k\, x(t - \tau), \]

where \(k > 0\) is the feedback gain and \(\tau > 0\) is the time delay. The unique steady state is \(x^* = 0\). To analyse stability, we seek solutions of the form \(x(t) = e^{\lambda t}\) (a characteristic solution). Substituting:

\[ \lambda e^{\lambda t} = -k e^{\lambda(t-\tau)} = -k e^{-\lambda\tau} e^{\lambda t}. \]

Dividing by \(e^{\lambda t} \neq 0\) gives the characteristic equation:

\[ \lambda = -k e^{-\lambda\tau}. \]

This transcendental equation has infinitely many complex roots \(\lambda_n\), in contrast to a polynomial characteristic equation from an ODE. For stability, we need all roots to satisfy \(\text{Re}(\lambda) < 0\).

---

Chapter 10: Stochastic Models and Intrinsic Noise

10.1 Why Stochasticity Matters in Biology

The deterministic ODE framework assumes that molecular concentrations are continuous variables evolving smoothly. In reality, molecules are discrete entities, and at the low copy numbers typical of gene regulatory proteins (tens to hundreds of molecules per cell), individual molecular events — a single transcription factor binding its promoter, a single mRNA being translated — can cause significant fluctuations. These fluctuations are intrinsic noise: they arise from the stochastic nature of chemical reactions even in a perfectly controlled environment.

Extrinsic noise, by contrast, arises from cell-to-cell variability in the concentrations of upstream regulators (RNA polymerases, ribosomes, signalling proteins). Both sources of noise are observable in single-cell experiments (fluorescence microscopy of individual cells) and contribute to the heterogeneity of cell populations — even genetically identical cells in identical conditions exhibit different protein levels.

Stochastic models are essential when: (1) copy numbers are low, (2) the system is near a bifurcation point (where noise can trigger transitions between attractors), (3) the timing of rare events (e.g., cell division, viral integration) matters, or (4) the qualitative behaviour of a deterministic model is noise-sensitive (e.g., bistability with a low barrier between states).

10.2 The Langevin Equation

The Langevin equation provides a bridge between deterministic and stochastic dynamics by adding a noise term to the ODE:

\[ \dot{x} = \beta - \gamma x + \sqrt{\beta + \gamma x}\,\xi(t). \]

The noise amplitude \(\sqrt{\beta + \gamma x}\) is the square root of the total propensity (a consequence of the Poisson statistics of individual reaction events in the linear noise approximation).

10.3 The Fokker-Planck Equation

Rather than tracking individual stochastic trajectories, one can describe the evolution of the probability distribution \(P(x, t)\) — the probability density of finding the system at state \(x\) at time \(t\).

\[ P_{stat}(x) \propto \exp\!\left(-\frac{2U(x)}{D}\right). \]

This is a Boltzmann-like distribution: the system is most likely to be found near the minima of the potential \(U(x)\) (stable fixed points), with probability exponentially decreasing toward the maxima (unstable fixed points). The ratio of the noise strength \(D\) to the potential barrier height determines the rate of noise-driven transitions between stable states.

10.4 Intrinsic vs. Extrinsic Noise in Gene Expression

10.5 The Two-State Promoter Model

A widely used stochastic model of gene expression is the two-state promoter (also called the telegraph model). The promoter switches between an active state (ON) and an inactive state (OFF):

\[ \text{OFF} \underset{k_{off}}{\stackrel{k_{on}}{\rightleftharpoons}} \text{ON} \xrightarrow{\beta} \text{mRNA} \xrightarrow{\gamma} \emptyset. \]

When the promoter is in the ON state, mRNA is produced at rate \(\beta\); in the OFF state, production is zero. mRNA degrades at rate \(\gamma\). This model accounts for bursty transcription: when \(k_{off} \gg k_{on}\) (the promoter spends most time OFF but occasionally fires), mRNA is produced in bursts, each burst lasting an exponentially distributed ON period of mean duration \(1/k_{off}\) and containing a geometrically distributed number of transcripts with mean \(\beta/k_{off}\).

---

Chapter 11: Cell Cycle Models

11.1 The Cell Cycle as a Biochemical Switch

The eukaryotic cell cycle — the ordered sequence of events by which a cell duplicates its genome (S phase) and divides (M phase), interspersed with gap phases (G1 and G2) — is controlled by a network of cyclin-dependent kinases (CDKs) and their regulatory partners. A central question for mathematical biology is: what ensures that the cell cycle is unidirectional and irreversible? Why does a cell committed to mitosis not spontaneously revert to interphase?

The answer lies in bistability. The transition from G1 to S phase (the “Start” transition in yeast) and from G2 to M phase (the “mitotic entry” transition) are both governed by bistable switches with hysteresis. Once the switch has been thrown in one direction, the cell remains committed even if the triggering signal is reduced — a property essential for robust cell cycle progression.

11.2 The Tyson–Novak Model: MPF Activation

The Tyson–Novak model (Tyson, 1991; Novak & Tyson, 1993) describes the mutual inhibition network controlling MPF (Maturation Promoting Factor), the complex of cyclin B and CDK1 that drives mitotic entry. The core circuit involves:

• Cyclin B (CycB) is synthesised constitutively and degraded in a CycB-dependent manner.
• CDK1 is always present but held inactive by phosphorylation on Tyr15 (by Wee1 kinase) and activated by dephosphorylation (by Cdc25 phosphatase).
• MPF = CycB/CDK1 complex. Active MPF activates Cdc25 (positive feedback) and inhibits Wee1 (double-negative = positive feedback).
• APC/C:Cdc20 (the Anaphase Promoting Complex) degrades cyclin B, providing negative feedback. Its activation requires MPF and additional signals.

The double-positive feedback loop (MPF activates Cdc25, which activates MPF; MPF inhibits Wee1, which normally inhibits MPF) creates a bistable switch.

11.3 Simplified Bistable MPF Model

A minimal model for the MPF switch uses a single variable \(m\) representing the fraction of CDK1 in the active (dephosphorylated) form, controlled by the balance of Cdc25 activation and Wee1 inactivation. The rate of CDK1 activation by Cdc25 is an increasing, ultrasensitive (sigmoidal) function of \(m\), while the rate of CDK1 inactivation by Wee1 is a decreasing sigmoidal function. The net equation:

\[ \frac{dm}{dt} = \frac{V_{Cdc25}\,m^n}{K_{Cdc25}^n + m^n} \cdot (1 - m) - \frac{V_{Wee1}}{1 + (m/K_{Wee1})^n} \cdot m, \]

where the first term is Cdc25-mediated activation of inactive CDK1 (fraction \(1-m\)) and the second term is Wee1-mediated inactivation of active CDK1 (fraction \(m\)).

---

Chapter 12: Spatial Models and Turing Instability

12.1 Reaction-Diffusion Systems

The models developed so far assume a well-mixed intracellular environment — all molecules are uniformly distributed in space and interact equally. This is a reasonable approximation for small cells (bacteria) or rapidly diffusing small molecules, but breaks down for larger cells (eukaryotes), tissue-level phenomena, and any situation where spatial gradients are the biological signal (morphogen gradients in development, calcium waves, neural pattern formation).

Reaction-diffusion (RD) systems couple chemical kinetics with molecular diffusion:

The remarkable discovery of Alan Turing (1952) was that diffusion — naively expected only to smooth out concentration differences — can actually drive pattern formation from a uniform state. A homogeneous steady state that is stable in the ODE (no diffusion) can become unstable when diffusion is included, if the two species diffuse at very different rates. This is diffusion-driven instability or Turing instability.

12.2 Turing’s Diffusion-Driven Instability

Suppose \((u^*, v^*)\) is a stable steady state of the well-mixed (ODE) system. To analyse the RD system, linearise about \((u^*, v^*)\) and write the perturbation as a Fourier mode \(\delta u = U e^{\sigma t + ikx}\), \(\delta v = V e^{\sigma t + ikx}\), where \(k\) is the spatial wavenumber. Substituting into the linearised PDE gives the dispersion relation — the growth rate \(\sigma(k)\) as a function of wavenumber \(k\).

\[ \sigma(k) = \frac{1}{2}\left[(J_{11} - D_u k^2) + (J_{22} - D_v k^2) \pm \sqrt{((J_{11} - D_u k^2) - (J_{22} - D_v k^2))^2 + 4J_{12}J_{21}}\right]. \]\[ h(k^2) := D_u D_v k^4 - (D_v J_{11} + D_u J_{22})k^2 + \det(J) < 0. \]\[ (D_v J_{11} + D_u J_{22})^2 > 4 D_u D_v \det(J). \]

---

Chapter 13: Calcium Oscillations

13.1 The Biology of Calcium Signalling

Calcium ions (Ca\(^{2+}\)) are ubiquitous second messengers: cytosolic Ca\(^{2+}\) concentration rises transiently in response to extracellular signals (hormones, neurotransmitters, growth factors) and triggers diverse cellular responses including muscle contraction, secretion, gene expression, and cell division. The resting cytosolic Ca\(^{2+}\) concentration is approximately \(100\,\text{nM}\), while the concentration in the endoplasmic reticulum (ER) lumen is approximately \(500\,\mu\text{M}\) — a 5000-fold gradient maintained by SERCA pumps. Stimulation causes a brief spike (to \(\sim 1\,\mu\text{M}\) cytosolic Ca\(^{2+}\)), sometimes followed by sustained oscillations with periods ranging from seconds to minutes.

The oscillation mechanism involves calcium-induced calcium release (CICR): Ca\(^{2+}\) released from the ER through IP3 receptors (IP3R) opens neighbouring IP3Rs, amplifying the release. This is a positive feedback loop. The subsequent negative feedbacks — Ca\(^{2+}\)-dependent inactivation of IP3Rs and reuptake by SERCA pumps — allow the system to reset, enabling oscillations.

13.2 The IP3 Receptor and CICR

The IP3 receptor is a ligand-gated Ca\(^{2+}\) channel on the ER membrane. It is activated by the second messenger inositol 1,4,5-trisphosphate (IP3, produced by phospholipase C in response to GPCR or receptor tyrosine kinase signalling) and by Ca\(^{2+}\) itself at low concentrations. At high Ca\(^{2+}\) concentrations, the receptor is inactivated — the bell-shaped dependence of IP3R open probability on cytosolic Ca\(^{2+}\) is the key biophysical feature enabling CICR-driven oscillations.

13.3 The Li–Rinzel Model

The Li–Rinzel model (Li and Rinzel, 1994) is a minimal two-variable reduction of the De Young–Keizer model of IP3R dynamics. The two variables are:

• \(c(t)\): cytosolic Ca\(^{2+}\) concentration (in \(\mu\)M)
• \(h(t)\): fraction of IP3 receptors not inactivated by Ca\(^{2+}\)

The equations are:

\[ \frac{dc}{dt} = J_{IP3R}(c, h) - J_{pump}(c) + J_{leak}, \]\[ \frac{dh}{dt} = \frac{h_\infty(c) - h}{\tau_h(c)}, \]\[ J_{IP3R} = v_1 m_\infty^3(p)\, n_\infty^3(c)\, h^3\, (c_{ER} - c), \]\[ J_{pump} = \frac{v_3 c^2}{k_3^2 + c^2}, \]\[ J_{leak} = v_2(c_{ER} - c), \]\[ m_\infty(p) = \frac{p}{p + d_1}, \quad n_\infty(c) = \frac{c}{c + d_5}, \quad h_\infty(c) = \frac{d_2}{c + d_2}, \]

and inactivation time constant \(\tau_h(c) = 1/(a_2(c + d_2))\). Here \(p\) is the IP3 concentration (treated as a fixed parameter), \(c_{ER}\) is the (approximately constant) ER Ca\(^{2+}\) concentration, and \(v_1, v_2, v_3, d_1, d_2, d_5, k_3, a_2\) are kinetic parameters.

---

Chapter 14: Metabolic Flux Analysis — Stoichiometry, Null Space, and Flux Modes

14.1 The Stoichiometric Framework Revisited

In Chapter 1 we introduced the stoichiometric matrix \(S\) as a compact representation of a reaction network. Here we extend this framework to metabolic flux analysis (MFA), a set of mathematical tools for understanding the distribution of reaction fluxes in a metabolic network at steady state. MFA is the quantitative foundation of metabolic engineering — the rational redesign of cellular metabolism for industrial biotechnology — and of the analysis of disease-associated metabolic reprogramming.

The steady-state constraint \(S\mathbf{v} = \mathbf{0}\) defines a linear subspace of flux space, the null space (or kernel) of \(S\). The dimension of this null space, \(\dim(\ker(S)) = n - \text{rank}(S)\), equals the number of degrees of freedom in the steady-state flux distribution — that is, the number of free flux variables that must be specified (by experiment or optimisation) to determine the complete flux distribution.

14.2 Null Space and Flux Cone

The set of all feasible steady-state flux vectors forms a flux cone: it is the intersection of the null space of \(S\) (the steady-state constraint) with the positive orthant (irreversibility constraints \(v_i \geq 0\) for irreversible reactions). The flux cone is a polyhedral cone, and its geometric structure encodes all possible metabolic behaviours.

The extreme rays of the flux cone — the generating vectors from which all feasible flux distributions can be composed as non-negative linear combinations — are the elementary flux modes (EFMs). Each EFM represents a minimal set of reactions that can carry a steady-state flux while satisfying stoichiometric and irreversibility constraints.

EFMs have a direct biological interpretation: each EFM corresponds to a distinct metabolic pathway through the network, carrying flux from specific inputs to specific outputs without violating mass balance at any internal metabolite.

14.3 Worked Example: A Four-Reaction Network

14.4 Connection to Flux Balance Analysis

Elementary flux modes provide the geometric foundation for flux balance analysis (FBA, Chapter 4). In FBA, one maximises a linear objective function \(c^T \mathbf{v}\) over the flux cone. By the theory of linear programming, the optimal solution always lies at a vertex of the feasible polytope (the intersection of the flux cone with capacity constraints \(v_i \leq v_i^{max}\)). Each vertex of this polytope is a non-negative linear combination of at most \(\rho\) EFMs — specifically, it is a basic feasible solution of the LP, corresponding to a minimal set of active reactions.

---

Appendix: Computational Tools

Python / Google Colab

Numerical integration of ODE systems uses scipy.integrate.solve_ivp (or the older odeint). Phase portraits are generated by plotting nullclines (solving \(f = 0\) and \(g = 0\) using numpy) and overlaying vector fields with matplotlib.pyplot.quiver. Stochastic simulations use custom Gillespie algorithm implementations or the gillespy2 library.

MATLAB

MATLAB’s ode45 (explicit Runge–Kutta 4/5 adaptive step) is the workhorse solver for non-stiff systems; ode15s is preferred for stiff systems (e.g., fast-slow systems or systems with widely varying time scales). The pplane package provides an interactive phase-plane graphical interface.

XPPAUT

XPPAUT (X-Windows Phase Plane plus AUTO) is the standard tool for bifurcation analysis in mathematical biology. AUTO, embedded in XPPAUT, performs numerical continuation: it follows a branch of steady states or limit cycles as a parameter varies, detecting bifurcation points automatically. XPPAUT .ode files specify the model symbolically, and the software handles numerical continuation, stability computation, and phase portrait generation interactively.

---

Appendix B: Summary of Key Formulas

Stoichiometric Matrix Framework

For a network with \(m\) species and \(n\) reactions, the ODE system is \(\dot{\mathbf{X}} = S\mathbf{v}(\mathbf{X})\) where \(S\) is the \(m \times n\) stoichiometric matrix. The key structural quantities are:

• Null space \(\ker(S)\): all flux vectors \(\mathbf{v}\) with \(S\mathbf{v} = \mathbf{0}\). At steady state, \(\mathbf{v}^*\) lies here.
• Left null space \(\ker(S^T)\): all conservation law vectors \(\mathbf{c}\) with \(\mathbf{c}^T S = \mathbf{0}\). Dimension = \(m - \text{rank}(S)\).
• Column space \(\text{im}(S)\): all achievable rates of change \(d\mathbf{X}/dt\).
• Row space: the space of stoichiometrically distinct combinations of species.

Enzyme Kinetics Summary

Model	Rate equation	Key parameters
Michaelis-Menten	\(v = V_{max}[S]/(K_M + [S])\)	\(K_M, V_{max}\)
Competitive inhibition	\(v = V_{max}[S]/(K_M(1+[I]/K_I) + [S])\)	\(K_I\) increases app. \(K_M\)
Non-competitive inhibition	\(v = V_{max}[S]/((K_M+[S])(1+[I]/K_I))\)	\(K_I\) decreases app. \(V_{max}\)
Hill equation	\(v = V_{max}[S]^n/(K_{half}^n + [S]^n)\)	\(n\): cooperativity; sigmoidal for \(n>1\)

Stability Classification at a 2D Fixed Point

Given Jacobian \(J\) with \(\tau = \text{tr}(J)\) and \(\Delta = \det(J)\):

\[ \begin{array}{ll} \Delta < 0: & \text{saddle (unstable)} \\ \Delta > 0,\; \tau < 0,\; \tau^2 > 4\Delta: & \text{stable node} \\ \Delta > 0,\; \tau < 0,\; \tau^2 < 4\Delta: & \text{stable spiral} \\ \Delta > 0,\; \tau > 0,\; \tau^2 > 4\Delta: & \text{unstable node} \\ \Delta > 0,\; \tau > 0,\; \tau^2 < 4\Delta: & \text{unstable spiral} \\ \Delta > 0,\; \tau = 0: & \text{centre (linear analysis; need nonlinear correction)} \end{array} \]

Metabolic Control Analysis

• Summation theorem: \(\sum_i C^J_{E_i} = 1\) (flux control coefficients sum to 1).
• Connectivity theorem: \(\sum_i C^J_{E_i}\,\varepsilon^{v_i}_S = 0\) for any intermediate metabolite \(S\).
• Response coefficient: \(R^J_p = C^J_E \cdot \varepsilon^v_p\) (how flux responds to parameter \(p\) via enzyme \(E\)).

Goldbeter-Koshland Function

\[ GK(u,v,J,K) = \frac{2uJ}{v-u+vJ+uK+\sqrt{(v-u+vJ+uK)^2 - 4(v-u)uJ}}. \]

When \(J, K \to 0\) (enzyme saturation), \(GK \to 0\) for \(u < v\) and \(GK \to 1\) for \(u > v\) — a perfect switch.