How Often Does a Good Gene Lose? Measuring the Reach of Genetic Drift

A Number That Should Bother You

Take a mutation that genuinely helps. Nothing spectacular. Call it a copy of a gene that leaves its carrier one percent more descendants per generation than the version everybody else is carrying, which by the standards of real adaptation is an ordinary and thoroughly plausible advantage. Give it room. A thousand diploid individuals, two thousand copies of the gene in circulation. Then let it run.

It vanishes. In our runs it vanished 98.01 times out of a hundred.

That figure rests on 50,000 independent simulated populations, of which 996 ended with the mutation at every copy of the locus and 49,004 ended with it gone. The measurement is \(u = 0.01992 \pm 0.00063\). A better gene loses almost always. The reason has nothing to do with its being better. The first carrier may simply fail to reproduce for reasons that have nothing to do with its gene, and one copy is all the lineage ever had to spend.

The 41.8 deserves a moment. Run the same population with \(s = 0\) and a new copy fixed 286 times out of 599,700, a rate of 0.000477, which is what a gene is worth when it is worth nothing. So the one percent advantage did something enormous: it multiplied the odds of survival by more than forty, and it still left failure as the overwhelmingly likely ending. Both readings belong to the same number, and most of the confusion people carry about natural selection lives in the gap between them, in the distance between a mutation that is favoured and a mutation that survives.

What follows measures that survival probability across a grid of selection strengths and population sizes, sets the measurement beside the two classical formulas that claim to predict it, cell by cell for all forty, reports the two places where those formulas stop working, and closes with an account of what we could not model.

The Model, and What It Leaves Out

The club has no laboratory. No fruit flies, no bacterial cultures, no sequencer, no freezer. We have a room with computers in it, so our research is computation. Everything below came out of a program we wrote, which runs a mathematical model of a population and knows nothing whatever about the organisms such a population is meant to contain, or about the deaths and the failures to breed that the binomial draw quietly stands in for. Nothing here was observed in anything alive. Saying measured, we mean measured from our own replicates, the way you would measure the frequency of heads by flipping a coin many times instead of reasoning about coins.

The model carries two names. Sewall Wright and Ronald Fisher built the machinery of theoretical population genetics through the 1920s and 1930s [1][2][3], and the chain below is the smallest working piece of what they left. It works like this. A population holds exactly \(2N\) copies of one gene, forever. Every generation the whole adult population dies and is replaced, and each of the \(2N\) new copies is drawn at random, independently, from the pool of parental copies. Nothing else happens.

Selection enters as a thumb on the scale. If the mutant sits at frequency \(p\) and each mutant copy is worth \(1+s\) as many descendants as a wild-type copy, the pool for the next generation is weighted to

