A Population That Grows on Average Can Still Die Out

What the question actually asks

Begin with the average. Every individual in this population leaves 1.15 offspring and then dies, so the population grows by 15 per cent a generation and doubles in five generations. An average, however, is a summary of many outcomes rather than a promise about any one of them. No individual ever leaves 1.15 offspring. Offspring arrive in whole numbers: zero, one, two, six. The average is only what those whole numbers come to once you have added them all up and divided, and nothing prevents a small population from drawing a run of low numbers long before the average has any chance to assert itself. Worse, the drawing repeats. Twenty individuals who between them happen to leave fourteen offspring have started downhill, and those fourteen must then take their own draw, from the same distribution, with no memory of the twenty that produced them. Nothing in the process remembers that the population has already had a bad generation, and no compensating good generation is owed to it, which is the whole difficulty with small numbers and the reason the arithmetic is worth taking slowly.

The question therefore has two halves. How large must a population with positive mean growth be before reproductive randomness loses the power to kill it, and how far does that size move when reproduction becomes more variable while the average stays exactly where it was?

The second half is what drew us in. Two populations can share an average and share nothing else about how reproduction is distributed among their members, and the difference hides entirely inside the spread. In one, almost everybody leaves one or two offspring. In the other, almost everybody leaves none, and a rare individual leaves thirty. The long-run arithmetic is identical. The short-run risk is not. How large that gap is, and whether a single number as crude as the variance can predict it, is the measurement this article was built to report.

The experiment is the program

The club has no laboratory. No field site, no licence to handle animals, no budget. We have a room, a table and a set of laptops. So when this article says experiment, it means a program that was written, seeded and run, and every number quoted below is printed output from that program rather than an observation of anything alive. Nobody counted a beetle. Nobody watched a population of anything.

The limitation is real, and it governs every claim this article is entitled to make. A simulation cannot tell you whether beetles reproduce the way our model says they do. It can tell you what follows, exactly and checkably, from a set of assumptions you are free to read in the source, and where those assumptions are wrong about the world our numbers are wrong about the world in precisely the same way. Our only defence is to write every assumption down in a place where you can find it and attack it.

The code lives at analysis/extinction-by-bad-luck.py. Every random number comes from a numpy PCG64 generator seeded from the master seed 21091847, with one spawned child stream per parameter cell, so the run reproduces exactly, and the output file linked at the top of this page is the run that produced every figure and table here.

The two models, and what they leave out

Model A is a Galton-Watson branching process, the oldest model in this subject: Watson and Galton built it in 1875 to work out how often aristocratic surnames die out [1]. Generations do not overlap. Every individual alive in a generation independently produces a random number of offspring, drawn from a fixed distribution with mean m, and then dies; the survivors of that draw face the identical rule in the generation after. We fixed m = 1.15 throughout and varied only the shape of the distribution, so that the five variance levels differ in nothing else:

The sum of n independent draws from each of these families is itself a draw from one named distribution, so a whole generation of any size is produced with a single vectorised call. Closure under addition is an exact identity rather than an approximation, and it is the only reason a study of this size finishes in under three minutes on a laptop.

Model B runs in continuous time. Generations overlap and nothing is synchronised. Each individual gives birth at per-capita rate λ and dies at per-capita rate μ. We held the growth rate r = λ − μ at 0.15 and raised the turnover λ + μ, which is the rate at which demographic variance accumulates, and simulated the process exactly through its embedded jump chain with exponential waiting times, the standard construction going back to Kendall [3] and used as a simulation method by Gillespie [5].

What both models leave out is a good deal longer than what they contain, and the omissions are not accidental. No carrying capacity constrains them. Survivors simply grow for ever at 15 per cent a generation, which no population in the world has ever done. Environmental variation is absent too: the offspring distribution never has a good year or a bad one, and that single omission matters more than anything else in this article, because for real populations environmental stochasticity usually dominates once numbers are above a few dozen [7, 9]. Individuals reproduce independently. Mate finding, crowding, the Allee effect and competition for a shared resource are all absent, as are age structure, sex, space and immigration, and the list could be made longer than the model. Every one of those omissions was chosen rather than overlooked, and every one of them makes the model kinder to the population than the world would be, so the numbers below are a best case and should be read as one. What remains is demographic stochasticity, alone. Isolating it is artificial, and isolation is nonetheless the only way to find out how large the term is on its own.

