Ecology & Population Dynamics — Predator-Prey Cycles, Food Webs, Trophic Cascades and the Carbon Cycle

The Mathematics of Living Systems

Ecological models occupy an unusual position in science: they are simple enough to write on one line yet generate behaviour — oscillations, chaos, extinction, coexistence — that surprises ecologists every year. The Lotka-Volterra predator-prey equations date from 1925 and 1926. A century later, they remain the foundation of population ecology, extended and refined but never replaced.

What the simulations in this collection share is a common mathematical language: coupled nonlinear ordinary differential equations (ODEs) in population variables. Wolves and elk, plankton and krill, trees and atmospheric CO₂ — all are described by the same form of equation, differing only in the number of interacting variables and the structure of the interaction terms.

Simulation 1: Predator-Prey Dynamics

Lotka-Volterra — oscillations without a clock

Alfred Lotka (1925) and Vito Volterra (1926) derived the same pair of equations independently. Volterra was trying to explain why the proportion of predatory fish in the Adriatic Sea rose during the First World War — when fishing was suspended, both predator and prey increased, but not proportionally. The equations explained why: the two species are locked in a cycle governed by a conserved quantity analogous to energy in classical mechanics. Stop the external perturbation (fishing) and the cycle resumes.

State variables:
  N = prey population (rabbits, fish, deer, ...)
  P = predator population (foxes, sharks, wolves, ...)

ODEs:
  dN/dt = α·N  - β·N·P     (prey: grow, get eaten)
  dP/dt = δ·N·P - γ·P      (predator: gain from prey, die)

Parameters:
  α = prey birth rate         [day⁻¹]
  β = predation rate          [day⁻¹ per predator]
  γ = predator death rate     [day⁻¹]
  δ = efficiency (prey→pred.) [dimensionless conversion]

Fixed points:
  Trivial:    (N, P) = (0, 0)   — extinction (unstable)
  Coexistence: (N*, P*) = (γ/δ, α/β)   — neutral centre

Conserved quantity (Lyapunov function):
  V(N,P) = δN - γ·ln N + βP - α·ln P = const along trajectories
  → closed orbits in phase plane (neutrally stable cycles)

Summary: The Lotka-Volterra system is "structurally unstable" —
any small modification (prey crowding, saturating predation)
changes neutral cycles to spirals in or out. Real systems need
the Holling-type II or III functional responses for stability.

The Lotka-Volterra simulation displays both the time series (population vs. time) and the phase portrait (prey population on the x-axis, predator on the y-axis) simultaneously. You can see the closed orbits in phase space — the mathematical signature of conservative oscillations. Slide the initial conditions around the fixed point and watch the amplitude change while the period barely shifts: a hallmark of the LV equations that any modification with density‐dependence in either species would destroy.

Simulation 2: Food Web Dynamics

Food Web — multi-species interaction networks

Real ecosystems contain not two species but hundreds, linked by a complex web of who-eats-whom relationships. A food web is a directed graph where each node is a species (or trophic guild) and each edge represents energy flow. The web's structure — its connectance, trophic levels, and the distribution of interaction strengths — determines whether the ecosystem is stable or fragile.

For n species with population vector x = (x₁, x₂, ..., xₙ):

  dxᵢ/dt = xᵢ (rᵢ + Σⱼ aᵢⱼ·xⱼ)

  rᵢ  = intrinsic growth rate (+ for producers, - for top predators)
  aᵢⱼ = interaction coefficient matrix:
         aᵢⱼ > 0: species j benefits species i (prey eaten by i)
         aᵢⱼ < 0: species j harms species i (competitor, predator)
         aᵢⱼ = 0: no direct interaction

Stability criterion (May 1972):
  A random ecosystem with n species, mean interaction strength s,
  and connectance C is stable if:
         s · √(n·C) < 1
  (stability decreases as n, C, or s increases)

This famous result suggested that diverse, highly connected
ecosystems should be unstable — contradicting ecologists' field
intuition. The resolution: real food webs are NOT random;
they have specific structures (weak links, omnivory, loops)
that confer stability.

Trophic levels:
  L1: Primary producers     (plants, phytoplankton) — rᵢ > 0
  L2: Primary consumers     (herbivores)
  L3: Secondary consumers   (predators of herbivores)
  L4: Apex predators        — rᵢ < 0, depend entirely on prey
  Detritivores: break down dead matter → recycle to L1

The Food Web simulation builds a multi-species ecosystem with up to 8 species across 3–4 trophic levels. Remove a species with a click and watch the cascade: removing an apex predator may cause a "mesopredator release" where mid-level predators explode in population. Add a nutrient subsidy to a producer and watch enrichment destabilise the web. The simulation computes the community matrix eigenvalues in real time and displays a stability indicator.

