Skip to content

Original ResearchVolume 1, Issue 1, Fall 2024

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

A mutation worth a one percent fitness gain still disappears 98 times out of 100. We ran 11.9 million Wright-Fisher replicates to measure the fixation probability across a grid of selection coefficients and population sizes, and to find out exactly where Haldane and Kimura stop working.

Written jointly by the Science Journaling Club

Study
Computational study
Review
Peer-edited by the club review board
Length
12 sectionsabout 31 min

Abstract

A mutation that makes its carrier one percent fitter still disappears almost every time it appears. We measured how often. The instrument was a Wright-Fisher simulation the club wrote and ran on a laptop, and no organism entered the study at any point, because here the computation is the experiment. Across a grid of eight selection coefficients (s=0s = 0 to 0.10.1) and five population sizes (N=50N = 50 to 10,00010{,}000) we ran 11,932,191 independent replicates from a single new mutant copy and counted 90,675 fixations, holding the relative standard error to no worse than 5.93% in any of the 40 cells. For N=1000N = 1000 and s=0.01s = 0.01, a single new copy fixed in 1.992% ± 0.063% of replicates. Loss took the other 98.01%. The neutral case recovered the exact answer 1/(2N)1/(2N) at all five population sizes, the one corner of the study where the theory leaves no room for approximation at all, and every 95% interval covered the true value, the largest z score being 1.33. Haldane's 2s2s is excellent while 4Ns≳44Ns \gtrsim 4. Below that it fails badly, falling short of the measured value by 488% at N=50,s=0.001N = 50, s = 0.001; its mean absolute error is 255.6% for 4Ns<14Ns < 1 and 2.5% for 4Ns4Ns between 4 and 20. Kimura's diffusion formula holds to a few percent across the whole grid, then breaks at the opposite edge: at s=0.1s = 0.1 it sits 3.0% above our measurement, a deviation confirmed independently at all five population sizes with a pooled z of −7.33. Haldane's own branching process, taken before the small-ss step that produces 2s2s and solved numerically instead, matches that same measurement to 0.2%, which puts the fault in the diffusion limit rather than in either man's reasoning. Conditional fixation times reproduced the 4N4N neutral prediction to within 5.1% at every size. Seed 20240917.

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±0.00063u = 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.

0.01992fixation prob., N=1000, s=0.01
98.01%of replicates lost the mutation
41.8×better than the neutral rate
11,932,191replicates run in total

The 41.8 deserves a moment. Run the same population with s=0s = 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 2N2N copies of one gene, forever. Every generation the whole adult population dies and is replaced, and each of the 2N2N 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 pp and each mutant copy is worth 1+s1+s as many descendants as a wild-type copy, the pool for the next generation is weighted to

p′=p(1+s)1+spp' = \frac{p(1+s)}{1+sp}

and the number of mutant copies in the next generation is a binomial draw, it+1∼Binomial(2N,p′)i_{t+1} \sim \mathrm{Binomial}(2N, p'). The mutation starts at i=1i = 1. The process runs until ii reaches 0, which we call lost, or 2N2N, 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+s1+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

π=1−e−(1+s)π.\pi = 1 - e^{-(1+s)\pi}.

Haldane then expanded this for small ss and obtained the result everybody quotes:

u≈2s.u \approx 2s.

A one percent advantage gives a two percent chance of establishment. Notice what the formula does not contain. No NN 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)M(p) = sp(1-p) and variance V(p)=p(1−p)/(2N)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 2N2N:

u=1−e−2s1−e−4Ns.u = \frac{1 - e^{-2s}}{1 - e^{-4Ns}}.

Kimura's answer does contain NN, and it contains it only in the combination 4Ns4Ns. 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−4Nsxe^{-4Nsx}, so selection reaches the answer only through how sharply that factor bends across the unit interval. Take 4Ns4Ns large and the denominator goes to 1, leaving u→1−e−2su \to 1 - e^{-2s}, which for small ss is Haldane's 2s2s. Take s→0s \to 0 and L'Hôpital gives u→1/(2N)u \to 1/(2N), the neutral answer, which the discrete chain also gives exactly by a symmetry argument: every one of the 2N2N 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)1/(2N).