Checking the simulator against answers we already know

We chose this model class for one reason: it has exact answers, and exact answers let a simulator be caught out rather than merely admired. A probability generating function is a bookkeeping device, and there is nothing deep about it at all. It packs the whole offspring distribution into a single function of a dummy variable s, and composing it with itself once corresponds exactly to letting one generation pass. Write f(s) for that function. The probability that a lineage descended from one individual is extinct by generation G is then exactly the G-fold composition of f with itself, evaluated at zero, and the probability of eventual extinction is the smallest non-negative root of f(s) = s [1, 4].

So every one of our 210 simulated cells in Model A has an exact partner computed a different way, and the comparison is a test the simulator can fail rather than a reassurance it cannot. We computed the root two ways: by bisection on f(s) − s, and, where a closed form exists, by solving that closed form independently. For geometric offspring the closed form is q = 1/m, and our bisection reproduces it with a difference of −4.4×10⁻¹⁶, which is a few times the spacing between adjacent doubles. For Poisson offspring the root is −W(−me^−m)/m, where W is the Lambert function, which we solved with our own Newton iteration; the difference from bisection is −1.2×10⁻¹⁵, and substituting each root back into its own generating function leaves a residual of at most 1.1×10⁻¹⁶.

The simulation was then compared with that analytic answer, cell by cell, across the whole grid. If the simulator is correct, the quantity (P_simulated − P_exact) divided by the binomial standard error is a draw from a standard normal, and the collection of them across cells should look like a sample from one. We separated out the cells where fewer than ten extinctions were expected, because the normal approximation is worthless there, and tested those instead with an exact Poisson tail probability.

Passing looks like that. The measured curve is the analytic curve, to within the noise that 20,000 replicates leave behind, at every starting size and under every offspring law we tried.

Which raises the question of whether 20,000 replicates is enough, and Figure 1 answers it for a single cell by running that cell out to 200,000 and recording the estimate as the trials accumulate. At 100 replicates the estimate was 0.050 with a 95 per cent interval of ±0.043, which is useless; at 20,000 it was 0.0620 ± 0.0033; at 200,000 it was 0.06139 ± 0.00105 against an exact value of 0.061100. The interval shrinks as the square root of the trial count. The estimate stays inside it. Convergence means nothing more than that.

The arithmetic

The last column is the flattest statement we can make without any model at all: the smallest size on our grid where the whole 95 per cent interval of the measured extinction probability sits below 5 per cent. For the lottery that size is 400 rather than 303, because our grid jumps from 300 straight to 400 and the interval at 300 straddles the line.

Notes from the table

How long it takes a population to die

Among the runs that ended at zero, we recorded the generation in which the last individual died, and the pattern was not the one the club expected. More variance kills faster. It kills earliest of all in the populations a reader would have picked out as the safest ones.

With Poisson offspring and a starting size of 2, the populations that died took 5.24 generations on average, with a median of 3; under the lottery distribution at the same starting size the mean is 1.21 generations and the median is 1. Nearly all of those populations died in the first generation, because in a lottery distribution the most likely thing for two individuals to do is leave nothing at all. An individual leaves 30 offspring with probability 1.15/30 = 0.0383 and otherwise leaves none, so two of them leave nothing between them with probability 0.9617² = 0.9248.

Starting size pushes extinction later. Geometric offspring starting at 2 give a mean of 4.17 generations. Starting at 20 gives 13.34. The lottery starting at 500 gives 14.94. Larger populations that fail take longer to fail, because they have to lose more ground, and the conditional distribution has a long right tail in every case: the 95th percentile is regularly two to three times the median. Figure 4 shows the means with their interquartile ranges.

