# Fluctuations of Fitness Distributions and the Rate of Muller’s Ratchet

^{*}Max-Planck Institute for Developmental Biology, 72070 Tübingen, Germany, and^{†}Kavli Institute for Theoretical Physics and Department of Physics, University of California, Santa Barbara, California 91306

- 1Corresponding author: Max-Planck Institute for Developmental Biology, Spemannstr 35, 72076 Tübingen, Germany. E-mail: richard.neher{at}tuebingen.mpg.de

## Abstract

The accumulation of deleterious mutations is driven by rare fluctuations that lead to the loss of all mutation free individuals, a process known as Muller’s ratchet. Even though Muller’s ratchet is a paradigmatic process in population genetics, a quantitative understanding of its rate is still lacking. The difficulty lies in the nontrivial nature of fluctuations in the fitness distribution, which control the rate of extinction of the fittest genotype. We address this problem using the simple but classic model of mutation selection balance with deleterious mutations all having the same effect on fitness. We show analytically how fluctuations among the fittest individuals propagate to individuals of lower fitness and have dramatically amplified effects on the bulk of the population at a later time. If a reduction in the size of the fittest class reduces the mean fitness only after a delay, selection opposing this reduction is also delayed. This delayed restoring force speeds up Muller’s ratchet. We show how the delayed response can be accounted for using a path-integral formulation of the stochastic dynamics and provide an expression for the rate of the ratchet that is accurate across a broad range of parameters.

BY weeding out deleterious mutations, purifying selection acts to preserve a functional genome. In sufficiently small populations, however, weakly deleterious mutations can by chance fix. This phenomenon, termed Muller’s ratchet (Muller 1964; Felsenstein 1974), is especially important in the absence of recombination and is thought to account for the degeneration of Y chromosomes (Rice 1987) and for the absence of long-lived asexual lineages (Lynch *et al.* 1993).

A click of Muller’s ratchet refers to the loss of the class of individuals with the smallest number of deleterious mutations. To understand the processes responsible for such a click, it is useful to consider a simple model of accumulation of deleterious mutations with identical effect sizes illustrated in Figure 1. Because of mutations, the population spreads out along the fitness axis, which in this model is equivalent to the number of deleterious mutations in a genome. The population can hence be grouped into discrete classes, each characterized by the number of deleterious mutations. Mutation carries individuals from classes with fewer to classes with more mutations, hence shifting the population to the left. This tendency is opposed by selection, which amplifies fit individuals on the right, while decreasing the number of unfit individuals on the left. These opposing trends lead to a steady balance, at least in sufficiently large populations. However, in addition to selection and mutation, the distribution of individuals among fitness classes is affected by fluctuations in the number of offspring produced by individuals of different classes, *i.e.*, by genetic drift. Such fluctuations are stronger (in relative terms) in smaller populations and in particular in classes that carry only a small number of individuals. When mutation rate is high and selection is weak, the class of individuals with the smallest number of mutations (*k* = 0 in Figure 1) contains only few individuals and is therefore susceptible to accidental extinction—an event that corresponds to the “click” of Muller’s ratchet.

Despite the simplicity of the classic model described above, understanding the rate of the ratchet has been a challenge and remains incomplete (Stephan *et al.* 1993; Gessler 1995; Higgs and Woodcock 1995; Gordo and Charlesworth 2000; Stephan and Kim 2002; Etheridge *et al.* 2007; Jain 2008; Waxman and Loewe 2010). Here, we revisit this problem starting with the systematic analysis of fluctuations in the distribution of the population among different fitness classes. We show that fitness classes do not fluctuate independently. Instead, there are collective modes affecting the entire distribution, which relax on different time scales. Having identified these modes, we calculate the fluctuations of the number of individuals in the fittest class and show how these fluctuations affect the mean fitness. Fluctuations in mean fitness feed back on the fittest class with a delay and thereby control the probability of extinction. These insights allow us to arrive at a better approximation to the rate of the ratchet. In particular, we show that the parameter introduced in the earlier work (Haigh 1978; Stephan *et al.* 1993; Gordo and Charlesworth 2000) to parameterize the effective strength of selection in the least loaded class is not a constant but depends on the ratio of the mutation rate and the effect size of mutations. We use the path-integral representation of stochastic processes borrowed from physics (Feynman and Hibbs 1965) to describe the dynamics of the fittest class and arrive at an approximation of the rate of Muller’s ratchet that is accurate across a large parameter range.