So 4Ns4Ns 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

Meeting 1
First version simulated one replicate at a time in a Python loop. Estimated finish for the full grid: about nine hours. Nobody wanted to wait nine hours.

Meeting 1, later
Rewrote it so every replicate in a cell steps forward together. One call to rng.binomial(2N, p) on a whole array. The array shrinks as replicates hit 0 or 2N2N. Late generations cost almost nothing. Whole grid now runs in 66 seconds.

Meeting 2
Argument about how many replicates. 20,000 per cell was the plan. Someone pointed out that for N=10,000N = 10{,}000 with s=0s = 0 the answer is 5×10−55 \times 10^{-5}, so 20,000 replicates would produce about one fixation and a standard error bigger than the estimate. Fix: compute the replicate count from the predicted answer, n≈3×100(1−u)/un \approx 3 \times 100(1-u)/u. Floored at 50,000, capped at 8,000,000.

Meeting 2
That cap matters. N=10,000,s=0N = 10{,}000, s = 0: 5,999,700 replicates, 323 fixations, 36 seconds. Most expensive cell in the study. It exists only to confirm an answer we already knew.

Meeting 3
First run: three cells at a relative standard error slightly over 10%, which was the target we set ourselves. Those estimates had fluctuated low, so the realised error came out worse than planned. Added a factor of 3 everywhere. Restating the target would have been the cheaper fix and the dishonest one. Worst cell is now 5.93%.

Meeting 3
Kimura off by about 3% at s=0.1s = 0.1. Same direction at all five population sizes. Spent most of the meeting hunting the bug. No bug. See §7.

Meeting 4
Seed fixed at 20240917 and never touched again. Every number printed here comes from one run of the file as it now stands.

The Arithmetic

The grid holds five values of NN and eight values of ss. 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 ss on log axes for three of the five population sizes, with Haldane's 2s2s as the dashed straight line and Kimura's curves drawn underneath the points.

Log-log plot of measured fixation probability against selection coefficient for N equal to 50, 1000 and 10000, with the Haldane 2s line and Kimura curves 0.001 0.002 0.005 0.01 0.02 0.05 0.1 selection coefficient s 0.002 0.005 0.01 0.02 0.05 0.1 0.2 fixation probability u Kimura, N = 50 Kimura, N = 10000 Haldane 2s Fixation probability against selection coefficient
  • N = 50
  • N = 1000
  • N = 10000
Figure 1. Measured fixation probability of a single new mutant copy against the selection coefficient, for three of the five population sizes in the grid. Points are club Monte Carlo estimates, bars are 95% binomial confidence intervals. The dashed line is Haldane's 2s2s. The solid curves are Kimura's diffusion result for N=50N = 50 (upper) and N=10,000N = 10{,}000 (lower, essentially on top of the dashed line). At N=50N = 50 the points flatten out on the left because drift, not selection, is setting the outcome: the mutation fixes at roughly the neutral rate 1/(2N)=0.011/(2N) = 0.01 no matter how good it is.

The left-hand end of the N=50N = 50 series carries the whole story. At s=0.001s = 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.