The practical reading follows. If demographic bad luck is going to take a small population, it takes it quickly, within a handful of generations, while the population is still small enough for anyone monitoring it to think the numbers are just noisy. A population that survives its first twenty generations at this growth rate has mostly escaped: for a single Poisson lineage the probability of extinction is 0.7413 by generation 20 and 0.7510 by generation 100, so the last eighty generations contribute about one per cent of the risk.

The strongest objection we could find

What follows is the case against this article, assembled by the members who wanted the lottery distribution thrown out, and put as strongly as we know how to put it. We have kept it in the body of the article rather than in a footnote, because a study whose strongest criticism appears only in small type at the bottom of a page is a study asking not to be criticised.

The result is arranged. We held the mean fixed and increased the variance by changing the distribution to one that is more and more skewed, and skewness is doing the work rather than variance. In the lottery distribution 96.2 per cent of individuals leave no offspring at all, which is not variance so much as sterility. Of course such a population dies: it dies because almost everybody is sterile and the survival of the lineage depends on a rare jackpot, which is a completely different biological situation from a population where reproduction is merely noisy. Putting the two on a single axis labelled variance invites the reader to believe that one quantity governs extinction risk, when two or three of them do.

We think this objection is partly right, and the pleasant thing about it is that we can measure how right. If variance alone governed the threshold, the ratio of threshold to variance would be constant across our five laws, and it is not quite: the ratio runs 9.10, 8.92, 8.67, 8.44, 9.13, a spread of about 8 per cent. Variance predicts the threshold well but not perfectly. What is left over depends on the rest of the distribution's shape. Kendall and Fox make a related and sharper point in the other direction: variation among individuals in their expected reproduction, as opposed to variation in outcomes, can reduce the effective demographic variance rather than raise it [8, 11], so the mapping from biological heterogeneity to the variance in our model is not automatic.

A second form of the objection holds that discrete, non-overlapping generations are an artefact of the model rather than a fact of biology. A population in which everybody reproduces at once and then dies is a fiction for most organisms, and the synchronisation itself could be what manufactures the risk. That one we could test, which is why Model B exists. In continuous time, with overlapping generations and no synchronisation at all, holding the growth rate at 0.15 and raising the turnover from 1.15 to 6.15 moves the threshold from 11.42 to 61.40, and the ratio of threshold to turnover is 9.93, 9.97, 9.98, 9.98. The scaling survives the assumption being dropped. It survives more cleanly here than in the discrete model. Figure 5 shows those curves.

The third form of the objection is the one we cannot answer, because environmental variation was left out of both models entirely, and for real populations it is usually the dominant term above a few dozen individuals, scaling with population size quite differently from demographic noise [7, 9]. Our thresholds are therefore lower bounds on how big is big enough, and they are lower bounds derived from a model that deliberately removed the thing that usually matters most.

Our numbers against the published ones

The first comparison is to century-old theory. Haldane, working on whether a new advantageous type survives its own rarity, gave the survival probability of a single lineage as approximately 2s/V, where s is the excess growth and V the offspring variance [2]. Haldane derived it for small s, and ours is 0.15, which is not especially small. Comparing it with our exact roots gives ratios of 1.0477, 0.9984, 0.9302, 0.8683 and 0.9187 across the five laws, so the century-old formula is right to within 0.2 to 13.2 per cent over a thirty-fold range of variance. That formula is why the linear scaling in Figure 3 was predictable before we ran anything, and what the simulation adds is the size of the error in the approximation and the shape of the whole curve rather than its rate alone.

The second comparison is to measured offspring distributions. The negative binomial with dispersion parameter around 0.1 to 0.25 is not a mathematical curiosity: it is roughly what is estimated for the number of secondary infections per case in several directly transmitted diseases, where a small minority of individuals cause most of the transmission [12], and estimates for COVID-19 put the dispersion near 0.1 [16]. Our r = 0.25 level sits inside that measured range, at variance 6.44. Its threshold is 54 individuals. The numbers do not transfer directly, since a disease with an R₀ of 2 or 3 is not a population growing at 15 per cent, but the shape of the distribution is real and measured rather than invented for this study.