Understanding the rate of the ratchet is important, for example, to estimate the number of beneficial mutations required to halt the ratchet and prevent the mutational meltdown of a population (Lynch *et al.* 1993; Goyal *et al.* 2012; Pfaffelhuber *et al.* 2012) (for an in-depth and up-to-date discussion of the importance of deleterious mutations we refer the reader to Charlesworth 2012). Furthermore, fluctuations of fitness distributions are a general phenomenon with profound implications for the dynamics of adaptation and genetic diversity of populations. Below we place our approach into the context of the recent studies of the dynamics of adaptation in populations with extensive nonneutral genetic diversity (Tsimring *et al.* 1996; Rouzine *et al.* 2003; Desai and Fisher 2007; Neher *et al.* 2010). The study of fluctuations in the approximately stationary state of mutation selection balance that we present here is a step toward more general quantitative theory of fitness fluctuations in adapting populations.

## Model and Methods

We assume that mutations happen at rate *u* and that each mutation reduces the growth rate of the genotype by *s* ≪ 1. Within this model, proposed and formalized by Haigh, individuals can be categorized by the number of deleterious mutations they carry. The equation describing the fitness distribution in the population, *i.e.*, what part *n _{k}* of the population carries

*k*deleterious mutations, is given by

*i.e.*, genetic drift, and has the properties of uncorrelated Gaussian white noise with 〈

*η*(

_{k}*t*)

*η*(

_{l}*t*′)〉 =

*δ*(

_{kl}δ*t*−

*t*′). In the infinite population limit, this equation has the well-known steady-state solution

*λ*=

*u*/

*s*. A time-dependent analytic solution of the deterministic model has been described in Etheridge

*et al.*(2007).

Note that we have deviated slightly from the standard model, which assumes that genetic drift amounts to a binomial resampling of the distribution with the current frequencies *N*^{−1}*n _{k}*. This choice would result in off-diagonal correlations between noise terms that stem from the constraint that the total population size is strictly constant. This exact population size constraint is an arbitrary model choice that we have relaxed to simplify the algebra. Instead, we control the population size by a soft constraint that keeps the population constant on average but allows small fluctuations of

*N*. The implementation of this constraint is described explicitly below. We confirmed the equivalence of the two models by simulation.

### Computer simulations

We implemented the model as a computer simulation with discrete generations, where each generation is produced by a weighted resampling of the previous generation. Specifically,*N*_{0}. This specific discretization is chosen because it has exactly the same stationary solution as the continuous-time version above (Haigh 1978). The simulation was implemented in Python using the scientific computing environment SciPy (Oliphant 2007). If the parameter of the Poisson distribution was larger than 10^{4}, a Gaussian approximation to the Poisson distribution was used to avoid integer overflows.

To determine the ratchet rate, the population was initialized with its steady-state expectation ^{4} generations, and then run for further *T* = 10^{8} generations. Over these 10^{8} generations, the number of clicks of the ratchet were recorded and the rate was estimated as clicks per generation.

The source code of the programs, along with a short documentation, is available as supporting information (File S2). In addition, we also provide some of the raw data and the analysis scripts producing the figures as they appear in the manuscript.

### Numerical determination of the most likely path

The central quantities in our path-integral formulation of the rate of Muller’s ratchet are (i) the most likely path to extinction and (ii) the associated minimal action *m* equidistant time points *ρ _{i}* between 0 and

*τ*, where

*x*

_{0}(

*ρ*) = 0. For a given set of

_{m}*x*

_{0}(

*ρ*), a continuous path

_{i}*x*

_{0}(

*ρ*) is generated by linear interpolation. For a given trajectory

*x*

_{0}(

*ρ*), we determine the mean fitness by solving the deterministic equations for

*x*(

_{k}*ρ*),

*k*≥ 1. Note that in this scheme, the only independent variable is the path

*x*

_{0}(

*ρ*), and all other degrees of freedom are slaved to

*x*

_{0}(

*ρ*). From

*x*

_{0}(

*ρ*) and the resulting

*S*

_{λ}({

*x*

_{0}(

*ρ*)}) as defined in Equation 28.

*S*

_{λ}({

*x*

_{0}(

*ρ*)}) is then minimized by changing the values of

*x*

_{0}(

*ρ*), 0 <

_{i}*i*<