Ns4Nsreplicatesfixationsu measuredSE2svs 2sKimuravs Kimura
N = 50
500050,0004750.0095000.000434n/an/a0.010000−5.0%
500.0010.250,0005880.0117600.0004820.002000+488.0%0.011022+6.7%
500.0020.450,0006030.0120600.0004880.004000+201.5%0.012109−0.4%
500.0051.050,0008010.0160200.0005610.010000+60.2%0.015741+1.8%
500.012.050,0001,1280.0225600.0006640.020000+12.8%0.022901−1.5%
500.024.050,0002,0070.0401400.0008780.040000+0.4%0.039942+0.5%
500.051050,0004,7380.0947600.0013100.100000−5.2%0.095167−0.4%
500.12050,0008,7330.1746600.0016980.200000−12.7%0.181269−3.6%
N = 200
20000119,7002910.0024310.000142n/an/a0.002500−2.8%
2000.0010.882,3842920.0035440.0002070.002000+77.2%0.003628−2.3%
2000.0021.659,6782830.0047420.0002810.004000+18.6%0.005002−5.2%
2000.0054.050,0004780.0095600.0004350.010000−4.4%0.010136−5.7%
2000.018.050,0009880.0197600.0006220.020000−1.2%0.019808−0.2%
2000.021650,0001,9330.0386600.0008620.040000−3.4%0.039211−1.4%
2000.054050,0004,6970.0939400.0013050.100000−6.1%0.095163−1.3%
2000.18050,0008,8180.1763600.0017040.200000−11.8%0.181269−2.7%
N = 1000
100000599,7002860.0004770.000028n/an/a0.000500−4.6%
10000.0014.0147,1002850.0019370.0001150.002000−3.1%0.002035−4.8%
10000.0028.074,8252990.0039960.0002310.004000−0.1%0.003993+0.1%
10000.0052050,0004620.0092400.0004280.010000−7.6%0.009950−7.1%
10000.014050,0009960.0199200.0006250.020000−0.4%0.019801+0.6%
10000.028050,0001,9880.0397600.0008740.040000−0.6%0.039211+1.4%
10000.0520050,0004,7350.0947000.0013090.100000−5.3%0.095163−0.5%
10000.140050,0008,8180.1763600.0017040.200000−11.8%0.181269−2.7%
N = 5000
5000002,999,7003160.0001050.000006n/an/a0.000100+5.3%
50000.00120149,8512870.0019150.0001130.002000−4.2%0.001998−4.1%
50000.0024074,8513270.0043690.0002410.004000+9.2%0.003992+9.4%
50000.00510050,0005020.0100400.0004460.010000+0.4%0.009950+0.9%
50000.0120050,0001,0070.0201400.0006280.020000+0.7%0.019801+1.7%
50000.0240050,0001,9730.0394600.0008710.040000−1.3%0.039211+0.6%
50000.05100050,0004,7130.0942600.0013070.100000−5.7%0.095163−0.9%
50000.1200050,0008,7830.1756600.0017020.200000−12.2%0.181269−3.1%
N = 10000
10000005,999,7003230.0000540.000003n/an/a0.000050+7.7%
100000.00140149,8513150.0021020.0001180.002000+5.1%0.001998+5.2%
100000.0028074,8513130.0041820.0002360.004000+4.5%0.003992+4.8%
100000.00520050,0005170.0103400.0004520.010000+3.4%0.009950+3.9%
100000.0140050,0001,0020.0200400.0006270.020000+0.2%0.019801+1.2%
100000.0280050,0002,0260.0405200.0008820.040000+1.3%0.039211+3.3%
100000.05200050,0004,7790.0955800.0013150.100000−4.4%0.095163+0.4%
100000.1400050,0008,7700.1754000.0017010.200000−12.3%0.181269−3.2%

The neutral rows are the code check. At s=0s = 0 the answer is 1/(2N)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 NN. So wherever NN 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 2s2s against 4Ns4Ns. Were Haldane right, every point would sit on the horizontal line at 1. The 35 cells with s>0s > 0 span four orders of magnitude in 4Ns4Ns, and they land almost exactly on the curve 1/(1−e−4Ns)1/(1 - e^{-4Ns}), which is what Kimura's formula reduces to when ss is small enough that 1−e−2s≈2s1 - e^{-2s} \approx 2s.