The third comparison is the awkward one. Conservation biology has argued about minimum viable population size since 1981, when Shaffer defined the term [6], and the published numbers dwarf anything in this article. A meta-analysis of 212 species put the median published estimate at 4,169 individuals with a 95 per cent interval of 3,577 to 5,129 [14], and the genetic rules of thumb now argued for are an effective size of 100 against inbreeding and 1,000 for evolutionary potential [15]. Our thresholds run from 10 to 303, and the distance between those numbers and theirs measures what we left out. Demographic bad luck alone stops being the binding constraint at a few hundred individuals, and everything above that number in the published estimates is environmental variance, catastrophes, genetics and habitat. Lande set out why the terms scale differently [7]. The published minimum viable population sizes and ours are answers to different questions, and setting them side by side is only useful once somebody says which question each one answers.

The fourth comparison is to field data on small populations. Matthies and colleagues followed 359 populations of eight threatened plants in northern Germany for ten years and found 27 per cent of them gone, with extinction strongly concentrated in the smallest populations [10]. Fagan and Holmes assembled ten vertebrate populations monitored to extinction and found that time to extinction scales logarithmically with abundance [13], and our conditional extinction times rise roughly linearly with the logarithm of starting size, at between three and five generations per e-fold. We record that the form matches while being clear that the settings do not: their populations were declining and ours are growing on average, and their units are years of real organisms while ours are generations of a model.

Where a different choice would have changed the answer

Four choices in this study move the headline number, and every one of them deserves to be named out loud rather than left in the code.

The mean growth rate matters at least as much as the variance does, and we held it fixed at 1.15 throughout. Recomputing the exact threshold at other means, for the same three laws: at m = 1.02 the Poisson threshold is 65.5 and the lottery threshold is 1,901. At m = 1.10 they are 15.5 and 448. At m = 1.50, 3.4 and 98.8. So any single number quoted from this study is a number about populations growing at 15 per cent per generation, and a population barely above replacement needs an order of magnitude more individuals for the same safety.

We picked the 5 per cent criterion for no better reason than that it is the convention. A 1 per cent criterion would multiply every threshold by ln(100)/ln(20), which is 1.54. A 10 per cent criterion would multiply them by 0.77. The threshold is proportional to the logarithm of the risk you are willing to accept, so it is not a cliff, and nothing changes character at 5 per cent.

The 100-generation horizon turned out not to matter. We know that only because we checked. Extinction probability within G generations, for a single lineage, is 0.6210 at G = 5 and 0.7510 at G = 100 for Poisson offspring, and the gap between 100 generations and eventual extinction is 7.5×10⁻⁸. At a horizon of 20 generations the answers would have been visibly different, and at 5 generations this would have been an altogether different study.

Labelling these five laws by variance alone hides an 8 per cent spread, described in section 8. If we had chosen five different distributions with the same five variances, the thresholds would have come out slightly different, and the fitted constant of 8.83 would have moved by something like that amount. The slope of 1.0014 is far more stable than the constant, because the slope is a statement about scaling while the constant is a statement about one particular set of distributions.

Reproducing this

Everything in this article comes from one file and one command. Python 3.12.3 with numpy 2.4.2 was what we used. From the project root:

python analysis/extinction-by-bad-luck.py > analysis/extinction-by-bad-luck-output.txt

Expected runtime is about three minutes on a laptop, and our own run printed 174.5 seconds before it stopped. The script imports only math, sys, time and numpy, writes nothing except what it prints, and with the master seed left at 21091847 the output should match ours line for line, including every z-score quoted above. If it does not, the seed handling or the numpy version has changed, and the analytic columns are the ones to trust, because they contain no randomness at all.

The raw output file is linked at the top of this page, and the last section of it is a machine-readable block that every figure here was plotted from. The interactive model runs a smaller version of the same two experiments in the browser, so you can set the mean, the variance and the starting size yourself and watch the lines fall to zero.

References