Lotka–Volterra Equations — Predator–Prey Population Cycles

1. The two-species model and assumptions

The classic Lotka–Volterra system involves two populations interacting in a predator–prey relationship:

The model rests on four key assumptions:

• Unlimited prey growth: in the absence of predators, prey grow exponentially at rate α.
• Predator-only mortality: in the absence of prey, predators die exponentially at rate δ.
• Mass-action encounters: predation events are proportional to the product x·y — every predator is equally likely to encounter every prey individual.
• Closed ecosystem: no immigration, emigration, or seasonal effects. Only prey and predator.

2. System of ODEs — derivation

Let x(t) be the prey population and y(t) be the predator population. The four parameters are:

The governing equations are:

The mass-action term β·x·y is the key nonlinearity. Without it the equations decouple into two independent exponentials. With it, the populations are coupled and oscillate indefinitely.

3. Equilibria and stability analysis

Setting dx/dt = dy/dt = 0 reveals two fixed points:

Linearisation at E*: The Jacobian matrix evaluated at the non-trivial equilibrium has purely imaginary eigenvalues:

Purely imaginary eigenvalues indicate a centre in the linearised system. For the fully nonlinear Lotka–Volterra equations, this centre is neutrally stable: perturbations do not grow or decay — they orbit the equilibrium forever on a closed trajectory whose shape depends on initial conditions.

4. Phase portrait and conservation law

The Lotka–Volterra system has a remarkable first integral (conserved quantity), analogous to energy conservation in mechanics:

Each closed orbit in the phase plane (x, y) corresponds to a constant value of V. The orbit closest to the equilibrium E* has the smallest V; the amplitude of oscillation grows as V increases.

In the phase portrait you will observe:

• Clockwise closed curves around E* = (γ/δ, α/β).
• Prey peaks before predator peaks — the lag is a quarter-period (≈ T/4).
• Outer curves reach dangerously low minima — near-extinction before recovery.
• The equilibrium itself is the innermost point (infinitesimal oscillation).

5. Parameter effects — what changes the cycle?

Each parameter changes a specific aspect of the oscillation:

Adding a logistic cap for prey (αx(1−x/K)) transforms the neutral centre into a stable spiral: trajectories spiral inward, reaching a unique limit cycle (if the equilibrium lies inside the feasible region) or a stable point (if K is small).

6. JavaScript implementation with RK4

The 4th-order Runge–Kutta method (RK4) gives excellent accuracy at reasonable cost. Each step requires four evaluations of the right-hand side function:

// Lotka–Volterra right-hand side
function lv(x, y, α, β, γ, δ) {
  return {
    dxdt:  α * x - β * x * y,
    dydt: -γ * y + δ * x * y,
  };
}

// RK4 single step
function rk4Step(x, y, h, α, β, γ, δ) {
  const k1 = lv(x,            y,            α, β, γ, δ);
  const k2 = lv(x + h/2*k1.dxdt, y + h/2*k1.dydt, α, β, γ, δ);
  const k3 = lv(x + h/2*k2.dxdt, y + h/2*k2.dydt, α, β, γ, δ);
  const k4 = lv(x + h*k3.dxdt,   y + h*k3.dydt,   α, β, γ, δ);

  return {
    x: x + (h/6) * (k1.dxdt + 2*k2.dxdt + 2*k3.dxdt + k4.dxdt),
    y: y + (h/6) * (k1.dydt + 2*k2.dydt + 2*k3.dydt + k4.dydt),
  };
}

// Simulate for T years with step h
function simulate(x0, y0, α, β, γ, δ, T, h = 0.01) {
  const history = [{ t: 0, x: x0, y: y0 }];
  let { x, y } = { x: x0, y: y0 };

  for (let t = h; t <= T; t += h) {
    ({ x, y } = rk4Step(x, y, h, α, β, γ, δ));
    // Guard against numerical blow-up near extinction
    if (x < 0) x = 0;
    if (y < 0) y = 0;
    history.push({ t, x, y });
  }
  return history;
}

// Example: classic parameter set
const α = 1.0;  // prey birth rate
const β = 0.1;  // predation rate
const γ = 0.4;  // predator death rate
const δ = 0.02; // predator conversion efficiency

const data = simulate(10, 5, α, β, γ, δ, 50);
// data[i].x → prey population at time data[i].t
// data[i].y → predator population at time data[i].t

Rendering the phase portrait

// Draw (x, y) trajectory on a 2D canvas
function drawPhasePortrait(ctx, data, W, H) {
  const maxX = Math.max(...data.map(p => p.x)) * 1.1;
  const maxY = Math.max(...data.map(p => p.y)) * 1.1;

  const cx = p => (p.x / maxX) * W;
  const cy = p => H - (p.y / maxY) * H;

  ctx.beginPath();
  data.forEach((p, i) => {
    i === 0
      ? ctx.moveTo(cx(p), cy(p))
      : ctx.lineTo(cx(p), cy(p));
  });
  ctx.strokeStyle = '#22c55e';
  ctx.lineWidth = 1.5;
  ctx.stroke();

  // Mark equilibrium E* = (γ/δ, α/β)
  const eq = { x: γ/δ, y: α/β };
  ctx.beginPath();
  ctx.arc(cx(eq), cy(eq), 5, 0, Math.PI * 2);
  ctx.fillStyle = '#f59e0b';
  ctx.fill();
}

7. Extensions: competition, three species, stochasticity

Logistic prey (Rosenzweig–MacArthur model)

Replacing the linear prey term with a logistic cap introduces a carrying capacity K:

This extension removes the structural instability and gives realistic closed limit cycles (Hopf attractors) for intermediate K values.

Competitive Lotka–Volterra (two prey species)

Two competing species sharing a resource:

Three-species food chains

Adding a top predator z (which eats y) creates a 3D ODE system. The trophic cascade can exhibit chaotic dynamics — the Hastings–Powell (1991) model showed that simple food chains can produce chaos without any stochastic driving, purely from the nonlinear interactions.

Stochastic extensions

Deterministic models give ensemble-average behaviour. Real populations are finite and discrete, making stochasticity essential for small populations:

• Demographic stochasticity: individual births and deaths are random events — use the Gillespie algorithm or τ-leaping for efficient simulation.
• Environmental stochasticity: parameters α, β vary randomly in time (e.g., seasonal forcing). Modelled by Itô stochastic differential equations.
• Spatial models: agents on a grid with local interactions. Can stabilise chaotic systems through spatial segregation of predator and prey.