Why Sex Exists: Watching Mutations Pile Up in a Population That Cannot Recombine

The Click

A ratchet turns one way. That is the whole of it, and the whole of the problem.

Consider a population where nobody recombines. Every individual carries some number of deleterious mutations, and you can sort the population into classes by that number: the ones carrying none, the ones carrying one, and so on up. Mutation moves individuals up. Selection pushes against the move, because carrying more mutations means leaving fewer offspring. Between them they settle into a balance, and the balance has a bottom rung, the class carrying the fewest mutations in the whole population.

Now lose that rung. Not to selection, which would never do it, but to plain bad luck: in some generation, by chance, none of the individuals in the best class leave offspring. The class is gone. And it cannot come back, because mutation only adds and there is no recombination to assemble a clean genome out of two dirty ones. The bottom rung is now one mutation higher than it was, permanently. The ratchet has clicked.

Muller wrote this down in 1932 and named it in 1964 [1][2]. Felsenstein gave it the name it carries now and showed what it implies about why recombination exists at all [3]. The question we took to the club table was narrow and, we thought, tidy: how fast does the thing actually click, and how much sex does it take to stop it?

We got both answers. We also got a third one that nobody asked for, which is that the formula everybody quotes for the click rate is wrong by fifty-five orders of magnitude in the worst place we looked, and wrong by a factor of twenty in the opposite direction somewhere else. That took over the study. This article is mostly about the wreckage.

What Is Actually In The Machine

A haploid Wright-Fisher population of constant size \(N\). One integer per individual, \(k\), the number of deleterious mutations it carries. Fitness is multiplicative, \(w(k) = (1-s)^k\). A generation is two steps. First reproduction: \(N\) offspring drawn with replacement from the parents with probability proportional to fitness. Then mutation: each offspring gains a further \(\mathrm{Poisson}(U)\) new deleterious mutations, with no back mutation and no beneficial mutation ever.

Because an individual's entire state is the integer \(k\), the population can be carried as a histogram over load classes rather than as a list of individuals. For the asexual model that representation is exact and not an approximation, and it is the reason 200 replicate populations across every asexual cell ran in under a minute on a school laptop. One generation in histogram form is: weight the classes by fitness, normalise, convolve with the Poisson mutation kernel, take one multinomial draw of \(N\) offspring.

Recombination cannot be done honestly that way. Two parents carrying \(k_1\) and \(k_2\) mutations produce a recombinant whose load depends on how many of those mutations sit at the same sites, and a load-only model has thrown that information away. So every run involving recombination carries real genomes: 1024 biallelic sites per individual, packed into 64-bit words. A fraction \(R\) of offspring are made by free recombination, each site taken independently from one of two fitness-weighted parents; the remaining \(1-R\) are clones of a single fitness-weighted parent. \(R\) is the rate of sex.

What is not in it

Haploid, so no dominance and no diploid masking. Constant \(N\), so no crashes and no bottlenecks. Non-overlapping generations. One panmictic pool with no geography. Every mutation has exactly the same effect \(s\), which means there is no distribution of fitness effects at all. Fitness is multiplicative, so there is no epistasis of any kind. \(N\) is both the census size and the effective size. And in the recombination runs, recombination is free, every site assorting independently, which is the strongest recombination there is.

Section 11 goes through which of those would have changed the answer and by how much. Two of them would have changed it a great deal.

The Arithmetic

Start with the balance. At deterministic mutation-selection balance under multiplicative fitness, the distribution of mutation counts is Poisson with mean \(\lambda = U/s\). Haigh proved it in 1978 [4]. So the expected size of the best class, the one carrying zero mutations, is

\(n_0 = N e^{-U/s}\).

Put numbers in. Take \(N = 500\), \(U = 0.20\), \(s = 0.05\). Then \(\lambda = 0.20/0.05 = 4\), and \(e^{-4} = 0.018316\), and \(n_0 = 500 \times 0.018316 = 9.158\). Nine individuals out of five hundred are carrying no mutations at all. That is the rung.

