Agent-Based Modelling — Flocking, Epidemics, Ant Colonies and Emergent Traffic

A murmuration of starlings contains no choreographer. A traffic jam forms without any single driver deciding to stop. A virus spreads through a population following equations invented by epidemiologists a century ago. Agent-based modelling is the discipline that explains how local rules between simple agents produce complex global behaviour — and six interactive simulations make every concept visible in real time.

Bottom-Up vs Top-Down Modelling

Classical mathematical models describe systems from the top down: write down equations for aggregate quantities (total infected, average velocity, mean temperature), solve them, and read off the behaviour. This is powerful and efficient. But it struggles when individual heterogeneity matters, when spatial structure is important, or when agents have memory or adaptive strategies.

Agent-based modelling (ABM) takes the opposite approach: define rules for each individual agent, run thousands of agents simultaneously, and let the collective behaviour emerge. The global pattern is never explicitly programmed — it arises from local interactions. Real ABM power comes from the gap between these two approaches: many systems produce aggregate behaviour that is mathematically equivalent to top-down ODEs under certain conditions, but diverges from them whenever the assumptions (homogeneous mixing, infinite population, no spatial structure) break down.

Part 1: Flocking

Boids Simulation — three rules, infinite complexity

In 1986 Craig Reynolds published a paper that changed computer graphics forever. He showed that realistic flocking behaviour — the murmuration patterns of starlings, the tight balls of fish schools — emerges from just three local rules applied to each agent: separation, alignment, and cohesion. No global choreography, no central control, no explicit flock structure.

Boids Algorithm — Three Steering Rules

For each boid i, consider neighbours within radius R:

1. Separation (avoid crowding):
   f_sep = -Σⱼ (p_j - p_i) / |p_j - p_i|²   for |p_j-p_i| < r_sep
   Steers away from nearby boids

2. Alignment (steer toward average heading):
   f_ali = (Σⱼ v_j / |N|) - v_i              for j ∈ neighbours
   Steers toward average velocity of neighbourhood

3. Cohesion (steer toward centre of mass):
   f_coh = (Σⱼ p_j / |N|) - p_i - v_i        for j ∈ neighbours
   Steers toward average position of neighbourhood

Combined acceleration:
   a_i = w_sep·f_sep + w_ali·f_ali + w_coh·f_coh

Velocity update:
   v_i += a_i · dt;   |v_i| = clamp(|v_i|, v_min, v_max)
   p_i += v_i · dt

Spatial hash to find neighbours in O(1) per boid:
  Naive: O(n²) comparisons; Spatial hash: O(n·k), k=avg neighbours
  Critical: 500 boids → 125,000 pairs → 5×: 15,625,000 (~unusable)

From these three rules, the Boids simulation produces every hallmark of real flocking: tight formation flight, split and merge around obstacles, wave-like ripples through the flock, and the elongated "column" shape during cruise. Adjust the weight of each rule independently — turn separation up high and the flock explodes; turn cohesion to maximum and it collapses into a tight ball. The sweet spot produces the spooky, coordinated behaviour seen in nature.

Extension — Birds Flock: The Birds Flock simulation adds predator avoidance — a fourth rule that produces the rapid, coordinated escape manoeuvres filmed in real starling murmurations. When the predator approaches, the flock contracts, rotates, and expands in ways that confuse pursuit — a beautiful emergent defence strategy that no single bird planned.

Part 2: Epidemic Spreading

SIR Model — the mathematics of outbreaks

The SIR model, formulated by Kermack and McKendrick in 1927, is the foundation of mathematical epidemiology. Every pandemic response model — COVID-19 SIR extensions, influenza SEIR, dengue metapopulation models — descends from these three equations. Understanding SIR means understanding herd immunity thresholds, the shape of epidemic curves, and why vaccination at 70% coverage can prevent outbreaks even though 30% of people remain susceptible.

SIR Model — Equations & Key Parameters

State variables (fractions of population):
  S = Susceptible   dS/dt = -β·S·I
  I = Infected      dI/dt =  β·S·I - γ·I
  R = Recovered     dR/dt =  γ·I
  S + I + R = 1

Parameters:
  β = transmission rate  [contacts⁻¹·day⁻¹]
  γ = recovery rate      [day⁻¹]  ≡ 1/D, D = infectious period

Basic reproduction number:
  R₀ = β/γ = β·D
  R₀ > 1   → epidemic grows       (each case infects >1 others)
  R₀ < 1   → epidemic dies out
  R₀ = 1   → endemic equilibrium

Herd immunity threshold:
  p_c = 1 - 1/R₀   (vaccination fraction needed)
  Measles (R₀≈15):    p_c = 93%
  COVID-19 (R₀≈2.5):  p_c = 60%