and the number of mutant copies in the next generation is a binomial draw, \(i_{t+1} \sim \mathrm{Binomial}(2N, p')\). The mutation starts at \(i = 1\). The process runs until \(i\) reaches 0, which we call lost, or \(2N\), which we call fixed, and both boundaries absorb, so a lineage that touches either one never comes back from it. Those two endings are the only ones. We run it a great many times and count.

The binomial draw is the entire source of randomness. Population geneticists call it drift, and the name flatters it. Drift acts on nothing. Drift is the arithmetic consequence of finitely many parents having finitely many offspring.

Two Formulas, Written Out

In 1927 Haldane asked what becomes of a single new beneficial mutation [4]. His population was large enough that the mutation's own frequency never matters. While the mutant stays rare, each copy leaves a number of copies next generation that is approximately Poisson with mean \(1+s\), independently of everything else, because a rare allele almost never meets another rare allele inside the same parent. A branching process, in other words. The probability \(\pi\) that such a lineage never goes extinct satisfies

Haldane then expanded this for small \(s\) and obtained the result everybody quotes:

A one percent advantage gives a two percent chance of establishment. Notice what the formula does not contain. No \(N\) appears in it. Haldane's population is effectively infinite, so the mutation escapes early or dies early, and nothing that happens afterwards can take back a lineage that has already escaped.

Kimura solved the finite problem in 1962 [5]. He treated the frequency as a diffusion with infinitesimal mean \(M(p) = sp(1-p)\) and variance \(V(p) = p(1-p)/(2N)\), which is the continuous limit of exactly the chain we simulate, and integrated the backward equation [7]. For a mutation starting at one copy out of \(2N\):

Kimura's answer does contain \(N\), and it contains it only in the combination \(4Ns\). Put in words rather than symbols: the diffusion says a mutation's fixation probability is the share of an exponentially warped frequency axis already lying behind its starting point, the warping factor being \(e^{-4Nsx}\), so selection reaches the answer only through how sharply that factor bends across the unit interval. Take \(4Ns\) large and the denominator goes to 1, leaving \(u \to 1 - e^{-2s}\), which for small \(s\) is Haldane's \(2s\). Take \(s \to 0\) and L'Hôpital gives \(u \to 1/(2N)\), the neutral answer, which the discrete chain also gives exactly by a symmetry argument: every one of the \(2N\) copies alive today is equally likely to be the ancestor of all copies in the distant future, so any particular one carries probability \(1/(2N)\).

So \(4Ns\) decides which regime a population lives in. It weighs selection against drift. Our grid sweeps it from 0.2 to 4000, four orders of magnitude, which covers Haldane's limit at one end, Kimura's neutral limit at the other, and the awkward middle where neither description is quite in charge.

Working Notes From the Club Table

The Arithmetic

The grid holds five values of \(N\) and eight values of \(s\). That makes 40 cells. Each cell runs an independent set of replicates. Its stream is spawned from the master seed 20240917 through numpy's SeedSequence. The generator is PCG64. Total replicates: 11,932,191. Total fixations: 90,675. Wall clock: 66.2 seconds on a laptop.

Figure 1 plots the measured fixation probability against \(s\) on log axes for three of the five population sizes, with Haldane's \(2s\) as the dashed straight line and Kimura's curves drawn underneath the points.

The left-hand end of the \(N = 50\) series carries the whole story. At \(s = 0.001\) the mutation has a real advantage. It fixes 1.176% of the time. Haldane says 0.2%. The neutral rate is 1.0%. Drift handles this mutation almost entirely; the advantage merely rides along.

The neutral rows are the code check. At \(s = 0\) the answer is \(1/(2N)\). No approximation enters anywhere. Our five measurements came out 0.009500, 0.002431, 0.000477, 0.000105 and 0.000054, against exact values of 0.010000, 0.002500, 0.000500, 0.000100 and 0.000050, so the agreement holds to the third significant figure in the smallest case and to the second in the largest. All five 95% intervals cover the exact value. The z scores run −1.12, −0.48, −0.80, +0.93 and +1.33.

Where Haldane's 2s Falls Apart

Haldane's formula carries no \(N\). So wherever \(N\) matters, Haldane must be wrong. The question is where that boundary sits and how sharply it cuts.

Figure 2 answers it by plotting our measurement divided by \(2s\) against \(4Ns\). Were Haldane right, every point would sit on the horizontal line at 1. The 35 cells with \(s > 0\) span four orders of magnitude in \(4Ns\), and they land almost exactly on the curve \(1/(1 - e^{-4Ns})\), which is what Kimura's formula reduces to when \(s\) is small enough that \(1 - e^{-2s} \approx 2s\).

The collapse is the point. \(N = 1000\) with \(s = 0.001\), \(N = 200\) with \(s = 0.005\) and \(N = 50\) with \(s = 0.02\) all share \(4Ns = 4\), and they gave 0.9687, 0.9560 and 1.0035 times Haldane's value. Inside the error bars, one number. Selection and drift never act separately; only their ratio reaches the answer.

Grouped by band, the mean absolute error of Haldane's \(2s\) across our cells was:

Where does Haldane break down? Below about \(4Ns = 4\), and the breakdown is fast. At \(4Ns = 2\) the error is 12.8%. At \(4Ns = 1\), 60.2%. At \(4Ns = 0.2\), 488.0%, and by then the formula has stopped describing the same phenomenon, because the mutation fixes at very nearly the neutral rate and its advantage has dropped out of the outcome altogether.

Note which direction the error runs. When \(4Ns\) is small, Haldane understates the fixation probability, because in a small population drift will sometimes carry a mediocre mutation all the way to fixation for no reason at all. The errors in the two right-hand bands run the other way. They are §7.

The Other Edge, at s = 0.1

The \(s = 0.1\) column misbehaves. Nothing else in the grid does. At all five population sizes our measurement sits below Kimura's prediction, by −3.6%, −2.7%, −2.7%, −3.1% and −3.2%, and the size of the population makes no difference at all to the sign of the miss. Five independent random streams produced those five signs. Noise would have scattered them.

Stouffer's method (sum of z divided by \(\sqrt{5}\)) pools the five z scores to −7.33. Every other column pools to something between −0.88 and +0.89, which is what five independent random streams look like when they have nothing in particular to report. So the deviation at \(s = 0.1\) is real. Nothing else in the grid deviates at all.

We spent a meeting looking for the bug. What we found was nothing. Kimura's formula solves a diffusion, and a diffusion is the limit of the discrete chain as \(s \to 0\) and \(N \to \infty\) with \(4Ns\) held fixed. It throws away terms of order \(s^2\). At \(s = 0.001\) those terms are invisible. At \(s = 0.1\) they are worth about three percent, which is exactly what we see.

The test sends you back to Haldane, past \(2s\), which is his answer after the small-\(s\) step, to the branching process he started from, solved without that step. For a lineage whose copy number is Poisson with mean \(1+s\) each generation, the escape probability \(\pi\) solves \(\pi = 1 - e^{-(1+s)\pi}\). We solve that numerically. Figure 3 sets all three predictions against our pooled measurements from the cells where \(4Ns \ge 100\), the region where finite \(N\) has stopped mattering and all three ought to apply.

Three numbers, one measurement. At \(s = 0.1\): \(2s = 0.2\) is 12.1% too high, \(1 - e^{-2s} = 0.181269\) is 3.0% too high, and \(\pi = 0.176134\) is 0.2% off. At \(s = 0.05\) the branching value is 1.2% low and Kimura is 0.3% low, so at that strength Kimura is the better of the two and the difference is inside our error bars anyway.

The honest summary: at \(s = 0.1\) the diffusion approximation is measurably worse than the branching one, and the gap is a real property of the limit Kimura took rather than a rounding artefact of ours. A three percent effect clears the noise only because we ran 11.9 million replicates. Neither formula was wrong when it was written. Haldane and Kimura were approximating different limits, and by \(s = 0.1\) those two limits have quietly stopped agreeing with each other by anything a reader could dismiss as rounding.

Watching an Estimate Settle

A Monte Carlo estimate claims a limit you never reach. So the defensible thing is to show how the estimate moved as the trials piled up, leaving the reader to judge whether it had settled or whether we stopped at a convenient moment.

The program keeps the full replicate-by-replicate record for three cells. It reports the running estimate at a geometric ladder of checkpoints. Because the replicates are independent and identically distributed, the first \(k\) of them form a valid smaller experiment, which makes the trace a real record of convergence rather than a redrawing of the final answer.

Sit with the neutral trace, which shows what an unlucky start looks like from the inside, and how long that start goes on resembling an answer. Five checkpoints in a row with zero fixations. Then one, then two, then a stretch at 0.47 times the truth, then a swing to 1.06, and only after about 200,000 replicates did the estimate look like it had a value. Stopping at 12,807 replicates would have had us report 0.000234 against a true 0.0005, wrong by a factor of two, with a straight face, and nothing on the screen at that moment would have looked the least bit wrong to us.

Hence the replicate count chosen from the expected answer instead of a round one. The rule we used is \(n \approx 300(1-u)/u\), which targets a relative standard error of 10% and then triples the count for safety, the tripling being the concession we made after the first run came back with three cells over target. The worst realised relative standard error anywhere in the 40 cells is 5.93%.

The neutral cells double as the check on the code, \(1/(2N)\) being exact. Figure 5 shows all five against the line.

How Long the Winners Take

Every fixing replicate also records how long it took. How long and how often are different questions with different answers. People confuse them constantly.

A neutral mutation that fixes takes a very long time. Kimura and Ohta gave the theory in 1969 [6]. A neutral mutation starting from a single copy takes on average about \(4N\) generations to fix, which for a population of ten thousand is a stretch of time nobody will ever watch. Our measurements, in generations:

The neutral rows land within 5.1% of \(4N\) at every size. The worst of the five, \(N = 200\) at −5.1%, sits 1.7 standard errors from the prediction, comfortably inside what three hundred fixing replicates can resolve. The selected rows are benchmarked against the deterministic logistic sweep, \((2/s)\ln(2N)\) generations, which comes from solving the noiseless version of the same dynamics, the trajectory a mutation would follow if the binomial draw were replaced at every step by its own expectation. In practice it is an upper bound. Our ratios climb from 0.43 towards 0.92 as \(4Ns\) grows, exactly as they should: the deterministic calculation only describes things well once the mutation has enough copies for randomness to stop mattering, and a successful mutation reaches that point faster than the deterministic path predicts, because the lineages that succeed are the ones that got lucky early.

The practical reading: in a population of ten thousand, a mutation worth 1% takes about 1300 generations to take over, if it takes over at all, with a spread from roughly 1100 to 1600. A neutral one needs about 40,000 generations, spread across 18,000 to 69,000. Drift is slow. Its slowness depends on nothing but the size of the population.

The Strongest Objection We Could Make

The objection runs like this.

Everything above compares two pieces of mathematics. Kimura derived his formula from the diffusion limit of the Wright-Fisher chain. We simulated the chain. The formula described it, which is what a correct derivation promises, and our forty cells could never have added anything beyond a confirmation of the promise already made. No discovery about populations lives in that. We checked that our implementation matches the model the theorem is about, and that numpy's binomial sampler draws the numbers it claims to draw. A student who wrote the model correctly could have predicted every number in the grid table to within our error bars before running anything, armed with a pocket calculator and the two formulas we were supposedly testing.

We think the objection largely correct. These parts of the study survive it.

The neutral check genuinely is only a code check, and we said so where it appears instead of dressing it up as a finding about populations, which it never was. Its value lies in what it would have caught. Publish 40 numbers and you should be able to show that at least one of them is exactly right, with the exactness coming from somewhere other than yourselves.

The Haldane breakdown in §6 says something quantitative about a formula still quoted in textbooks and in the primary literature with its condition left off, and a formula quoted without its condition is a formula waiting to be turned on a small population. "The probability of fixation is \(2s\)" holds when \(4Ns\) is above about 4. At \(4Ns = 0.2\) it is wrong by a factor of five. Our contribution is the number 4, measured rather than asserted, together with the shape of the failure on either side of it, which is a modest thing to add to a literature eighty years old, and it is ours. Patwa and Wahl's review makes the same point in a more general setting [8]. Our grid is a small checkable version of it.

The \(s = 0.1\) result in §7 is the part we would defend hardest. Nobody handed us the prediction that Kimura would be 3% high there. We found it as an anomaly, went looking for a bug, failed to find one, and then identified the correct explanation and tested it against a third formula. That sequence is the useful part. Any club with a laptop can run it.

A second objection resists us. Our \(N\) is at once the census size and the effective size. Real populations keep those apart, often by a large factor [11][16]. A species of ten million individuals may have an effective size of ten thousand, and the effective number is the one that governs, because old bottlenecks and unequal reproductive success both subtract from the census count. So the row of our table that fits a real organism is rarely the row you would guess by counting organisms. Everything we report about \(4Ns\) is correct. Getting \(N\) for a real population lies outside what we did.

What This Does and Does Not Tell You

A beneficial mutation is a lottery ticket. The odds are splendid by lottery standards and dreadful by any standard a person would apply to something they cared about. At \(s = 0.01\) the ticket wins about one time in fifty. Most adaptations that ever appeared in any lineage were lost immediately, and the ones still around to be studied are survivors of a filter that discards 98% of equally good candidates, so the history of adaptation we are able to read is a record of luck at least as much as of merit, and the far larger history of the unlucky left nothing behind to be read.

\(4Ns\) is the number to compute first. Below about 1, a population cannot tell a beneficial mutation from a neutral one, and its fate is decided by drift; the consequences of this for small and endangered populations are the subject of a large literature [11][16]. Above about 4, selection is in charge and Haldane's \(2s\) is a good working number.

The study says nothing about a population changing size, about mutations competing, about a shifting environment, about dominance or recessivity, and nothing about how to obtain \(N_e\) for a real species. About any particular organism it says nothing whatsoever. We did not look at one.

We would rather a reader took away the shape of the check than any single number, since the numbers belong to our grid and to the seed we happened to choose, while the shape belongs to anybody willing to set a prediction beside a measurement and look at the two of them honestly. A prediction, a measurement with a stated uncertainty, the two printed side by side for all 40 cells, and when one column disagreed, a search for the reason instead of a rounding away of the disagreement.

Where a different choice would change the answer

Several choices could have gone otherwise. Saying which ones would have moved the answer seems the least we owe a reader.

Genic against dominant. Each mutant copy was worth \(1+s\), so the heterozygote sits exactly halfway. A fully dominant mutation would show its whole advantage at once while rare. The fixation probability would roughly double. A fully recessive one would be invisible to selection while rare, since almost every carrier is a heterozygote, and its probability would collapse towards the neutral value. Of every choice in the model this one moves the headline number most, and we did not vary it, so read the 98% as a statement about genic selection and about nothing else.

Wright-Fisher against Moran. The Moran model replaces one individual at a time instead of the whole population, and produces a fixation probability of the same form with \(4Ns\) replaced by \(2Ns\) [10]. We would have got the same picture on a shifted axis.

Constant against changing size. Otto and Whitlock showed that in a growing population a beneficial mutation fixes more often than \(2s\) and in a shrinking one less often [9], which matters for anything recovering from a bottleneck or sliding towards one. Population size enters only through \(4Ns\), so a changing \(N\) makes the answer depend on the whole trajectory rather than on any one value along it. A selection coefficient that varies over time does something similar [15]. Our constant-\(N\) result is the baseline both of those departures start from.

One locus against many. Of our assumptions it fails most often, the supposition that nothing else in the genome is under selection at the same moment and in the same organisms. In a large asexual population, beneficial mutations arise faster than they fix and compete with each other, which Gerrish and Lenski called clonal interference [12], a regime our one-locus grid cannot represent at all. Under interference the fixation probability of a given mutation falls below \(2s\), sometimes far below, and the mutations that do fix are drawn from the strong tail of the distribution rather than from all of it. Lineage-tracking experiments in yeast have since watched it happen, lineage by barcoded lineage [14]. None of that is in our model.

Reproducing this

Everything in this article comes from one file and one command. You need Python 3.12 and numpy; no other package is used.

python drift-versus-selection.py > drift-versus-selection-output.txt

Expect about 70 seconds on a modern laptop. Ours took 66.2 s, on numpy 2.4.2 and Python 3.12.3. Peak resident memory measured 240 MB, set by the \(N = 10{,}000\), \(s = 0\) cell and its 5,999,700 replicates, that being the cell holding the largest live array of replicate states. The master seed is 20240917, hard-coded at the top of the file, and each of the 40 cells draws an independent stream from it through numpy's SeedSequence.spawn, so the whole output is deterministic. Numbers differing from ours by more than the printed standard errors mean something is wrong, and we would like to hear about it. The interactive model runs a smaller version of the same simulation in your browser.

References