Lotka–Volterra Predator–Prey Model

The equations

With prey x and predators y (numbers or densities), four positive rates define the model:

• α — prey growth rate without predators (exponential, no food limit).
• β — rate at which encounters remove prey; encounters are proportional to xy (mass action).
• δ — rate at which encounters add predators.
• γ — predator death rate without prey.

Setting both derivatives to zero gives the coexistence equilibrium x* = γ/δ, y* = α/β, where the dashed lines (nullclines) cross. The prey level at equilibrium is set by the predator's parameters and vice versa. Averaged over a full cycle, the populations equal these equilibrium values exactly: integrating d(ln x)/dt = α − βy over one period gives ⟨y⟩ = α/β, and likewise ⟨x⟩ = γ/δ.

Closed orbits and neutral stability

Dividing one equation by the other separates the variables and yields a quantity that never changes along a solution:

V has a single minimum at the equilibrium and rises in every direction, so its level curves are nested closed loops. Each solution is trapped on the loop it starts on: the populations return exactly to their initial values, again and again. Move the initial prey or predators and the lab draws a different loop; the gray curves are neighboring loops for the same parameters.

This is neutral stability, and it is not the same thing as a limit cycle. A limit cycle is a single loop that nearby trajectories approach, so the amplitude is a property of the system. Here there is no preferred loop: a disturbance moves the system to another orbit and it stays there, with a new amplitude, forever.

That makes the model structurally fragile. The closed orbits exist only because V is exactly conserved, and almost any change to the equations breaks the conservation law. Give the prey a carrying capacity and the orbits spiral in to a steady state; add predator saturation and they can spiral out or settle on a genuine limit cycle. Real cycles, such as those of the agent-based wolf–sheep model, are noisy and have amplitudes set by the system, not by the starting point.

Near the equilibrium the orbits are nearly ellipses and the period approaches 2π/√(αγ), from linearizing the equations. Wider orbits take longer: they spend a long time near the axes where one population is scarce and changes slowly.

How the simulation solves it

The lab integrates the equations with the classical fourth-order Runge–Kutta method (RK4). The number of steps per output interval is chosen from the current per-capita rates, so the fast crashes of a wide orbit get small steps and a quiet stretch gets large ones. Because V should be exactly constant, its largest relative drift over the run is a direct measure of numerical error. With RK4 it stays tiny; the simple Euler method would instead let V grow every step and turn the closed loops into outward spirals — see Euler vs Runge–Kutta.

The simulated period is the mean time between upward crossings of the prey population through γ/δ, interpolated within steps.

Example result

With α = 1, β = 0.05, δ = 0.01, γ = 0.5 and 40 prey with 9 predators at the start, the equilibrium is (50.00, 20.00). The prey swing between 14.66 and 119.62, the predators between 8.81 and 38.08. One cycle takes 9.472 time units — longer than the small-oscillation period 8.886, because this orbit is wide. Over 50 time units the conserved quantity V = −3.1917 drifts by only 3.3e-13 of its value.

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator. The excerpts show the core of each.

// dx/dt = αx − βxy (prey),  dy/dt = δxy − γy (predators)

export function derivatives(x, y, p) {
  return [p.alpha * x - p.beta * x * y, p.delta * x * y - p.gamma * y];
}

export function rk4Step(x, y, h, p) {
  const k1 = derivatives(x, y, p);
  const k2 = derivatives(x + (h / 2) * k1[0], y + (h / 2) * k1[1], p);
  const k3 = derivatives(x + (h / 2) * k2[0], y + (h / 2) * k2[1], p);
  const k4 = derivatives(x + h * k3[0], y + h * k3[1], p);
  return [x + (h / 6) * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]), y + (h / 6) * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1])];
}

// Along every exact solution this quantity is constant, so each orbit is a closed curve.

export const conserved = (x, y, p) => p.delta * x - p.gamma * Math.log(x) + p.beta * y - p.alpha * Math.log(y);

def derivatives(x, y, p):
    """dx/dt = alpha*x - beta*x*y (prey), dy/dt = delta*x*y - gamma*y (predators)"""
    return (p["alpha"] * x - p["beta"] * x * y, p["delta"] * x * y - p["gamma"] * y)


def rk4_step(x, y, h, p):
    k1 = derivatives(x, y, p)
    k2 = derivatives(x + h / 2 * k1[0], y + h / 2 * k1[1], p)
    k3 = derivatives(x + h / 2 * k2[0], y + h / 2 * k2[1], p)
    k4 = derivatives(x + h * k3[0], y + h * k3[1], p)
    return (x + h / 6 * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]), y + h / 6 * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]))


def conserved(x, y, p):
    """Constant along every exact solution, so each orbit is a closed curve."""
    return p["delta"] * x - p["gamma"] * math.log(x) + p["beta"] * y - p["alpha"] * math.log(y)

Run the full model from a terminal:

node simulation.js --alpha 1 --beta 0.05
python simulation.py --alpha 1 --beta 0.05

Download the model

Assumptions & limitations

• Unlimited prey food. Without predators the prey grow exponentially forever; there is no carrying capacity.
• Mass-action predation. Encounters are proportional to xy, and a predator never gets full or spends time handling prey.
• Continuous and deterministic. Populations are real numbers with no randomness, so they can dip to tiny values and recover; real populations that small would likely go extinct.
• No structure. No age, space, seasons or time delays; every individual is identical.
• Structurally fragile. The closed orbits depend on the exact form of the equations; small changes to the model change the qualitative behavior.

Related simulations

• OS-027Predator–preyWolves, sheep and grass on a grid: noisy population cycles, a phase loop and the risk of extinction.
• OS-028Logistic growthThe S-curve of a population meeting its carrying capacity, against unchecked exponential growth.
• OS-029SIR modelSusceptible, infected, recovered: the three-compartment ODE behind epidemic curves and R₀.
• OS-035Harmonic oscillatorA mass on a spring with friction and a periodic push: damping regimes, phase portraits and resonance.

Background: Euler vs Runge–Kutta: Time Stepping in Simulations

References

1. Lotka, A. J. (1925). Elements of Physical Biology. Williams & Wilkins.
2. Volterra, V. (1926). Fluctuations in the abundance of a species considered mathematically. Nature, 118, 558–560. doi.org/10.1038/118558a0
3. Murray, J. D. (2002). Mathematical Biology I: An Introduction, 3rd ed. Springer.
4. Strogatz, S. H. (2015). Nonlinear Dynamics and Chaos, 2nd ed. Westview Press.