Now the classical rate argument, which is one line long. If the best class sits at its equilibrium size \(n_0\), the number of its offspring next generation is \(\mathrm{Binomial}(N, n_0/N)\), which for \(n_0 \ll N\) is \(\mathrm{Poisson}(n_0)\). The class dies out this generation if that draw comes up zero. So

click rate per generation \(\approx e^{-n_0}\).

For our cell that gives \(e^{-9.158} = 1.05 \times 10^{-4}\), one click per 9,488 generations. Our measurement was 0.01909 ± 0.00030, one click per 52 generations. The prediction is too slow by a factor of 181.

Do the same at \(N = 2000\), \(U = 0.10\), \(s = 0.05\). Then \(\lambda = 2\), \(n_0 = 2000 \times e^{-2} = 270.7\), and \(e^{-270.7} = 2.8 \times 10^{-118}\). A click once every \(10^{117}\) generations, which is to say never, in any universe with a finite age. We measured \(6.00 \times 10^{-6} \pm 3.45 \times 10^{-6}\) clicks per generation. One click per 167,000 generations. Rare, but not the same kind of rare at all.

Both arithmetic lines are correct. Both predictions are wrong. The rest of the study is about which step in that one-line argument fails, and the answer turns out to be the first one, the innocuous-looking assumption that the best class is sitting at its equilibrium size when it dies.

Working Notes From The Club Table

The Formula Everyone Quotes, Autopsied

Figure 1 is the whole problem in one panel. The horizontal axis is \(n_0\), the predicted size of the best class. The dashed line is \(e^{-n_0}\). The points are what we measured, at 200 replicate populations each, with error bars that are mostly smaller than the markers.

The dashed line leaves the page near \(n_0 = 13\). The measurements do not. At \(n_0 = 30\) the classical estimate is \(9.4\times10^{-14}\) clicks per generation and the measurement is \(0.00692 \pm 0.00011\), which is a ratio of \(7.4\times10^{10}\). At \(n_0 = 135.3\) the estimate is \(1.7\times10^{-59}\) and the measurement is \(2.34\times10^{-4} \pm 2.18\times10^{-5}\), one click per 4,273 generations, sitting eleven standard errors above zero and \(1.4\times10^{55}\) times the prediction. That is the worst cell in the study and it is worth saying plainly what it means: the formula says this population is safe forever, and it is not. Gordo and Charlesworth said as much in 2000, warning that the ratchet can turn at a biologically significant rate even when the deterministic mutation-free class is well above a hundred individuals [9]. Our \(n_0 = 135.3\) cell is a direct instance.

And the formula fails in the other direction too. Where \(n_0\) is small, below about 1, the classical estimate is too fast. At \(N = 1000\), \(U = 0.20\), \(s = 0.02\), \(n_0\) works out to 0.0454, and \(e^{-0.0454} = 0.9556\), which says the ratchet should click almost every generation. We measured 0.04651 ± 0.00048. That is a twentieth of the prediction, and it sits 1873 standard errors below it. Our Monte Carlo error is nowhere near large enough to cover a gap like that.

The "n₀ observed" column is the autopsy finding. It is the average size of the least-loaded class as the runs actually saw it, measured in the moving frame that rides up with the minimum. Compare it to the predicted column. Where the prediction says 0.00454 individuals, the runs contain 6.3. Where it says 135.3, the runs contain 130.6, which is close. The deterministic equilibrium is a good description of a population whose ratchet is quiet and a fantasy about a population whose ratchet is turning, because a turning ratchet never gets to sit at equilibrium. That is the step that breaks.

Haigh's equilibrium itself is fine. Restricting to replicate-generations in which the ratchet has not clicked at all, where a deterministic equilibrium is at least defined, the mean load came to 1.0288 times \(\lambda\), the variance to 1.0249 times \(\lambda\), and the size of class zero to 0.9533 ± 0.0106 times \(N e^{-\lambda}\) over fourteen cells. Three predictions, all within about 5%. The theory of the balance survives. The theory of the click does not, and it was built on the balance.

The Number That Was Missing