Log-log plot of measured fixation probability divided by 2s against 4Ns for all 35 positive-s grid cells, showing collapse onto the curve one over one minus exp minus 4Ns 0.2 1 10 100 1000 4Ns (selection strength relative to drift) 1 2 3 5 7 measured u divided by 2s 1 / (1 - e^-4Ns) 2s exactly right s = 0.05 and 0.1 sit low 35 grid cells, bars = 95% CI All 35 cells collapse onto one curve in 4Ns, except where s is large
  • s ≤ 0.02
  • s = 0.05, 0.1
Figure 2. Every positive-ss cell in the grid, plotted as the measured fixation probability divided by Haldane's 2s2s, against 4Ns4Ns. The solid curve is 1/(1−e−4Ns)1/(1-e^{-4Ns}), the small-ss limit of Kimura's formula. Filled circles are cells with s≤0.02s \le 0.02; open squares are s=0.05s = 0.05 and s=0.1s = 0.1, which sit systematically below the line for a reason that has nothing to do with drift (§7). Cells at very different NN and ss that share a value of 4Ns4Ns land on top of each other, which is the practical meaning of the statement that 4Ns4Ns is the parameter that matters.

The collapse is the point. N=1000N = 1000 with s=0.001s = 0.001, N=200N = 200 with s=0.005s = 0.005 and N=50N = 50 with s=0.02s = 0.02 all share 4Ns=44Ns = 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 2s2s across our cells was:

4Ns bandcellsmean |error| of 2smean |error| of Kimura
below 13255.6%3.1%
1 to 4330.5%2.8%
4 to 2072.5%1.9%
20 to 100106.2%4.0%
100 and above124.9%1.9%

Where does Haldane break down? Below about 4Ns=44Ns = 4, and the breakdown is fast. At 4Ns=24Ns = 2 the error is 12.8%. At 4Ns=14Ns = 1, 60.2%. At 4Ns=0.24Ns = 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 4Ns4Ns 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.1s = 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 5\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.1s = 0.1 is real. Nothing else in the grid deviates at all.

−7.33pooled z, s = 0.1
+0.89largest pooled z elsewhere
−3.0%Kimura's error at s = 0.1
−0.2%branching process error

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→0s \to 0 and N→∞N \to \infty with 4Ns4Ns held fixed. It throws away terms of order s2s^2. At s=0.001s = 0.001 those terms are invisible. At s=0.1s = 0.1 they are worth about three percent, which is exactly what we see.

The test sends you back to Haldane, past 2s2s, which is his answer after the small-ss step, to the branching process he started from, solved without that step. For a lineage whose copy number is Poisson with mean 1+s1+s each generation, the escape probability π\pi solves π=1−e−(1+s)π\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≥1004Ns \ge 100, the region where finite NN has stopped mattering and all three ought to apply.

Fixation probability against selection coefficient at large 4Ns, comparing 2s, one minus exp minus 2s, and the branching process fixed point against pooled club measurements 0.00 0.02 0.04 0.06 0.08 0.10 selection coefficient s 0.00 0.05 0.10 0.15 0.20 fixation probability u 2s 1 - e^-2s branching pi club measurements At s = 0.1 the club number lands on Haldane's own branching process
Figure 3. The large-NN end of the grid, where 4Ns≥1004Ns \ge 100 and population size has stopped affecting the answer. Three predictions: Haldane's published 2s2s (dashed), Kimura's large-NN limit 1−e−2s1 - e^{-2s} (solid), and the escape probability of Haldane's own branching process π=1−e−(1+s)π\pi = 1 - e^{-(1+s)\pi} (dotted), all three in the same grey. Points pool all cells in the grid at that ss with 4Ns≥1004Ns \ge 100; bars are 95% intervals and are narrower than the symbols at the left. At s=0.1s = 0.1 the pooled measurement is 0.175807 ± 0.000983, which misses 2s2s by 12.1%, misses Kimura by 3.0% and misses the branching process by 0.2%.