Simulation 3: Trophic Cascades

Trophic Cascade — top-down control and keystone species

A trophic cascade occurs when a predator indirectly benefits plants by suppressing herbivores. The Yellowstone wolf reintroduction is the most famous terrestrial example; sea otters, wolves of kelp forests, keeping sea urchins in check, are the marine equivalent. Understanding trophic cascades is central to ecological restoration and conservation biology.

Three species: Plant (V), Herbivore (H), Carnivore (C)

dV/dt = r·V·(1 - V/K)  - f(V)·H
dH/dt = e₁·f(V)·H      - g(H)·C  - d_H·H
dC/dt = e₂·g(H)·C      - d_C·C

Functional responses:
  f(V) = a·V / (1 + a·h·V)   — Holling Type II (saturating predation)
  g(H) = b·H / (1 + b·k·H)   — Holling Type II for carnivores
  a, b = attack rates
  h, k = handling times (max ingestion rate = 1/h)

Cascade mechanism:
  Increase C → decrease H → increase V   (positive indirect effect)
  This is a negative-positive-negative chain → net positive (cascade)

Keystone species criterion:
  A species is a "keystone" if its removal has disproportionately
  large effects relative to its biomass (Paine 1969).
  Keystones are typically apex predators with low biomass
  but high per-capita interaction strength.

Behavioural suppression (the "landscape of fear"):
  Predators need not eat prey to control them.
  Fear of predation restricts herbivore foraging locations,
  reducing grazing impact in predator-frequented areas.
  This "landscape of fear" doubles the realized top-down effect.

Simulation 4: The Ocean's Twilight Zone

Aphotic Zone — life without light

Below 200 metres the last photon of sunlight is absorbed. No photosynthesis is possible. Yet the deep ocean is not lifeless: the "biological pump" — the sinking of organic matter from the sunlit photic zone — delivers a flux of energy that supports one of Earth's largest and least-explored ecosystems. Marine snow, fecal pellets, dead organisms, and migrating zooplankton carry carbon from surface to abyss, sequestering CO₂ on timescales of centuries.

Carbon export flux (Martin 1987 power law):
  F(z) = F₀ · (z / z₀)^(-b)
  F₀ = export flux at reference depth z₀ (typically 100 m)
  b  ≈ 0.86 (global average, varies 0.4–1.4)
  z  = depth [m]

Example:
  F(100m)  = 10 g C/m²/yr  → surface export
  F(500m)  = 10 × (5)^-0.86 ≈ 2.3 g C/m²/yr
  F(4000m) = 10 × (40)^-0.86 ≈ 0.5 g C/m²/yr  (only 5% reaches bottom)

Organisms at depth:
  Mesopelagic zone (200–1000 m):
    Myctophids, squid, copepods — diel vertical migration
    Travel 200-300m nightly to feed at surface, return at dawn
    This migration transports ~1–5 Gt C/year globally

  Bathypelagic zone (1000–4000 m):
    Anglerfish, viperfish — chemoreception, bioluminescence
    No sunlight → no photosynthesis → all energy from marine snow

Diel vertical migration energy cost:
  A copepod migrating 200 m/night expends ~15% of its daily energy budget
  in swimming, but gains ~30% more food access — net positive

Simulation 5: The Global Carbon Cycle

Carbon Cycle — Earth's most important ecological flow

Carbon cycles through five major reservoirs: atmosphere, terrestrial biosphere (plants and soils), ocean surface layer, deep ocean, and lithosphere (fossil fuels and carbonate rocks). Natural fluxes between these reservoirs were approximately balanced before the industrial era. Since 1750, the burning of fossil fuels and land-use change have added ~640 Gt of carbon to the atmosphere — about half has been absorbed by the ocean and land biosphere; the rest has accumulated, raising atmospheric CO₂ from 280 to over 420 ppm.

Reservoirs (Gt C, approximate 2024 values):
  Atmosphere:          850 Gt C     (CO₂: ~420 ppm)
  Land plants:        550 Gt C
  Soil + detritus:   2500 Gt C
  Surface ocean:       900 Gt C
  Deep ocean:       38,000 Gt C
  Fossil fuels:    ~1200 Gt C remaining recoverable

Annual fluxes (natural, approximate):
  GPP (gross primary prod.):  ~120 Gt C/yr   (photosynthesis)
  Ecosystem respiration:      ~118 Gt C/yr   (back to atm.)
  Ocean uptake (air-sea):      ~2.3 Gt C/yr  (net into ocean)
  Weathering into ocean:        ~0.4 Gt C/yr

Human perturbation (2023):
  Fossil fuel + industry:     +9.7 Gt C/yr
  Land use change:            +1.1 Gt C/yr
  Total emissions:           +10.8 Gt C/yr
  Atmospheric accumulation:   +5.1 Gt C/yr  (residual)
  Ocean uptake:               +2.9 Gt C/yr
  Land uptake:                +2.8 Gt C/yr