*m*, using the simplex minimization algorithm implemented in SciPy (Oliphant 2007). To speed up convergence, the minimization is first done with a small number of pivot points (

*m*= 4), which is increased in steps of 2 to

*m*= 24. The total time

*τ*= 20 (in units of

*s*

^{−1}) was used, which is sufficiently large to make the result independent of

*τ*. The code used for the minimization is provided in File S2.

## Results and Discussion

Fluctuations of the size *n*_{0} of the least loaded class can lead to its extinction. In the absence of beneficial mutations this class is lost forever (Muller 1964), and the resulting accumulation of deleterious mutations could have dramatic evolutionary consequences. Considerable effort has been devoted to understanding this process, and it has been noted that the rate at which the fittest class is lost depends strongly on the average number of individuals in the top class *et al.* 2008). However, a quantitative understanding of the

Here, we present a systematic analysis of the problem by first analyzing how selection stabilizes the population against the destabilizing influences of mutation and genetic drift, and later we use this insight to derive an approximation to the rate of Muller’s ratchet. Before analyzing Equation 1, it is useful to realize that it implies a common unit of time for the time derivative, the mutation rate, and the selection coefficient that is of our choosing (days, months, generations, etc.). We can use this freedom to simplify the equation and reveal what the important parameters are that govern the behavior of the equation. In this case, it is useful to use *s*^{−1} as the unit of time and work with the rescaled time *τ* = *ts*. Furthermore, we formulate the problem in terms of frequencies *x _{k}* =

*N*

^{−1}

*n*rather than numbers of individuals, and obtain

_{k}*λ*=

*u*/

*s*is the dimensionless ratio of mutation rate and selection strength. In other words,

*λ*is the average number of mutations that happen over a time

*s*

^{−1}(our unit of time). Note that

*λ*uniquely specifies the deterministic part of this equation and its steady-state solution

*Ns*)

^{−1}has a simple interpretation as the variance of the stochastic effects accumulated over time

*s*

^{−1}. Other than through a prefactor determining the unit of time, any quantity governed by Equation 3 can depend only on

*λ*and

*Ns*. Hence it is immediately obvious that the ratchet rate cannot depend on

Before turning to the ratchet rate, we analyze in greater detail the interplay of deterministic and stochastic forces in Equation 3. A full time-dependent analytic solution of the deterministic model was found in Etheridge *et al.* (2007). Below, we present an analytic characterization of the stochastic properties of the system in a limit where stochastic perturbations are small.

### Linear stability analysis

In the limit of large populations, the fluctuations of *x _{k}* around the deterministic steady state

*L*. A quick calculation shows that

_{km}*L*has eigenvalues

_{km}*κ*= −

_{i}*i*with

*i*= 0, 1, 2.... The right eigenvector corresponding to

*κ*

_{0}= 0 is simply

*i*> 0 are given by

*k*numbers the coordinate of the vector. This is readily verified by direct substitution (note that

*i*< 0).

The eigenvector *Model and Methods*). All other eigenvalues are negative, which is to say that

The eigenvectors for *i* > 0 have an intuitive interpretation: Eigenvector *i* corresponds to a shift of a fraction of the population by *i* fitness classes downward. Since such a shift reduces mean fitness, the fittest classes start growing and undo the shift. More generally, any small perturbation of the population distribution can be expanded into eigenvectors *a _{j}*(

*τ*) will decay exponentially in time with rate

*j*(remember that the unit of time is

*s*

^{−1}). Since the amplitudes are projections of

*δx*onto the left eigenvectors of

_{k}*L*, we need to know those as well. For

_{mk}*κ*

_{0}= 0, the left eigenvector is simply

*k*>

*i*. With the left and right eigenvectors and the eigenvalue spectrum of the deterministic system on hand, we now reinstantiate the stochastic part of the dynamics:

*x*-dependent noise term is reintroduced later when we turn to Muller’s ratchet. Substituting the representation of

_{k}*η*contributes to every

_{k}*a*and induces correlations between the

_{j}*a*, but each amplitude can be integrated explicitly:

_{j}*δx*

_{0}(

*τ*) and

### Fluctuations of x_{0} and the mean fitness

_{0}