If \(e^{-n_0}\) fails, the first question is whether any function of \(n_0\) could work. Maybe the exponential is too aggressive and the right answer is a power law in the same variable. That is a testable statement, and we tested it the cheap way first, by sweeping \(N\) at fixed \(s = 0.05\) and \(U = 0.20\) so that \(n_0\) walks from 0.49 up to 30 while nothing else moves.

Fit both forms over the thirteen sweep cells with \(n_0 \ge 1\). The exponential fit gives \(\ln(\text{rate}) = -3.1031 - 0.0707\,n_0\) with an rms residual of 0.1957. The power law gives \(\ln(\text{rate}) = -2.6635 - 0.6293\ln n_0\) with an rms residual of 0.0819. The power law wins by more than a factor of two on residual. Over this window the ratchet slows roughly as \(n_0^{-0.63}\), so the interclick time grows as \(n_0^{0.63}\), which is close to the linear growth in \(n_0\) that Stephan, Chao and Smale derived for intermediate \(N\) and \(s\) [6] and nothing like exponential. Note also the exponential decay constant we measure: 0.0707. The classical form demands 1. It is fourteen times too weak.

But that was the cheap test, and it was the wrong question. Here is the expensive one.

Build three cells that share the same \(n_0 = 5\) and the same \(\lambda = 4\), by choosing \(N\), \(U\) and \(s\) together, with \(s\) varying five-fold across them. If the rate is a function of \(n_0\) alone, they must agree.

They do not. At \(s = 0.02\) the rate is \(0.01788 \pm 0.00030\). At \(s = 0.10\) it is \(0.03544 \pm 0.00054\). The difference is \(0.01756 \pm 0.00062\), which is 28.5 standard errors from zero. Raising the selection coefficient at fixed \(n_0\) makes the ratchet turn faster, by very nearly a factor of two. Stephan, Chao and Smale predicted that sign in 1993, writing that for intermediate \(N\) and \(s\) at fixed \(n_0\), increasing \(s\) accelerates the ratchet [6]. We did not know that when we built the cells. We found it afterwards, which is a worse way to do science and an honest description of what happened.

So the rate needs \(s\) as a separate ingredient. Jain put the interclick time into a scaling form set by the size of the least-loaded class together with the selection coefficient [11]. Neher and Shraiman are blunter, writing that the ratchet rate "cannot depend on \(n_0\) alone, but must depend on \(n_0 s\) instead" [13]. That is a prediction with a sharp edge, and we designed eighteen cells to break it: five values of \(s \times n_0\) crossed with three values of \(s\), all at \(\lambda = 4\), 200 fresh replicate populations each. If the claim holds, then rate divided by \(s\) should be a function of \(s \times n_0\) alone, and cells matched on that product should agree even when \(s\) and \(n_0\) each differ five-fold.

Figure 2 is the answer, and it is the cleanest result in the study. Panel A matches on \(n_0\) and gets three curves. Panel B matches on \(s \times n_0\) and gets one. The numbers behind the picture: matched on \(s \times n_0\), the spread in rate over \(s\) runs from 0.3% to 6.9%, median 3.7%, which for the tightest of them is 0.1 standard errors. Matched on \(n_0\) alone, the spread runs from 24.1% to 78.4%, median 67.1%, the worst at 41.8 standard errors. The collapse variable is tighter by a factor of 18.1 on the medians.

This is the one place in the study where a modern prediction was put at risk and came through. It also explains why the classical formula could never have been patched. No constant in front of \(e^{-n_0}\), no adjusted exponent, nothing that is a function of \(n_0\) by itself, can reproduce a rate that moves by a factor of two when \(s\) moves and \(n_0\) does not.

Watching An Estimate Settle

Every number in this article is a mean over 200 independently seeded replicate populations, so every number has a standard error, and the honest thing is to show those errors shrinking rather than assert them. Figure 3 does that for three cells: it plots the running mean click rate against the number of replicates folded in, with a band at plus and minus one standard error of the running mean.

The shape is the one you want. A single replicate is worthless: the \(N = 100\) cell opens at 0.0600 from one population, wanders to 0.0440 by five, and is inside half a percent of its final value by thirty. The band closes as \(1/\sqrt{k}\) and nothing surprising happens at the end. Final values and relative standard errors: 0.04765 ± 0.00063 (1.32%), 0.01909 ± 0.00030 (1.57%), 0.01135 ± 0.00019 (1.67%).