CO₂ airborne fraction ≈ 47–50% (rest absorbed by ocean + land)
Ocean acidification:  pH ≈ 8.2 (1850) → 8.1 (now) → ~7.9 (2100 BAU)

Simulation 6: Flocking and Collective Behaviour

Birds Flock — emergent structure from local rules

A murmuration of starlings — thousands of birds moving in fluid, coordinated waves — has no conductor. Each bird follows three local rules (separation, alignment, cohesion) that together produce coordinated group motion robust to perturbation. From an ecological perspective, flocking behaviour evolved primarily for predator avoidance: a flock of 10,000 birds is 10,000 times harder to ambush than a single bird, and the confusion effect makes individual targeting almost impossible.

Standard Boids (Reynolds 1987):
  f_sep = -Σⱼ (pⱼ - pᵢ) / |pⱼ - pᵢ|²   [separation]
  f_ali = (Σⱼ vⱼ / |N|) - vᵢ            [alignment]
  f_coh = (Σⱼ pⱼ / |N|) - pᵢ - vᵢ      [cohesion]

Predator avoidance (4th rule):
  f_pred = -k_p · (pᵢ - p_pred) / |pᵢ - p_pred|³   if |pᵢ - p_pred| < R_alert
  Inverse-cube repulsion produces sharp, sudden escape manoeuvres

Collective anti-predator behaviours:
  Vacuole: flock splits, leaves empty space around predator
  Flash expansion: all birds simultaneously flee outward
  Coordinated turn: wave of turning propagates at ~15 ms/bird (starling data)
  Confusion effect: predator lock-on time ∝ n (flock size)

Information propagation:
  A "startle" response in one bird propagates to neighbours in ~13 ms.
  With 10 nearest-neighbour interactions, a flock of 1000 updates in
  ~125 ms — fast enough to appear instantaneous to a predator.

Hamilton's selfish herd (1971):
  Each animal minimises its "domain of danger" (area closer to it
  than to any other animal). This drives animals to move toward
  the centre of the group — producing aggregation without cohesion-force.

Simulation 7: Crop Growth and Agronomy

Crop Growth — from seed to harvest via growing degree days

Agriculture is applied ecology: managing the population dynamics of crop plants to maximise food production. Crop growth models — from simple degree-day counters to full process-based models like APSIM and DSSAT — predict phenological stages (germination, flowering, maturity), biomass accumulation, and yield from weather inputs. They are used everywhere: optimising planting dates, estimating yield under climate scenarios, and scheduling irrigation.

Growing Degree Days (GDD) accumulated:
  GDD_day = max(T_base, (T_max + T_min)/2) - T_base
  GDD_cumulative += GDD_day if > 0

Phenology thresholds (wheat, approximate):
  Germination:   ~150 GDD₅  (base 5°C)
  Tillering:     ~450 GDD₅
  Heading:      ~1300 GDD₅
  Physiological maturity: ~2000 GDD₅

Biomass accumulation (Monteith 1977):
  ΔB = ε_c · PAR_intercepted · FPAR
  ε_c = radiation use efficiency [g/MJ]  (wheat ≈ 1.2 g/MJ)
  FPAR = fraction of PAR intercepted by canopy = 1 - e^(-k·LAI)
  k    = extinction coefficient ≈ 0.5
  LAI  = leaf area index [m² leaf / m² ground]

Yield:
  Y = B_total · HI
  HI = harvest index (grain / total biomass) ≈ 0.45 wheat, 0.50 maize

Water stress factor:
  ks = min(1, ETa / ETc)
  Multiplies ΔB: limited water → limited growth
  ETa = actual evapotranspiration, ETc = potential

The Structure of All These Models

Looking across all seven simulations a pattern emerges. Every model, from the two-species Lotka-Volterra to the five-reservoir carbon cycle, has the same form: a vector of state variables x and a right- hand side F(x, parameters):

dx/dt = F(x, p)

The models differ in what goes in x (species populations, carbon reservoir sizes, crop biomass), what structure F has (bilinear product terms for predation, saturating Holling-II terms for predator satiation, logistic terms for density-dependent competition), and what the parameters mean. But the mathematical task — integrate these ODEs forward in time, find fixed points, assess stability via the Jacobian eigenvalues, explore bifurcations as parameters vary — is always the same.

The full catalogue of ecology simulations — browse all ecology simulations — also includes models of the global ocean's aphotic zone, the soil-erosion RUSLE model, and the plant-growth GDD model. The Spotlight #17 on Ecology & Biology covers earlier simulations in the collection. For the underlying mathematics, the Learning #16 on ODEs in Biology walks through each equation type with derivations.