For the variance of the fittest class *τ* = 0 in the above expression) is therefore *λ*, it simplifies to *x*_{1} ≈ *λ* = *u*/*s*. The opposite limit of large *λ* corresponds to a broad fitness distribution where the top class represents only a very small fraction of the entire population. In this limit, the leading behavior of the variance *σ*^{2} is *λ* and compared to simulation results, which agree within measurement error. In our rescaled units, the correlation functions decay over a time of order 1, corresponding to a time of order 1/*s* in real time. More precisely, the decay time (in scaled units) increases with increasing *λ* as log *λ*.

In a similar manner, we can calculate the autocorrelation of the mean fitness*λ* at *τ* = 0. It is hence inversely proportional to the size of the fittest class *x*_{0}. For large *λ*, *x*_{0} represents only a tiny fraction of the population and fluctuations of the mean can be substantial even for very large *N*. This emphasizes the importance of fluctuations of the size of the fittest class for properties of the distribution.

If fluctuations of the mean fitness *x*_{0}, we expect a strong correlation between those fluctuations (Etheridge *et al.* 2007). Furthermore, fluctuations of *x*_{0} should precede fluctuations of the mean. These expectations are confirmed by the analytic result*λ*. The cross correlation *λ*, the peak of the correlation function moves slowly (logarithmically) to larger delays. This result is intuitive, since we expect that fluctuations in the fittest class will propagate to less and less fit classes and that the dynamics of the entire distribution is, at least partly, slaved to the dynamics of the top class.

In all of these three cases, the magnitude of the fluctuations is governed by the parameter *Ns*, while the shape of the correlation function depends on the parameter *λ*. Only the unit in which time is measured has to be compared to the strength of selection directly.

### The rate of Muller’s Ratchet

The ratchet clicks when the size of the fittest class hits 0, and the rate of the ratchet is given by the inverse of the mean time between successive clicks of the ratchet. Depending on the average size *et al.* (2008). Conversely, if *x*_{0} prior to the click, while Figure 3B shows the realized trajectory that ends at *x*_{0} = 0. Prior to extinction, *x*_{0}(*τ*) fluctuates around its equilibrium value and large excursions are rare and short. The final fluctuation that results in the click of the ratchet is zoomed in on in Figure 3C. Compared to the time the trajectory spends near *s* = 0.01 in this example).

In rescaled time, the equation governing the frequency of the top class is*x*_{0}, as well as on the size of the other classes *x _{k}*. For sufficiently large

*λ*,

*x*

_{0}is much smaller than

*x*with

_{k}*k*≥ 1, such that the stochastic force is most important for

*x*

_{0}. The dynamics of

*x*,

_{k}*k*≥ 1, is approximately slaved to the stochastic trajectory of

*x*

_{0}(

*τ*). We can therefore try to find an approximation of

*x*

_{0}(

*τ*) only. The linear stability analysis of the mutation–election balance has taught us that the restoring force exerted by the mean fitness on fluctuations in

*x*

_{0}is delayed with the delay increasing ∝ log

*λ*. The latter observation implies that the restoring force on

*x*

_{0}will depend mainly on the values of

*x*

_{0}some time of order log

*λ*in the past. Such history dependence complicates the analysis, and this delay has been ignored in previous analysis, which assumed that

*x*

_{0}(

*τ*) via

*et al.*1993; Gordo and Charlesworth 2000; Jain 2008). The parameter

*α*was chosen ad hoc between 0.5 and 0.6. This restoring force is akin to a harmonic potential centered around

*P*(

*x*

_{0},

*τ*),

*x*

_{0}by

*z*for simplicity. The fact that the fittest class is lost whenever its size hits 0 corresponds to an absorbing boundary condition for

*P*(

*z*,

*τ*) at

*z*= 0. For such a one-dimensional diffusion problem, the mean first passage time can be computed in closed form (Gardiner 2004) and this formula has been used in Stephan

*et al.*(1993) and Gordo and Charlesworth (2000) to estimate the rate of the ratchet. An accurate analytic approximation to that formula has been presented by Jain (2008). For completeness, we present an alternative derivation of these results that helps interpret the more general results presented below. In the limit of interest,

*z*. We can therefore approximate the distribution as

*P*(

*z*,

*τ*) ≈

*e*

^{−}

*(*

^{γτ}p*z*), where

*γ*is the rate of the ratchet. In this factorization,

*p*(

*z*) is the quasi-steady distribution shown in Figure 3A, while

*γ*is the small rate at which

*P*(

*z*,

*τ*) loses mass due to events like the one shown in Figure 3C. Inserting this ansatz and integrating Equation 18 from

*z*to ∞, we obtain

*γ*, we solve this equation in a regime of small

*z*≫ (

*Ns*)

^{−1}, where the term on the left can be neglected. For the general discussion below, it is useful to solve this equation for a general diffusion constant

*D*(

*z*) (here equal to

*z*/2