Final size (fraction eventually infected):
  r_∞ satisfies: ln(1-r_∞)/r_∞ = -R₀

SEIR extension:
  Add E (exposed/latent): dE/dt = β·S·I - σ·E
  σ = rate from exposed to infectious (1/latent period)

The SIR simulation displays both the ODE solution (aggregate curves) and an agent-based version where individual people move through a grid and infect neighbours on contact. The two representations diverge when spatial clustering matters: in the agent-based version, dense clusters of susceptibles create localised outbreaks that persist after the aggregate SIR curve has flattened. Slide R₀ below 1 to see the epidemic extinguish itself before reaching the entire population.

Part 3: Stigmergy — Ant Colony Optimisation

Ants Simulation — pheromone communication

An ant colony has no central planner, yet it solves the shortest-path problem between nest and food source with remarkable efficiency. The mechanism is stigmergy: agents modify their shared environment (by depositing pheromone), and other agents respond to those environmental modifications. This indirect communication, mediated by the environment rather than direct contact, produces globally optimal solutions from locally simple rules.

Ant Colony Optimisation (ACO) — Pheromone Dynamics

Pheromone update rule:
  τᵢⱼ ← (1-ρ)·τᵢⱼ + Σᵏ Δτᵢⱼᵏ
  ρ     = evaporation rate ∈ (0,1)    (forgetting)
  Δτᵢⱼᵏ = Q/Lᵏ if ant k used edge (i,j), else 0
  Q     = pheromone strength constant
  Lᵏ    = total path length of ant k

Probabilistic path selection:
  P(edge j | at node i) = [τᵢⱼ]^α · [η᷊ⱼ]^β / Σₗ [τᵢₗ]^α · [ηᵢₗ]^β
  η᷊ⱼ = 1/dᵢⱼ (heuristic: prefer short edges)
  α   = pheromone importance weight
  β   = heuristic importance weight

Positive feedback (shortening loop):
  Shorter path → more trips per time → more pheromone deposited
  → higher probability chosen by future ants →
  → even shorter effective path chosen → ...

Convergence:
  Too high ρ: ants don't follow established paths (random)
  Too low ρ:  stale pheromone traps ants on suboptimal paths
  Optimum: ρ ≈ 0.1–0.3 for most problems

In the Ants simulation, watch pheromone trails form in real time as ants explore randomly, find food, and return to the nest. The first ant to find food takes a random, winding path. Subsequent ants that stumble across its pheromone trail are biased to follow it. Over time, a positive feedback loop selected a near-optimal route — one that has been experimentally confirmed to approximate the shortest path to within a few percent. Introduce a wall mid-trail and watch the colony route around it.

Part 4: Population Cycles — Predator-Prey Dynamics

Predator-Prey Simulation — Lotka-Volterra in two dimensions

The lynx-snowshoe hare population data collected by the Hudson's Bay Company between 1845 and 1935 shows regular oscillations with a period of about 10 years. Vito Volterra and Alfred Lotka, working independently in the 1920s, derived a simple pair of ODEs that produces exactly these oscillations — without knowing about the hare data. The Lotka-Volterra equations are one of the most famous models in ecology and remain the foundation of modern population dynamics.

Lotka-Volterra Equations — Predator-Prey Dynamics

Classic equations:
  dN/dt = αN - βNP       (prey: grow, die when eaten)
  dP/dt = δNP - γP       (predator: grow when fed, die naturally)

  N = prey population;  P = predator population
  α = prey growth rate
  β = predation rate (per predator-prey encounter)
  δ = predator efficiency (prey converted to predators)
  γ = predator death rate

Equilibrium:
  N* = γ/δ   (prey at equilibrium)
  P* = α/β   (predator at equilibrium)

Conserved quantity (first integral):
  V = δN - γ·ln(N) + βP - α·ln(P) = constant
  → Closed orbits in phase plane (neutral stability)

Lotka-Volterra paradox:
  Increasing prey growth rate α → raises predator eq, not prey eq
  "Paradox of enrichment": adding nutrients destabilises the cycle

Real-world extensions:
  Logistic growth: α → α(1-N/K)     (carrying capacity)
  Holling types II/III: saturating functional response
  Rosenzweig-MacArthur: leads to limit cycles, chaos

The Predator-Prey simulation runs the Lotka-Volterra equations alongside an agent-based spatial version where individual prey and predators move on a grid. The ODE solution shows perfect closed orbits in phase space; the agent-based version shows irregular oscillations with occasional extinctions — exactly the difference between idealised models and real ecological systems. Introduce spatial structure (patchy habitat, corridors) and watch the metapopulation behaviour that allows survivors to recolonise extinct patches.

