Reaction–Diffusion and Spatial Spread
An invading species, a spreading epidemic, and an advantageous gene all share the same problem: something grows locally while also moving through space. Reaction–diffusion equations couple these two ingredients — local dynamics and random movement — and their signature result is a front that advances at a predictable speed.
as the population invades empty space." />
The general form
A reaction–diffusion model tracks a density that changes for two reasons. Locally it grows or decays according to a reaction term , and it spreads out by diffusion, the tendency of random movement to flatten gradients. Combining them gives where is the diffusion coefficient (units of area per time) and is the Laplacian. The partial derivative is the rate of change in time at a fixed point in space, while (in one dimension) measures the local curvature of the density. Where dips below its neighbors the curvature is positive and diffusion fills the dip in; where peaks the curvature is negative and diffusion erodes the peak.
The diffusion term alone, , is the heat equation: it smooths any initial profile toward uniformity. The reaction term alone, , is an ordinary differential equation describing growth at a single location. Their combination is what lets a locally growing population push outward into empty territory.
The Fisher–KPP equation
The most famous choice takes to be logistic growth, giving the Fisher–KPP equation Here is scaled so the carrying capacity is , is the intrinsic growth rate, and is the diffusion coefficient. R. A. Fisher introduced it in 1937 to model the spread of an advantageous gene through a population distributed along a line, and Kolmogorov, Petrovsky, and Piskunov analyzed it the same year.
The equation has two spatially uniform equilibria: the empty state and the invaded state . Treating the reaction term as a local ODE, , so makes the empty state unstable while makes the invaded state stable. Any small seed of population therefore grows, and diffusion carries the growth outward.
The traveling wave
The key result is that the solution settles into a traveling wave: a fixed front profile that connects behind it to ahead of it and moves at constant speed without changing shape. Writing for a front moving right at speed , the asymptotic speed selected from a compact (localized) initial condition is This minimal speed comes from linearizing at the leading edge, where is tiny and the equation reduces to ; requiring a non-oscillating (non-negative) front picks out . Faster growth or faster movement both speed the invasion, but only through their geometric mean under the square root, so quadrupling or merely doubles the speed.
A worked example
Suppose a beetle invades with growth rate and diffusion coefficient . The asymptotic front speed is Over a decade the invasion front therefore advances roughly km. If a management program halves the effective growth rate to , the new speed is — a reduction by a factor of , not a factor of two, because the speed depends on .
Simulation
We integrate the 1D Fisher–KPP equation by the method of lines: discretize space on a grid, approximate by the second difference , and step forward with explicit Euler. Starting from a small patch at the left, the front should advance at speed . For stability the time step must satisfy .
R
r <- 0.4; D <- 0.25
dx <- 0.5; dt <- 0.1 # D*dt/dx^2 = 0.1 <= 0.5, stable
L <- 200; n <- L / dx + 1
x <- seq(0, L, by = dx)
u <- as.numeric(x < 5) # invaded patch on the left
steps <- 2000 # t = 200
for (s in 1:steps) {
lap <- c(0, diff(diff(u)), 0) / dx^2 # Neumann (zero-flux) ends
u <- u + dt * (D * lap + r * u * (1 - u))
}
front <- max(x[u > 0.5]) # position of u = 0.5 level set
front / (steps * dt) # ~0.63 vs 2*sqrt(r*D) = 0.632
Python
import numpy as np
r, D = 0.4, 0.25
dx, dt = 0.5, 0.1 # D*dt/dx^2 = 0.1 <= 0.5, stable
L = 200.0
x = np.arange(0, L + dx, dx)
u = (x < 5).astype(float) # invaded patch on the left
steps = 2000 # t = 200
for _ in range(steps):
lap = np.zeros_like(u)
lap[1:-1] = (u[2:] - 2 * u[1:-1] + u[:-2]) / dx**2
u = u + dt * (D * lap + r * u * (1 - u))
front = x[u > 0.5].max() # position of u = 0.5 level set
print(front / (steps * dt)) # ~0.63 vs 2*sqrt(r*D) = 0.632
0.6075
Julia
r, D = 0.4, 0.25
dx, dt = 0.5, 0.1 # D*dt/dx^2 = 0.1 <= 0.5, stable
L = 200.0
x = collect(0:dx:L)
u = Float64.(x .< 5) # invaded patch on the left
steps = 2000 # t = 200
lap = zeros(length(u))
for _ in 1:steps
@views lap[2:end-1] .= (u[3:end] .- 2 .* u[2:end-1] .+ u[1:end-2]) ./ dx^2
u .+= dt .* (D .* lap .+ r .* u .* (1 .- u))
end
front = maximum(x[u .> 0.5]) # position of u = 0.5 level set
println(front / (steps * dt)) # ~0.63 vs 2*sqrt(r*D) = 0.632
Why it matters
Reaction–diffusion is the workhorse model for anything that spreads while it grows. Skellam used it to explain the spread of muskrats across Europe from a handful of escapees, and the same law predicts how fast an invasive species colonizes new range. In epidemiology, adding diffusion to an SIR model turns a well-mixed outbreak into a traveling wave of infection sweeping across a landscape. When space is patchy rather than continuous, metapopulation models play the analogous role, and when the reaction couples two interacting species the same framework produces the stationary Turing patterns rather than moving fronts.