*Ns*), force field

*A*(

*z*) (here equal to

*C*

*C*=

*γ*and

*β*is fixed by the boundary condition that

*p*(

*z*) is finite at

*z*= 0. Note that

*γ*=

*p*(0)/2

*Ns*relates the rate of extinction to

*p*(0). To determine the latter we need to match the

*C*is unimportant in this regime (

*e*

^{2}

*≫ 1). Setting*

^{Nsαz}*C*= 0 in Equation 21, we find

*β*in Equation 21 is fixed by the normalization. Since

*p*(

*z*) is concentrated around

*D*(

*z*) term, scaled so that it equals 1 in the vicinity of

*p*(

*z*) earlier, Equation 13, and that consistency with this result would require that

*α*is determined by Equation 13.

The two approximate solutions, Equation 22 and Equation 23, are both accurate in the intermediate regime *γ* in Equation 22 by matching the two solutions. This matching implies that*λ* and *Ns* of the rescaled model. Since rates have units of inverse time, this expression has to be multiplied by *s* to obtain the rate in units of inverse generations.

However, Equation 24 does not describe the rate accurately, as is obvious from the comparison with simulation results shown in Figure 4A. The plot shows the rescaled ratchet rate *αNse*^{−λ}], indicated by the black line. The plot shows clearly that the simulation results often differ from the prediction of Equation 24 by a large factor. It seems as if α needs to depend on *λ*, as we already noted above when comparing the variance of *p*(*z*) to Equation 13. In fact, fixing α via Equation 13 improves the agreement substantially, but still does not describe the simulations quantitatively.

The reason for the discrepancy is the time delay between *z*, which we quantified by calculating the correlation between *z*(*τ* + Δ*τ*). Hence we cannot use an approximation where *z*, but must calculate *z*. If the fittest class is the only one that is strongly stochastic, we can calulate *z*(*ρ*), *ρ* ≤ *τ,* by integrating the deterministic evolution equations for *x _{k}* with

*k*≥ 1 with

*z*(

*ρ*) as an external forcing.

Equation 17 now depends not only on *z*(*τ*), but on all *z*(*ρ*) with *ρ* ≤ *τ* and cannot be mapped to a diffusion equation. Nevertheless, it corresponds to a well-defined stochastic integral, known as a path integral in physics (Feynman and Hibbs 1965), which is amenable to systematic numerical approximation. To introduce path integrals, it is useful to discretize Equation 17 in time and express *z*(*ρ _{i}*) in terms of the state at time

*ρ*

_{i}_{−1}=

*ρ*− Δ

_{i}*τ*and the earlier time points. For simplicity, we use the notation

*z*for

_{i}*z*(

*ρ*),

_{i}*ρ*with

_{j}*j*<

*i*. In the limit Δ

*τ*→ 0, this difference equation converges against Equation 17 interpreted in the Itô sense since the

*z*-dependent prefactor of the noise term is evaluated at

*ρ*

_{i}_{−1}rather than at an intermediate time point between

*ρ*

_{i}_{−1}and

*ρ*. We can express this transition probability

_{i}*P*(

_{τ}*z*|

*z*

_{0}) between the initial state

*z*

_{0}and the final state

*z*=

*z*as a series of integrals over all intermediate states

_{m}*z*for 0<i<

_{i}*m*,

*η*drawn from a standard Gaussian (Lau and Lubensky 2007). Hence

_{i}*τ*, the transition probability can therefore be written as

*z*is the limit of

_{ρ}*z*

_{0}and

*z*=

_{τ}*z*. The functional

*S*({

_{λ}*z*}) in the exponent closely corresponds to the “action” in physics (Feynman and Hibbs 1965), which is minimized by classical dynamics. Here minimization of the “action” defines the most likely trajectory. Note that

_{ρ}*S*({

_{λ}*z*}) depends on the entire path {

_{ρ}*z*} with 0 ≤

_{ρ}*ρ*≤

*τ*, while the functional itself depends only on

*λ*. The strength of genetic drift appears as a prefactor of

*S*({

_{λ}*z*}) in the exponent.

_{ρ}The most likely path *z*_{0} and *z _{τ}* in time

