Why ODEs Are Biology's Language
A differential equation asks: how fast does this quantity change? Biology is full of rates — rates of cell division, rates of predation, rates of ion channel opening, rates of enzyme catalysis. Whenever you have a rate that depends on the current state of the system, you have a differential equation. When multiple rates depend on each other, you have a system of coupled ODEs — and that's where interesting dynamics emerge.
Unlike physics, biological models are rarely derived from first principles alone. They blend first-principles constraints (conservation of mass, thermodynamics) with empirical rate laws (Hill kinetics, Michaelis-Menten, logistic growth) that are justified by their fit to data. The result is a pragmatic mathematical biology: models that are usefully predictive without being exact.
Part 1: Population Dynamics
Exponential Growth and the Logistic Equation
The simplest population model starts with exponential growth: dN/dt = rN. Every individual reproduces at rate r, so the growth rate is proportional to the current size. This gives N(t) = N₀ e^(rt) — unbounded growth. Real populations are bounded by resources, so Verhulst (1838) added a saturation term: dN/dt = rN (1 − N/K), where K is the carrying capacity. The logistic equation has an S-shaped (sigmoidal) solution that levels off at K.
Logistic Growth
dN/dt = rN(1 − N/K)
r = intrinsic growth rate (yr⁻¹)
K = carrying capacity (individuals)
Solution: N(t) = K / (1 + ((K − N₀)/N₀) · e^(−rt))
Fixed points: N* = 0 (unstable), N* = K (stable)
Maximum growth rate at N = K/2 (inflection point)
Doubling time at low density ≈ ln(2) / r
Lotka-Volterra: Two-Species Predator–Prey
Add a second species — a predator — and the equations become a 2D dynamical system. The Lotka-Volterra model consists of two coupled nonlinear ODEs. The prey grows exponentially in the absence of predators; predators die exponentially in the absence of prey. Predation terms couple them: βNP reduces prey, δNP increases predators.
Lotka-Volterra System
dN/dt = αN − βNP (prey)
dP/dt = δNP − γP (predator)
Nullclines: N* = γ/δ (vertical), P* = α/β (horizontal)
The fixed point (N*, P*) is a centre — neutral stability
Trajectories in phase space are closed curves (conserved quantity H)
H = δN − γ ln N + βP − α ln P = const.
Real systems: damped spirals toward equilibrium (carrying capacity + noise)
Prey–Predator Simulation
Agent-based Lotka-Volterra: rabbits graze and reproduce, foxes hunt and starve. Watch the classic counterclockwise phase-space cycle. Adjust birth and death rates to explore different dynamical regimes — stable cycles, extinction, coexistence.
Food Web (6-species)
Extended ODE system with six coupled species: Grass, Shrubs, Rabbit, Deer, Fox, Wolf. RK4 integration. Network visualisation where node size tracks population. Trigger trophic cascades by removing species and observe the system-wide response.
Phase plane analysis is the key tool. Instead of plotting N(t) and P(t) separately, plot P against N. In this phase plane, the trajectory reveals the attractor structure: is the equilibrium a stable spiral (populations converge), a centre (neutral oscillation), or a saddle (one population collapses)? The shape of the trajectory tells you everything about long-term behaviour without solving the ODE analytically.
Part 2: Neuroscience — The Hodgkin-Huxley Model
In 1952, Alan Hodgkin and Andrew Huxley published the first quantitative model of the action potential — the electrical impulse that neurons use to communicate. They inserted microelectrodes into the giant axon of a squid, voltage-clamped the membrane at different potentials, and measured the currents. From these experiments they extracted conductance curves for sodium and potassium channels, then wrote four coupled ODEs that reproduced the action potential with remarkable fidelity. They received the Nobel Prize in Physiology or Medicine in 1963.
The Hodgkin-Huxley model is a masterpiece of mechanistic biological modelling. The membrane is treated as an electrical circuit: membrane capacitance C_m stores charge, and conductances g_Na, g_K, and g_L (leak) carry current. The voltage-dependent conductances are described by gating variables — m, h (sodium) and n (potassium) — each obeying their own first-order ODE.
Hodgkin-Huxley System (4 coupled ODEs)
C_m · dV/dt = I_ext − g_Na·m³h·(V − E_Na) − g_K·n⁴·(V − E_K) − g_L·(V − E_L)
dm/dt = α_m(V)·(1−m) − β_m(V)·m
dh/dt = α_h(V)·(1−h) − β_h(V)·h
dn/dt = α_n(V)·(1−n) − β_n(V)·n
V = membrane voltage (mV), m = Na activation, h = Na inactivation
n = K activation; all gating vars ∈ [0, 1]
g_Na = 120, g_K = 36, g_L = 0.3 mS/cm² (Hodgkin-Huxley 1952 values)
E_Na = +50 mV, E_K = −77 mV, E_L = −54.4 mV (reversal potentials)
Hodgkin-Huxley Neuron
Full conductance-based neuron: V(t) with Na⁺, K⁺, leak currents. Gating variables m, h, n plotted in real time. Presets: single spike, repetitive firing, sub-threshold, high-frequency burst. Euler integration at dt = 0.01 ms. Live spike rate and threshold detection.
Cardiac Action Potential
Ventricular action potential: five phases (resting, upstroke, early rapid repolarisation, plateau, final repolarisation). Na⁺/Ca²⁺/K⁺ channel conductances, refractory period, pacemaker automaticity. Compare ventricular, atrial, and SA node waveforms.
The power of the Hodgkin-Huxley framework is its generality. The same formal structure — capacitive membrane, voltage-gated conductances, first-order gating kinetics — describes neurons across phyla, cardiac cells, pancreatic beta cells, and skeletal muscle fibres. Modern computational neuroscience uses multi-compartment extensions with dozens of ion channel types, but the mathematical skeleton is always Hodgkin-Huxley.
Part 3: Enzyme Kinetics — Michaelis-Menten
Enzymes are biological catalysts. They bind a substrate, convert it to a product, and release the product — all without being consumed themselves. The rate of an enzyme-catalysed reaction was described mathematically by Leonor Michaelis and Maud Menten in 1913. Their model involves a two-step mechanism: binding (reversible) and catalysis (irreversible), with rapid equilibrium between the enzyme–substrate complex and free enzyme.
Michaelis-Menten Kinetics
Reaction: E + S ⇌ ES → E + P
Reaction rate: v = V_max · [S] / (K_m + [S])
V_max = k_cat · [E_total] (maximal velocity)
K_m = (k₋₁ + k_cat) / k₁ (Michaelis constant)
At [S] = K_m: v = V_max / 2 (half-maximal rate)
At [S] ≪ K_m: v ≈ (V_max / K_m) · [S] (first-order in S)
At [S] ≫ K_m: v ≈ V_max (zeroth-order, saturated)
Lineweaver-Burk: 1/v = (K_m / V_max) · (1/[S]) + 1/V_max
Enzyme Kinetics
Michaelis-Menten rate curve v vs [S], Lineweaver-Burk double reciprocal plot, competitive / noncompetitive / uncompetitive inhibition, substrate depletion curve. V_max and K_m sliders. See how inhibitors shift the apparent K_m and V_max differently.
Chemical Reaction Kinetics
Arrhenius rate constant k = A·e^(−Ea/RT), sequential reactions A→B→C, catalyst effect on activation energy. First and second order kinetics. Watch concentration vs time curves change as you tune activation energy and temperature.
Hill kinetics generalise Michaelis-Menten to cooperative systems. Many biological responses are not hyperbolic but sigmoidal — the Hill equation v = V_max · [S]ⁿ / (K_d + [S]ⁿ) where n > 1 describes cooperative binding. Haemoglobin binds oxygen cooperatively (n ≈ 2.8): once one O₂ is bound, subsequent binding is easier. This creates a steep switch-like response that makes haemoglobin an efficient oxygen carrier. The same cooperative logic appears in transcription factor binding, ion channel gating, and circadian clock repression.
Part 4: Epidemiology — The SIR Model
The SIR model divides a population into three compartments: Susceptible (can catch the disease), Infected (currently infectious), and Recovered (immune). Kermack and McKendrick published this framework in 1927, motivated by the 1918 influenza pandemic. Despite its simplicity, the model captures the core dynamics of infectious disease outbreaks — including the critical concept of herd immunity.
SIR Compartmental Model
dS/dt = −β · S · I / N
dI/dt = β · S · I / N − γ · I
dR/dt = γ · I
N = S + I + R = const. (closed population, no births/deaths)
β = transmission rate (contacts × probability per contact)
γ = recovery rate (= 1 / infectious period)
R₀ = β / γ (basic reproduction number)
Epidemic threshold: R₀ > 1 (outbreak occurs)
Herd immunity threshold: fraction immune p_c = 1 − 1/R₀
Final size equation: S_∞ = N · e^(−R₀(1 − S_∞/N))
SIR Epidemic Model
Live S/I/R curves and 2D agent grid. Sliders for β (transmission) and γ (recovery). Watch epidemic threshold: below R₀ = 1, the infection dies out; above it, an epidemic wave sweeps through. Enable vaccination to see herd immunity in action.
Circadian Rhythm Oscillator
Goodwin oscillator producing 24-hour melatonin, cortisol, and temperature cycles. Classic example of a biological ODE feedback loop. Demonstrates sustained oscillation from a negative feedback loop with time delay — the principle powering all circadian clocks.
Part 5: ODE Solvers — Getting the Numerics Right
Biological ODEs are rarely solvable analytically. We solve them numerically: starting from initial conditions, we step forward in time using an approximation rule. The simplest is Euler's method: x(t + dt) = x(t) + dt · f(x, t). It is easy to implement but accumulates error rapidly and can be unstable for stiff systems.
The 4th-order Runge-Kutta (RK4) method is the workhorse of biological simulation. It evaluates the right-hand side four times per step and combines the results with optimal weights — giving 4th-order accuracy at moderate cost. For stiff systems (dramatically different timescales, as in enzyme cascades), implicit methods such as VODE or numerical differentiation formulas (NDFs) are required to avoid numerical blowup with reasonable step sizes.
RK4 Integration Step
k₁ = f(xₙ, tₙ)
k₂ = f(xₙ + dt/2 · k₁, tₙ + dt/2)
k₃ = f(xₙ + dt/2 · k₂, tₙ + dt/2)
k₄ = f(xₙ + dt · k₃, tₙ + dt)
x_{n+1} = xₙ + (dt/6)·(k₁ + 2k₂ + 2k₃ + k₄)
Local truncation error: O(dt⁵); global error: O(dt⁴)
Euler method error: O(dt²) per step, O(dt) global
Stiff systems: use dt ≪ 1/λ_max (largest eigenvalue)
Stiffness is the hidden enemy of biological simulations. The Hodgkin-Huxley sodium channel activates in ~0.2 ms, but a neuron simulation might run for 1 000 ms. The ratio of timescales is ~5 000. With Euler's method, dt must be smaller than the fastest timescale (0.01 ms), requiring 100 000 steps. RK4 at dt = 0.01 ms is stable and accurate. But for enzyme cascades with timescales spanning milliseconds to hours, adaptive implicit solvers become essential — otherwise simulations are either unstable or prohibitively slow.
Putting It All Together: Reading a Phase Diagram
The phase plane is the fundamental visualisation tool for 2D ODE systems. Plot one variable against another; the trajectory of the system traces a curve. Equilibria appear as fixed points; their stability is determined by the eigenvalues of the Jacobian matrix at the fixed point. Negative real parts = stable (spiral or node); positive real parts = unstable; zero real parts = centre (neutral).
With more than two variables, we cannot directly visualise the full phase space; instead, we use Poincaré sections (a plane through the trajectory), bifurcation diagrams (equilibrium vs parameter), and numerical continuation. These tools reveal how qualitative behaviour changes as a parameter crosses a bifurcation point: a stable equilibrium becomes unstable and gives birth to a limit cycle (Hopf bifurcation), or two equilibria collide and annihilate each other (saddle-node bifurcation). The circadian clock, the heartbeat, and the neural action potential are all best understood through this lens.