The Collapse Bench

Both models below are ours. Neither of them appears in Prentice et al. (2025), which is a microbiology paper and contains no population model at all. Every parameter except one is invented but plausible, chosen so that a healthy reef looks roughly like a real Northeast Pacific kelp forest. The exception is the sea star attack rate of 0.68 urchins per star per day, which comes from a published feeding study. Nothing on this page is a measurement.

Model 1. One reef, and what it takes to get it back

Kelp grows and gets grazed. Urchins graze and get eaten. Sea star density is imposed by hand rather than solved for, because sunflower stars spawn into the plankton and their local numbers are set by regional larval supply. Disease enters only as a collapse in that imposed density. The timeline runs baseline, then the wasting shock at year 5, then sixty years of whatever the reef becomes, then a restocking programme at year 65.

The slider that actually decides the outcome is the third one, the drift-algae subsidy. It is how much food an urchin can find on a bare reef when there is no kelp. Set it to zero and the barren starves itself back into a forest. Set it high and the barren is permanent no matter how many sea stars you add.

Recovery here means kelp back above half its pre-outbreak value inside the forty years after restocking. The threshold is found by bisecting the restocking multiple, the same method the Python script uses. On the default settings this page reproduces the script exactly: a steady state of 9.2190 kg/m2 and 0.5094 urchins/m2, and a threshold of 1.551.

Model 2. An epidemic walking down a coastline

Three thousand kilometres of coast, chopped into 100 patches of 30 km. Each patch holds susceptible, infectious and dead sea stars. Transmission between patches falls off exponentially with distance, and the rows of that kernel are normalised, so the basic reproduction number is exactly the transmission rate divided by the recovery rate. Nothing recovers in any meaningful sense: the infectious period is the time from exposure to death.

The point of the model is the gap it opens up. Tune the transmission rate until the front moves at the observed 8.22 km per day and the reproduction number lands near 1.57, which predicts about 62% mortality. The real epidemic killed 90% or more. Watch what you have to do to the other two sliders to close that gap.

One honest difference from the script. To stay responsive in a browser this version uses 100 patches of 30 km and a second-order integrator, where the Python script uses 150 patches of 20 km and fourth-order Runge-Kutta. At the script's fitted transmission rate of 0.1307 per day, the coarser grid here reports 7.96 km per day rather than 8.22, and pressing Fit lands on about 0.1328 rather than 0.1307. Final epidemic size agrees to within a twentieth of a percentage point either way.

The wave speed pins down the product of the transmission rate and the kernel length. Many different reproduction numbers reproduce the same 8.22 km per day, which is a real identifiability limit rather than a bug. Try holding the speed fixed while you move the kernel length: the reproduction number will move and the final epidemic size will move with it, and the model will not tell you which combination is the true one.