*τ*can be determined either by solving the Euler–Lagrange equations or by numerical minimization, see below. Along with the functional,

*λ*. Given this extremal path, we can parameterize every other path connecting

*z*

_{0}and

*z*as

_{τ}*δz*vanishes at both endpoints (

_{ρ}*δz*

_{0}=

*δz*= 0). Denoting the minimal action associated with

_{τ}^{−1}factor is equal to the integral over the fluctuations, which in general depends on

*Ns*in

*δS*({

*δz*},

_{ρ}*z*,

_{τ}*z*

_{0}) is independent of the final point

*z*,

_{τ}*P*(

_{τ}*z*|

_{τ}*z*

_{0}) with respect to

*z*. In the general case, calculating the fluctuation integral is difficult, and we determine it here by analogy to the history-independent solution presented above Equations 18–24.

_{τ}If the stochastic dynamics admits an (approximately) stationary distribution, *P _{τ}*(

*z*|

_{τ}*z*

_{0}) becomes independent of

*τ*and

*z*

_{0}and coincides with the steady-state probability distribution

*p*(

*z*). It therefore becomes the analog of Equation 23, which for arbitrary diffusion equations is given by the inverse diffusion constant (the prefactor

*z*

^{−1}), multiplied by an exponential quantifying the trade-off between deterministic and stochastic forces. In this path-integral representation, the exponential part is played by

*z*and

*λ*only. The prefactor is independent of the selection term and can hence be determined through the analogy to the Markovian case discussed above

*σ*

^{2}is given by Equation 13. Note that this solution is not valid very close to the absorbing boundary since this boundary is not accounted for by the path integral, at least not without some special care. As in the history-independent case discussed above, this approximate distribution should be thought of as the time-independent “bulk” distribution in

*P*(

*z*,

*τ*) =

*e*

^{−γτ}

*p*(

*z*). To determine the rate extinction rate

*γ*, we again need to understand how probable it is that a trajectory actually hits

*z*= 0, given that it has come pretty close.

To this end, we need a local solution of Equation 17 in the boundary layer *z* = 0, its fate, *i.e.*, whether it goes extinct or returns to *α*, and hence *γ*, by matching of the boundary solution to Equation 30 in the regime *α* and *γ* by the matching requirement *σ*^{2} is given by Equation 13 and depends on *λ* and *Ns* as *γ* is in units of *s* and needs to be multiplied by *s* for conversion to units of inverse generations. In contrast to Markovian case above, the variance of the “bulk” is no longer simply related to strength of selection near the *z* = 0 “boundary.”

Since we don’t know how to calculate *Model and Methods*. Examples of numerically determined most likely path and the corresponding trajectory of the mean fitness are shown in Figure 5 for different values of *λ*. Generically, we find a rapid reduction of the size *z* of the size of the fittest class such that the mean fitness has only partially responded. The inset of Figure 5A shows how changes in mean fitness *z* for different *λ*. For large *λ*, the mean fitness changes only very slowly with *z*, which increases the probability of large excursions and hence the rate of the ratchet.