Three numbers, one measurement. At s=0.1s = 0.1: 2s=0.22s = 0.2 is 12.1% too high, 1−e−2s=0.1812691 - e^{-2s} = 0.181269 is 3.0% too high, and π=0.176134\pi = 0.176134 is 0.2% off. At s=0.05s = 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.1s = 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.1s = 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 kk of them form a valid smaller experiment, which makes the trace a real record of convergence rather than a redrawing of the final answer.

Running Monte Carlo estimate divided by the Kimura prediction against number of replicates for three grid cells, with 95 percent confidence bands 1k 10k 100k 600k replicates used, k 0.0 0.5 1.0 1.5 running estimate / Kimura first fixation by k = 3795 Three cells converging on the theoretical value
  • N = 10000, s = 0.002
  • N = 1000, s = 0.01
  • N = 1000, s = 0
Figure 4. Convergence. The running estimate after the first kk replicates, divided by the theoretical value for that cell, so all three fit on one axis. Shaded bands are 95% confidence intervals at that kk. The neutral cell (N=1000N = 1000, s=0s = 0) is the slow one: it recorded no fixation at all in its first 2,530 replicates, and its running estimate was still wandering between 0.88 and 1.06 times the true value after a quarter of a million. The s=0.01s = 0.01 cell was inside 2% of the answer by roughly 13,000 replicates and stayed there.

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≈300(1−u)/un \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)1/(2N) being exact. Figure 5 shows all five against the line.

Log-log plot of measured neutral fixation probability against number of gene copies, with 95 percent Wilson intervals, lying on the line one over two N 100 400 2000 10000 20000 gene copies in the population, 2N 0.0001 0.001 0.01 neutral fixation probability 475 291 286 316 323 the line u = 1/(2N) small number above each point = fixations counted Neutral case: 9.77 million replicates against an exact answer
Figure 5. The neutral check. Measured fixation probability of a single copy with s=0s = 0, against the number of gene copies 2N2N, on log axes. The line is the exact answer 1/(2N)1/(2N); it is not fitted to anything. Bars are 95% Wilson score intervals, which behave properly when the fixation count is small. The number above each point is how many fixations that estimate rests on: 475, 291, 286, 316 and 323, out of 9,768,800 neutral replicates in total. Largest z score across the five: 1.33.

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 4N4N generations to fix, which for a population of ten thousand is a stretch of time nobody will ever watch. Our measurements, in generations:

Ns4Nsfixationsmean timeSEmedian10th–90th pctbenchmarkratio
5000475197.94.717792–329200 (4N)0.99
20000291759.523.4646360–1320800 (4N)0.95
1000002864038.7113.036582008–63824000 (4N)1.01
50000031619668.6580.8173529896–3257020000 (4N)0.98
100000032339786.21212.13500818401–6904240000 (4N)0.99
selected, N = 10000
100000.001403158468.6111.981836245–11150198070.43
100000.002803134927.152.848253782–612899030.50
100000.0052005172375.316.523361934–284239610.60
100000.014001,0021319.76.012941095–156919810.67
100000.028002,026727.42.0716623–8469900.73
100000.0520004,779330.70.5326289–3783960.83
100000.140008,770182.50.2180161–2071980.92

The neutral rows land within 5.1% of 4N4N at every size. The worst of the five, N=200N = 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)(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 4Ns4Ns 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 2s2s" holds when 4Ns4Ns is above about 4. At 4Ns=0.24Ns = 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.1s = 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 NN 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 4Ns4Ns is correct. Getting NN 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.01s = 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.

4Ns4Ns 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 2s2s 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 NeN_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+s1+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 4Ns4Ns replaced by 2Ns2Ns [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 2s2s 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 4Ns4Ns, so a changing NN 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-NN 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 2s2s, 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,000N = 10{,}000, s=0s = 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

Reproduce these numbers

From the analysis folder, run:

python drift-versus-selection.py > drift-versus-selection-output.txt
Seed
20240917
Analysis code
drift-versus-selection.py46 kB