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.
Deleterious mutation–selection balance. The population is distributed among classes of individuals carrying k deleterious mutations. Classes with few mutations grow due to selection (red arrows), but lose individuals through mutations (green arrows), while classes with many mutations are selected against but replenished by mutations.
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 nk of the population carries k deleterious mutations, is given by
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−1nk. 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,
To determine the ratchet rate, the population was initialized with its steady-state expectation
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
Results and Discussion
Fluctuations of the size n0 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
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 xk = N−1 nk rather than numbers of individuals, and obtain
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 xk around the deterministic steady state
The eigenvector
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
Fluctuations of x0 and the mean fitness
For the variance of the fittest class
(A) The covariance of the size of the fittest class x0(0) with x0(τ) a time τ later. The normalized autocorrelation of x0 increases with λ. (B) The covariance of x0(0) with the mean fitness at time τ in the past or future. One observes a pronounced asymmetry, showing that fluctuations of the fittest class propagate toward the bulk of the fitness distribution and results in delayed fluctuations. Simulation results are shown as dashed lines; theory curves are solid. In all cases, s = 0.01 and Nx0s = 100. Note that time is measured in units of 1/s, which is the natural time scale of the dynamics.
In a similar manner, we can calculate the autocorrelation of the mean fitness
If fluctuations of the mean fitness
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
An example of a click of the ratchet with N = 5 × 107, s = 0.01, and λ = 10, corresponding to an average size of the fittest class . (A) The distribution of x0 averaged over the time prior to extinction. (B and C) The trajectory of x0(τ), with the part of the trajectory that ultimately leads to extinction magnified in C. The final run toward x0 = 0 takes a few time units, as expected from the results for the correlation functions, which suggest a (rescaled) correlation time of ∼log λ. Note this time corresponds to a few hundred generations since s = 0.01.
In rescaled time, the equation governing the frequency of the top class is
The two approximate solutions, Equation 22 and Equation 23, are both accurate in the intermediate regime
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
The ratchet rate from simulation vs. prediction. Both A and B show the ratchet rate γ, rescaled with a prefactor to isolate the exponential dependence predicted by analytic approximations; λ is color-coded. (A) Comparison of simulation results with the prediction of Equation 24, which is shown as a straight line. The approximation works only for a particular value of λ, for otherwise the exponential dependence on is not predicted correctly. (B) Comparison of simulation results with the prediction of Equation 32, again indicated by the straight line. The exponential dependence of rate on
is well confirmed by simulation results.
The reason for the discrepancy is the time delay between
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 zi for z(ρi),
The most likely path
If the stochastic dynamics admits an (approximately) stationary distribution, Pτ(zτ|z0) becomes independent of τ and z0 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
To this end, we need a local solution of Equation 17 in the boundary layer
Since we don’t know how to calculate
(A) The most likely path to extinction of the fittest class and the concomitant reduction of the mean fitness for different λ are plotted against time. Times are shifted such that
. The inset shows the mean fitness
plotted against
for different values of λ. (B) Haigh’s factor
as a function of λ determined numerically.
This numerically determined minimal action
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
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 n0 of the top bin is proportional to n0/s and increases as log λ with the ratio λ of the mutation rate u and the mutational effect s. Fluctuations of n0 perturb the mean after a time ∼s−1log λ. These two observations have a straightforward connection: Sampling fluctuations can accumulate without a restoring force for a time s−1log λ. During this time, the typical perturbation of the top bin by drift is
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 n0 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.