This numerically determined minimal action *u*, *s*, and *N* with *λ* = *u*/*s* ranging from 1 to 30. Note that the vertical shift of the black line relative to the data points depends on the prefactor, which we have approximated. Hence we should not expect agreement better than to a factor of ∼2. The important point is that the exponential dependence of the rate on

Previous studies of Muller’s ratchet suggested that the rate depends exponentially on *αNse*^{−}* ^{λ}* (Jain 2008). To relate this to our results, we determined “Haigh’s factor”

*α*numerically from

*α*(

*λ*) drops from around 0.8 to 0.3 as

*λ*increases from 1 to 30. The previously used values 0.5–0.6 for

*α*correspond to

*λ ≈*6. Using

*s*. In units of generations, the mean time between clicks is given by

*ζ*(

*λ*) is determined by Equation 13 and the factor 2.5 is introduced to approximate the part of the prefactor that is independent of

*N*,

*s*, or λ. A direct comparison of this expression with simulation results is shown in File S1.

## Conclusion

The main difficulty impeding better understanding of even simple models of evolution is the fact that rare events involving a few or even single individuals determine the fate of the entire population. The important individuals are those in the high fitness tail of the distribution. Fluctuations in the high fitness tail propagate toward more mediocre individuals, which dominate a typical population sample.

We have analyzed the magnitude, decay, and propagation of fluctuations of the fitness distribution in a simple model of the balance between deleterious mutations and selection. In this model, individuals in the fittest class evolve approximately neutrally. Fluctuations in the size of this class propagate to the mean, which in turn generates a delayed restoring force opposing the fluctuation. We have shown that the variance of the fluctuations in the population *n*_{0} of the top bin is proportional to *n*_{0}/*s* and increases as log *λ* with the ratio *λ* of the mutation rate *u* and the mutational effect *s*. Fluctuations of *n*_{0} perturb the mean after a time ∼*s*^{−1}log *λ*. These two observations have a straightforward connection: Sampling fluctuations can accumulate without a restoring force for a time *s*^{−1}log *λ*. During this time, the typical perturbation of the top bin by drift is *n*_{0}*s*^{−1}log *λ*. We have used these insights in the coupling between *n*_{0}, the mean fitness, and the resulting delayed restoring force on fluctuations of *n*_{0} to approximate the rate of Muller’s ratchet.

The history dependence of the restoring force has not been accounted for in previous analysis of the rate of Muller’s ratchet (Haigh 1978; Stephan *et al.* 1993; Gordo and Charlesworth 2000; Jain 2008) who introduced a constant factor to parameterize the effective strength of the selection opposing fluctuations in the top bin, or Waxman and Loewe (2010), who replaced all mutant classes by one effective class and thereby mapped the problem to the fixation of a deleterious allele. We have shown that to achieve agreement between theory and numerical simulation one must account for the delayed nature of selection acting on fluctuations. Comparing our final expression for the ratchet rate with that given previously (Stephan *et al.* 1993; Gordo and Charlesworth 2000; Jain, 2008), the history dependence manifests itself as a decreasing effective strength of selection with increasing *λ* = *u*/s. This decrease is due to a larger temporal delay of the response of the mean fitness to fluctuations of the least loaded class. History dependence is a general consequence of projecting a multidimensional stochastic dynamics onto a lower dimensional space (here, the size *n*_{0} of the fittest class). Such memory effects can be accounted for by the path-integral formulation of stochastic processes, which we used to approximate the rate of Muller’s ratchet.

Even though the model is extremely simplistic and the sensitive dependence of the ratchet rate on poorly known parameters such as the effect size of mutations, population size, and mutation rate precludes quantitative comparison with the real world, we believe that some general lessons can be learned from our analysis. The propagation of fluctuations from the fittest to less fit individuals is expected to be a generic feature of many models and natural populations. In particular, very similar phenomena arise in the dynamics of adapting populations driven by the accumulation of beneficial mutations (Tsimring *et al.* 1996; Rouzine *et al.* 2003, 2008; Cohen *et al.* 2005; Desai and Fisher 2007; Neher *et al.* 2010; Hallatschek 2011). The speed of these traveling waves is typically determined by stochastic effects at the high fitness edges. We expect that the fluctuations of the speed of adaptation can be understood and quantified with the concepts and tools that we introduced above.

Populations spread out in fitness have rather different coalescence properties than neutral populations, which are described by Kingman’s coalescent (Kingman 1982). These differences go beyond the familiar reduction in effective population size and distortions of genealogies due to background selection (Charlesworth *et al.* 1993; Higgs and Woodcock 1995; Walczak *et al.* 2011). The most recent common ancestor of such populations most likely derives from this high fitness tail and fluctuations of this tail determine the rate at which lineages merge and thereby the genetic diversity of the population (Brunet *et al.* 2007; Rouzine and Coffin 2007; Neher and Shraiman 2011). Thus, quantitative understanding of fluctuations of fitness distributions is also essential for understanding nonneutral coalescent processes.

Generalizing the analysis of fluctuations of fitness distributions to adapting “traveling waves” and the study of their implications for the coalescent properties of the population are interesting avenues for future research.

## Acknowledgments

We are grateful for stimulating discussions with Michael Desai, Dan Balick, and Sid Goyal. R.A.N. is supported by the European Research Council through Stg-260686, and B.I.S. acknowledges support from National Institutes of Health under grant GM-086793. This research was also supported in part by the National Science Foundation under grant no. NSF PHY11-25915.

## Footnotes

*Communicating editor: W. Stephan*

- Received April 18, 2012.
- Accepted May 11, 2012.

- Copyright © 2012 by the Genetics Society of America

Available freely online through the author-supported open access option.