Those percentages are the reason we can say anything at all about a discrepancy of 1873 standard errors, and they are also the reason we cannot say much about the \(N = 2000\), \(s = 0.1\), \(U = 0.1\) cell, which clicked zero times in 500,000 replicate-generations. Zero clicks is a measurement with an error bar of zero and an interpretation of "slower than about \(6\times10^{-6}\) per generation, the rule-of-three bound on seeing no events in 500,000 replicate-generations", which is not the same as "never" and should not be written down as if it were. We report it as zero in the table and mean by it only that we did not see one.

Cells were also run with an early stop: a cell halts once it has banked 2000 clicks, so a fast cell measures for 400 generations and a slow one for the full 2500. That keeps the relative error roughly constant across the grid instead of spending all the compute where the answer was already obvious.

Two Thousand Generations Of Downhill

The click rate is an abstraction. Here is what it does to a population.

Every click costs the population one mutation's worth of fitness at the bottom of the distribution, and because the whole distribution rides above the minimum, it costs the mean the same factor. So if the load distribution is otherwise stationary in the moving frame, the bookkeeping is forced:

\(\dfrac{d\ln \bar{W}}{dt} = (\text{click rate}) \times \ln(1-s)\).

The lines are straight and they keep going. That is the grim part. There is no floor in this model, no compensatory mutation, no point at which the population has lost enough that selection starts winning again; mean fitness is an exponential decay with a rate set by how often the ratchet turns. The worst regime we ran, \(N = 2000\), \(U = 0.80\), \(s = 0.10\), lost a factor of 5.583×10⁻⁶ in mean fitness over 1500 generations while its mean load climbed from 29.2 to 144.3 mutations per individual. Lynch and colleagues called the endpoint of this a mutational meltdown, where declining fitness shrinks the population, which speeds the ratchet, which shrinks it further [7]. Our \(N\) is pinned, so we never see the feedback. We only see the first half of it, running forever in a straight line.

The identity holds tightly. Across six regimes spanning a ten-fold range in decay rate, the measured slope divided by (click rate) \(\times \ln(1-s)\) averaged 0.9961 with a standard deviation of 0.0031. The worst single regime was 0.9920. That is not a validation of biology, it is a validation of the code: two quantities computed by completely different routes through the same simulation, one by counting integer steps in the minimum and one by fitting a line through the log of a mean, agree to four parts in a thousand. If they had not, something would have been wrong with the accounting.

The flat line at the top of Figure 4 is the obligate sexual control from the recombination runs, which is where the story finally turns.

How Much Sex Is Enough

The base cell for this half is deliberately nasty: \(N = 200\), \(U = 0.30\), \(s = 0.05\), so \(\lambda = 6\) and \(n_0 = 0.496\), meaning that at equilibrium there is on average half an individual carrying no mutations. Genomes are explicit, 1024 sites, 200 replicates, 800 generations after a 150-generation burn-in. Asexually, that population's minimum load drifts upward at 0.06158 ± 0.00046 per generation, one click every 16.2 generations.

Then we turn on sex, a fraction \(R\) of offspring at a time.

Read the third column before the sixth, because it is the joke at the heart of the result. Recombination makes the minimum load go up more often, not less. At \(R = 0\) the minimum rises 0.0616 times per generation; at \(R = 1\) it rises 0.2688 times per generation, four times as often. Sex is not preventing clicks. It is un-clicking them. The fourth column from the right counts downward steps, generations in which a recombinant turned up carrying fewer mutations than anything alive the generation before, and it goes from exactly zero to 42,985. The ratchet still has a pawl. Recombination has simply made it a two-way pawl, and the net drift is the difference between two large and nearly equal numbers.

The thresholds, interpolated in \(\log R\) between the bracketing grid points:

At \(R = 0.20\) the net drift is \(1.62\times10^{-4} \pm 1.08\times10^{-4}\), 1.5 standard errors from zero. At \(R = 1\) it is \(1.06\times10^{-4} \pm 8.21\times10^{-5}\), 1.3 standard errors from zero. Within our resolution the ratchet is stopped rather than merely slowed, and the mean load stops climbing and sits at 6.29, which is \(\lambda = 6\) plus the finite-sites correction. Panel B of Figure 5 shows what that buys: the asexual line falls through 2.5 natural log units of mean fitness in 800 generations while the obligate sexual line is flat to the eye and flat to the statistics.

The headline is that a small amount of sex does most of the work. Fifteen percent of offspring, roughly one in seven, buys a hundred-fold reduction. That is much cheaper than the answer we expected going in, and it is the strongest argument in the study for why facultative sex is common in nature: you do not need much.

The Strongest Objection We Could Make

Here is the objection, stated as strongly as we can put it.

You have a simulation whose answer disagrees with a textbook formula by fifty-five orders of magnitude. The overwhelmingly likely explanation is that your simulation is wrong. Simulations written by school students are wrong all the time; closed-form results published in Theoretical Population Biology are wrong considerably less often. You have a window of load classes that you slide upward, which is exactly the kind of index arithmetic that goes wrong silently. You have a multinomial draw whose probabilities you renormalise twice. You have an early-stopping rule that halts different cells after different numbers of generations, which could easily correlate stopping time with the thing you are measuring. Any one of those could produce a spuriously fast ratchet. Why should anybody believe you rather than Haigh?

We take that seriously, and the answer is not a rebuttal, it is four pieces of evidence.

The null. Set \(U = 0\). Now there is no mutation, so there is no way for the minimum to move, and any bug in the window-sliding or the multinomial normalisation has a clean shot at producing a spurious click. Six cells, 200 replicates each, 1000 generations each: 1,200,000 replicate-generations, zero clicks, and mean fitness printing as 1.000000000000000 to fifteen decimal places. The explicit-genome engine passes the same test separately. A window that slid wrongly would not survive that.

Three implementations. The fast histogram engine, a deliberately dumb individual-by-individual version written from the model description rather than from the fast code, and the packed-genome engine run at \(R = 0\). Three people wrote three different things. At \(N = 100\), \(U = 0.20\), \(s = 0.05\) they return 0.04757 ± 0.00061, 0.04540 ± 0.00133 and 0.04422 ± 0.00085. The histogram engine and the dumb one differ by 1.5 standard errors. The packed-genome engine came in 3.2 standard errors low, which we chased and found: it loses 1.84% of its mutations to sites that are already mutated, so its realised \(U\) is 0.1963 rather than 0.20. Run the histogram engine at that realised \(U\) and the gap falls to 1.9 standard errors. At \(N = 200\), \(U = 0.30\) the same correction brings a 1.8 standard error gap down to 0.1.

The conservation check. The fitness identity in §8 ties the click rate to the slope of log mean fitness. Those two quantities are computed by entirely separate paths: one counts integer increments of \(k_{\min}\), the other fits a line to \(\ln \bar{W}\). They agree to 0.9961 ± 0.0031. A bug that inflated the click count would break that identity, because it would not simultaneously steepen the fitness slope by the matching amount.

Haigh still passes. If the code were broken in a way that fabricated clicks, the load distribution would be wrong too. It is not. Mean load, variance of load and size of the unloaded class all land within about 5% of Haigh's predictions in the cells where those predictions are defined [4]. The code reproduces the theory exactly where the theory applies and departs from it exactly where the theory's assumption (a class sitting at equilibrium) is violated.

The remaining part of the objection stands, and we want to be clear about it: the early-stopping rule does correlate measurement length with click rate, because that is what it is for. We checked the direction it could bias things by running the convergence cells to a fixed 400, 600 and 900 generations and comparing against the same cells stopped on clicks, and the rates agree, but we did not run the full grid both ways. If somebody re-runs this with fixed-length cells throughout and gets a different number, that is a real finding and we would like to hear it.