Part 5: Urban Dynamics — Bus Bunching & City Growth

Bus Bunching — why buses travel in pairs

Bus bunching is one of the most elegant examples of emergent dysfunction in urban systems. Even when buses are scheduled at perfectly equal intervals, small random delays cause the trailing bus to pick up more passengers (since the gap grew), slowing it further and causing the leading bus to travel faster (since it has fewer stops). The result is inevitable: buses cluster into pairs. No driver is behaving badly — the bunching is a property of the system, not the agents.

Bus Bunching — Headway Dynamics

Headway (time gap between buses):
  h_i = time since previous bus at same stop

Boarding time depends on headway:
  T_board(i) = λ · h_i   (λ = passenger arrival rate per unit time)

Bus dwell time at stop s:
  d_s(i) = d_min + T_board(i)/μ   (μ = boarding rate per second)

Perturbation amplification:
  If bus i is delayed by δ:
    - h_i increases → more passengers → longer dwell at each stop
    - h_{i+1} decreases → fewer passengers → shorter dwell
  Next loop: h_i increases further → unstable runaway

Bunching criterion:
  System bunches if ∂h_i/∂(delay) > 1  (perturbations amplify)
  This holds whenever λ/μ is large enough

Control strategies:
  Headway-based holding: hold leading bus if h < h_min
  Optimal control (Daganzo 2009): minimise total passenger wait
  Real-time GPS dispatch: modern transit systems use this

The Bus Bunching simulation puts you in control of a circular bus line. Watch perfectly spaced buses inevitably converge into pairs over a few loops. Enable the headway-holding control strategy and see buses stabilise — but at the cost of occasional intentional delays that feel counter-intuitive to frustrated passengers.

City Growth Simulation — spatial economics

Cities grow according to a self-reinforcing logic: economic activity attracts workers, workers create demand, demand attracts more businesses, business density raises land value, and land value concentrates development. This feedback loop, formalised in urban economics, produces the power-law scaling of city size distributions (Zipf's law), the fractal structure of road networks, and the concentric ring structure of land use.

City Growth — Zipf's Law & Urban Scaling

Zipf's law (rank-size rule):
  Population of city ranked r: P(r) = P₁ / r^α
  Empirically: α ≈ 1 for most countries
  US (2020): NYC≈8.3M, LA≈3.9M, Chicago≈2.7M → rank×pop ≈ const

Kleiber-scaling of urban metabolism:
  Y ∝ N^β     (Y = economic output, N = population)
  β > 1:      superlinear → cities become disproportionately productive
  β ≈ 1.15    (patents, GDP, wages)
  β ≈ 0.85    (infrastructure: roads, pipes)

Alonso-Muth-Mills monocentric model:
  Rent gradient: R(x) = R(0) - t·x  (R falls with distance x from CBD)
  t = commuting cost per unit distance
  Result: density peaks at centre, decays exponentially outward

Cellular automata land use (SLEUTH model):
  Urban spread probability:
    P = f(slope, exclusions, urbanised neighbours, distance to roads)
  Calibrated to real cities using historical satellite imagery

The City Growth simulation implements a cellular-automata land-use model where cells transition from rural to suburban to urban based on neighbourhood density, road access, and slope. Place an initial economic seed (a harbour, a rail junction) and watch the city nucleate, sprawl along transport corridors, and self-organise into the characteristic ring structure — all from local rules without any global plan.

Agent-Based Modelling Collection

Cross-Collection Connections

Agent-based thinking connects the entire platform. The Boids algorithm is the direct ancestor of the Crowd Simulation — the same three steering rules applied to pedestrians navigating corridors and exits. The SIR epidemic model shares its mathematical form with the Predator-Prey ODEs — both are two-variable, nonlinear systems with limit cycle behaviour. Ant colony pheromones are a biological implementation of the gradient descent that powers neural network training: both systems follow concentration gradients to find optima. City growth emergence connects to the Percolation simulation — urban spreading is a directed percolation process across the land-use lattice, and the city boundary in the growth simulation is precisely the percolation cluster frontier. All of these systems illustrate the same deep principle: local rules, iterated over many agents and many time steps, reliably produce global order that was never explicitly designed.

Algorithms & Methods in This Collection

Reynolds steering forces Spatial hash neighbour query Velocity clamping Euler ODE integration Runge-Kutta 4 SIR compartmental model SEIR extension Pheromone evaporation map ACO roulette selection Lotka-Volterra phase portrait Agent movement on grid Statistical extinction detection Headway-based bus holding CA land-use transition rules Zipf distribution fitting Monocentric rent gradient