The honest summary is that we did not set out to find \(e^{-n_0}\) failing and we spent two meetings trying to make it stop failing. Stephan, Chao and Smale said in 1993 that the approximations available then fit simulation well only in restricted regimes [6]; Gessler documented the constraints finite size imposes on exactly this quantity [8]; Gordo and Charlesworth built a diffusion approximation specifically because the one-generation argument was inadequate, and stated its validity condition as \(N e^{-U/s} \gg 1\) [9][10]. We are not overturning anything. We are measuring the size of a hole that the literature has known about for thirty years and that the teaching version of the formula does not mention.

Where A Different Model Would Have Given A Different Answer

Four modelling choices carry most of the weight, and two of them would move the headline numbers a long way.

Equal effects. Every mutation in this study costs exactly \(s\). Real genomes have a distribution of fitness effects spanning many orders of magnitude, and that changes the character of the question rather than its answer. With a distribution, the ratchet runs in the nearly neutral tail while strongly selected sites hold still, and "the" click rate stops being a single number. Our §6 result says the rate is governed by \(s \times n_0\); with a spread of \(s\) values there is no single \(s\) to put in it. We would expect the effective ratchet to be faster than our equal-effects model at the same mean \(s\), because the slow-selection tail is where \(n_0\) is smallest. This is the choice most likely to have changed our answer, and we did not test it.

Free recombination. Our \(R = 0.152\) threshold is the amount of sex needed when every site assorts independently, which is the most efficient recombination physically possible. A real chromosome has linkage, and linked sites recombine only across the map distance between them. So 15.2% is a lower bound on the rate of sex a real genome would need, and the true requirement for a genome with realistic map lengths would be higher, possibly much higher. Anyone quoting our threshold as "one in seven does it" without that caveat is misquoting us.

No beneficial or compensatory mutation. Our populations have no road back, which is why the fitness lines in Figure 4 run straight down forever. Real asexual lineages have compensatory mutation, gene conversion, occasional sex, and selection on population size. Any one of those can hold a ratchet that this model runs. Howe and Denver documented compensatory change doing exactly that job in the non-recombining mitochondrial genome of Caenorhabditis briggsae [12]. Our model can say nothing about when that rescue arrives and when it does not, because we did not put the mechanism in. What our populations get instead is the version with the rescue deleted, which is why they fall forever.

Non-overlapping generations. Everybody in our populations reproduces at once and then dies. Metzger and Eule showed that letting generations overlap changes the distribution of the fittest individuals and therefore the ratchet rate [14]. We did not test the size of that effect and we should not guess at it.

One choice that we checked and that did not matter: the load-class window. We were worried that truncating the distribution at a fixed number of classes above the minimum would quietly discard the tail that matters. The worst probability mass lost off the top of the window anywhere in the grid was \(4.44\times10^{-16}\), which is floating-point noise.

What the study does say, stated narrowly: in a haploid Wright-Fisher population with equal-effect multiplicative deleterious mutations and no back mutation, the rate at which the least-loaded class is lost is governed by \(s \times n_0\) rather than \(n_0\), the one-generation estimate \(e^{-n_0}\) is unusable across our whole grid, mean fitness declines at exactly (click rate) \(\times \ln(1-s)\), and free recombination in about one offspring in seven reduces the drift to a hundredth of the asexual rate. Everything past that sentence is extrapolation, and the extrapolation to any actual organism is somebody else's experiment.

Reproducing This

One file, no arguments, no data files, no configuration. It needs Python 3.12 and numpy 2.0 or later (for numpy.bitwise_count, used by the packed-genome engine); nothing else is imported beyond math, sys and time.

python mullers-ratchet.py > mullers-ratchet-output.txt

Expect under three minutes on a modern laptop. Ours took 169.7 s on Python 3.12.3 and numpy 2.4.2. The section timings printed at the end of the output show where it goes: the 60-cell grid takes 46 s and the recombination sweep takes 91 s, which is more than half the total, because explicit 1024-site genomes cost roughly a hundred times what the load-class histogram costs per generation. The master seed is 20250921, hard-coded at the top of the file, and every parameter cell draws an independent stream from it through numpy's SeedSequence.spawn, so the output is deterministic: we re-ran the file after writing this article and it reproduced the earlier run byte for byte. 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 both engines in your browser, and its default settings reproduce the asexual click rate and the recombination threshold quoted above.

References