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

Linear stability analysis

In the limit of large populations, the fluctuations of x_k around the deterministic steady state x−k can be analyzed in linear perturbation theory. In other words, we express deviations from the steady state as δxk=xk−x−k and expand the deterministic part of Equation 3 to order δxk2. This expansionddτδxk=−kδxk+λδxk−1+x−k∑m=0(m−λ)δxm=∑mLkmδxm (4)defines a linear operator L_km. A quick calculation shows that L_km has eigenvalues κ_i = −i with i = 0, 1, 2.... The right eigenvector corresponding to κ₀ = 0 is simply ψk(0)=x−k, while the right eigenvectors for i > 0 are given byψk(i)=x−k−i−x−k, (5)where k numbers the coordinate of the vector. This is readily verified by direct substitution (note that x−i=0 for i < 0).

The eigenvector ψk(0) corresponds to population size fluctuations that in our implementation are controlled by a carrying capacity. The eigenvalue associated with this mode in the computer simulation is large and negative and need not be considered here (see Model and Methods). All other eigenvalues are negative, which is to say that x−k is a stable solution.

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 δxk(τ)=∑jψk(j)aj(τ) and the associated amplitudes a_j(τ) will decay exponentially in time with rate j (remember that the unit of time is s⁻¹). Since the amplitudes are projections of δx_k onto the left eigenvectors of L_mk, we need to know those as well. For κ₀ = 0, the left eigenvector is simply φk(0)=1, while the other left eigenvectors are given byφk(i)=(−1)k−ieλλi−k(i−k)!, 0≤k≤i, (6)and φk(i)=0 for 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:ddτδxk=∑mLkmδxm+x−kNsηk. (7)Note that we approximated the strength of noise by its value at equilibrium. This approximation is justified as long as we consider only small deviations from the equilibrium. The full x_k-dependent noise term is reintroduced later when we turn to Muller’s ratchet. Substituting the representation of δxk(τ)=∑iψk(i)ai(τ) and projecting onto the left eigenvector φk(j), we obtain the stochastic equations for the amplitudesddτaj(τ)=−jaj(τ)+∑kφk(j)x−kNsηk(τ) . (8)Each noise term η_k contributes to every a_j and induces correlations between the a_j, but each amplitude can be integrated explicitly:aj(τ)=∫−∞τdτ′e−j(τ−τ′)∑kφk(j)x−kNsηk(τ′). (9)The covariances of different amplitudes are evaluated in File S1 and found to be〈ai(τ)aj(τ+Δτ)〉=e−jΔτi+j∑kφk(i)φk(j)x−kNs. (10)However, we are not primarily interested in the covariance properties of the amplitudes of eigenvectors, but expect that the fluctuations of the fittest class and fluctuations of the mean fitness are important for the rate of Muller’s ratchet and other properties of the dynamics of the population. To this end we express δx₀(τ) and k−(τ) asδx0(τ)=∑j>0ψ0(j)aj(τ)=−e−λ∑j>0aj(τ) (11)δk−(τ)=−∑j>0,kkψk(j)aj(τ)=−∑j>0jaj(τ) . (12)With Equation 10, we can now calculate the desired quantities. The calculations required to break down the multiple sums to interpretable expressions are lengthy but straightforward and detailed in File S1. Below, we present and discuss the results obtained in File S1.

Fluctuations of x₀ and the mean fitness

For the variance of the fittest class 〈δx02〉 and more generally its autocorrelation, we find〈δx0(0)δx0(τ)〉=e−λNs∫01dθθGλ(θ,τ)Gλ(θ,τ)=eλθ2e−τ−λ(1+e−τ)θ−e−λθ−e−λθe−τ+ 1. (13)The variance of the fittest class (τ = 0 in the above expression) is therefore σ2=(x−0/Ns)ζ2(λ), where ζ2(λ)=∫01dθθ−1Gλ(θ,0) is the standardized variance of the top bin, which depends only on λ. For small λ, it simplifies to ζ2(λ)≈12λ+O(λ2). This limit corresponds to x−0 close to 1 with only a small fraction of the population carrying deleterious mutations x₁ ≈ λ = 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 σ² is ∼x−0logλ/Ns. The full autocorrelation function is shown in Figure 2A for different values of λ 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 λ.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2 

(A) The covariance of the size of the fittest class x₀(0) with x₀(τ) a time τ later. The normalized autocorrelation of x₀ increases with λ. (B) The covariance of x₀(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 Nx₀s = 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〈δk−(0)δk−(τ)〉=λeλNs∫01dθIλ(θ,τ)Iλ(θ,τ)=e−τeλθ2e−τ−λ(1+e−τ)θ(θ+λθ(θe−τ−1)(θ−1)), (14)which asymptotes to (4Nsx−0)−1 for large λ at τ = 0. It is hence inversely proportional to the size of the fittest class x₀. For large λ, x₀ 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 k− are driven by fluctuations of the fittest class x₀, we expect a strong correlation between those fluctuations (Etheridge et al. 2007). Furthermore, fluctuations of x₀ should precede fluctuations of the mean. These expectations are confirmed by the analytic result〈δx0(0)δk−(τ)〉=λNs∫01dθHλ(θ,τ), (15)whereHλ(θ,τ)={(θ−1)e−τ+e−τλθ2−λ(1+e−τ)θ+e−τ−e−τλθ,τ>0(eτθ−1)eeτλθ2−λ(1+eτ)θ+e−λθ,τ<0. (16)This expression is shown in Figure 2B for different values of λ. The cross correlation 〈δx0(0)δk−(τ)〉 is asymmetric in time: With larger λ, 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 x−0 of the fittest class, the model displays very different behavior. If Nsx−0 is comparable to or smaller than 1, the ratchet clicks often without settling to a quasi-equilibrium in between clicks. This limit has been studied in Rouzine et al. (2008). Conversely, if Nsx−0≫1, ratchet clicks are rare and the system stays a long time close to its quasi-equilibrium state x−k. Such a scenario, taken from simulations, is illustrated in Figure 3. Figure 3A shows the distribution of x₀ prior to the click, while Figure 3B shows the realized trajectory that ends at x₀ = 0. Prior to extinction, x₀(τ) 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 x−0, the final large excursion away from the steady is short and happens in a few units of rescaled time. Translated back to generations, the final excursion took a few hundred generations (s = 0.01 in this example).

• Download figure
• Open in new tab
• Download powerpoint

Figure 3 

An example of a click of the ratchet with N = 5 × 10⁷, s = 0.01, and λ = 10, corresponding to an average size of the fittest class . (A) The distribution of x₀ averaged over the time prior to extinction. (B and C) The trajectory of x₀(τ), with the part of the trajectory that ultimately leads to extinction magnified in C. The final run toward x₀ = 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 isddτx0(τ)=δk−(τ)x0(τ)+x0(τ)Nsη0(τ), (17)where δk−=λ−k−. The restoring force δk− depends on x₀, as well as on the size of the other classes x_k. For sufficiently large λ, x₀ is much smaller than x_k with k ≥ 1, such that the stochastic force is most important for x₀. The dynamics of x_k, k ≥ 1, is approximately slaved to the stochastic trajectory of x₀(τ). We can therefore try to find an approximation of δk−(τ) in terms of x₀(τ) 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₀ is delayed with the delay increasing ∝ log λ. The latter observation implies that the restoring force on x₀ will depend mainly on the values of x₀ 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 δk−(τ) depends on the instantaneous value of x₀(τ) via δk−(τ)=α(1−x0(τ)/x−0) (Stephan 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 x−0 and the stochastic dynamics is equivalently described by a diffusion equation for the probability distribution P(x₀, τ),∂∂τP(z,τ)=12Ns∂2∂z2zP(z,τ)−α∂∂z(1−z/z−)zP(z,τ), (18)where we have denoted x₀ 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, Nsz−≫1, clicks of the ratchet occur on much longer time scales than the local equilibration of 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−γP(X>z)=12Ns∂∂zzp(z)−αz(1−z/z−)p(z), (19)where P(X>z)=∫z∞dz′p(z′), which is ≈1 for z<z− and rapidly falls to 0 for z>z−. To obtain the rate γ, we solve this equation in a regime of small z≪z−, where the term on the left is important but constant, and in a regime z ≫ (Ns)⁻¹, 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/2Ns), force field A(z) (here equal to −αz(1−z/z−)), and a constant C−C=∂∂zD(z)p(z)+A(z)p(z) (20)with solutionp(z)=1D(z)e−∫0zdyA(y)D(y)[β−C∫0zdye∫0ydy′A(y′)D(y′)]. (21)Note that this solution is inversely proportional to the diffusion constant, while the dependence on selection is accounted for by the exponential factors. For z≪z−, C = γ and 2αNs∫0zdy (1−z/z−)≈2αNsz, such thatp(z)≈γe2Nsαz−1αz, z≪z−, (22)where β is fixed by the boundary condition that p(z) is finite at z = 0. Note that γ = p(0)/2Ns relates the rate of extinction to p(0). To determine the latter we need to match the z≪z− regime to the bulk of the distribution z≈z−. As can be seen from Equation 22, the constant term C is unimportant in this regime (e^2Nsαz ≫ 1). Setting C = 0 in Equation 21, we findp(z)≈z−zNsαz−πexp[−αNs(z−z−)2z−], zNsα≫1. (23)The integration constant β in Equation 21 is fixed by the normalization. Since p(z) is concentrated around z=z− and has a Gaussian shape around z−, the normalization factor is simply 1/2πσ2, where σ2=z−/2αNs is the variance of the Gaussian. The factor z−/z corresponds to the 1/D(z) term, scaled so that it equals 1 in the vicinity of z−. Note that we have already calculated the variance 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 (Nsα)−1≪z≪z−, which allows us to determine the rate γ in Equation 22 by matching the two solutions. This matching implies thatγ=z−Nsα3πe−αNse−λ, (24)which agrees with result obtained previously (Jain 2008). Note that this rate depends only on the parameters λ 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 γ×π/z−Nsα3, which according to Equation 24 should be simply exp[−α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.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4 

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 δk− and z, which we quantified by calculating the correlation between δk−(τ) and z(τ + Δτ). Hence we cannot use an approximation where δk− depends on the instantaneous value of z, but must calculate δk− from the past trajectory of z. If the fittest class is the only one that is strongly stochastic, we can calulate δk−(τ) for a given trajectory 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_i for z(ρ_i),zi−zi−1=Δτδk−i−1zi−1+zi−1ΔτNs ηi−1, (25)where δk−i−1 depends on all previous time points ρ_j with 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 ρ_i. We can express this transition probability P_τ(z|z₀) between the initial state z₀ and the final state z = z_m as a series of integrals over all intermediate states z_i for 0<i<m,Pτ(z|z0)=∫∏i=1m−1dziPΔτ(z|{zj}j<m)PΔτ(zm−1|{zj}j<m−1)⋯PΔτ(z1|z0). (26)Each of these infinitesimal transitions corresponds to solutions of Equation 25 with η_i drawn from a standard Gaussian (Lau and Lubensky 2007). HencePΔτ(zi|{zj}j<i)=Ns2πΔtzi−1exp[−Ns(zi−zi−1−Δτzi−1δk−i−1)22Δτzi−1]. (27)In the limit of many intermediate steps and small Δτ, the transition probability can therefore be written asPτ(z|z0)=∫Dzρ exp[−Ns∫0τdρ[z˙ρ−zρδk−ρ]22zρ]=∫Dzρ e−NsSλ({zρ}), (28)where Dz_ρ is the limit of ∏i=1mdzi(2πΔτzi−1/Ns)−1/2 known as the path-integral measure, and we have replaced the discrete time index by its continuous analog. The path integral extends over all continuous path connecting the endpoints z₀ 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ρ* connecting the endpoints points z₀ and z_τ in time τ can be determined either by solving the Euler–Lagrange equations or by numerical minimization, see below. Along with the functional, zρ* depends only on λ. Given this extremal path, we can parameterize every other path connecting z₀ and z_τ as zρ=zρ* + δzρ, where δz_ρ vanishes at both endpoints (δz₀ = δz_τ = 0). Denoting the minimal action associated with zρ* by Sλ*(zτ,z0), we havePτ(zτ|z0)=e−NsSλ*(zτ,z0)∫Dδzρ e−Ns δSλ({δzρ},zτ,z0)=N−1e−NsSλ*(zτ,z0), (29)where N⁻¹ factor is equal to the integral over the fluctuations, which in general depends on zρ*. The prefactor Ns in e−Ns δS({δzρ},zτ,z0) implies that deviations from the optimal path are suppressed in large populations. If δS({δz_ρ}, z_τ, z₀) is independent of the final point z_τ, N can be determined by the normalizing P_τ(z_τ|z₀) 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₀) becomes independent of τ and z₀ 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⁻¹), multiplied by an exponential quantifying the trade-off between deterministic and stochastic forces. In this path-integral representation, the exponential part is played by e−NsSλ*(z), where Sλ*(z). is a function of the final point z and λ only. The prefactor is independent of the selection term and can hence be determined through the analogy to the Markovian case discussed abovep(z)≈z−z12πσ2e−NsSλ*(z), Nsz≫1, (30)where the normalization is obtained by assuming an approximately Gaussian distribution around the steady-state value z− and the variance σ² 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≪z− as already obtained for the history independent case in Equation 22. Once a trajectory comes close to z = 0, its fate, i.e., whether it goes extinct or returns to z≈z−, is decided quickly. Hence we can make an instantaneous approximation for δk−, which does depend on the past trajectory, but for the time window under consideration it is simply a constant, α, yet to be determined. Having reduced the problem to Equation 22 we can determine α, and hence γ, by matching of the boundary solution to Equation 30 in the regime (Ns)−1≪z≪z−, where both are accurate. The matching condition isz−z12πσ2e−Ns(Sλ*(0)+z∂zSλ*(z)|z=0)=γαze2αNsz, (31)which determines α and γ by the matching requirement 2α=−∂zSλ*(z)|z=0 andγ=|∂zSλ*(z)|z−e−NsSλ*(0)22πσ2≈Sλ*(0)8πσ2e−NsSλ*(0), (32)where we approximated |∂zSλ*(z)|z=0≈Sλ*(0)/z−. The variance σ² is given by Equation 13 and depends on λ and Ns as σ2=z−ζ2(λ)/Ns. Note that γ 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 Sλ*(0) or the most likely path zρ* analytically, we determined discrete approximations to zρ* numerically as described in 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 δk− are related to 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.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5 

(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 Sλ*(0) together with the approximation Equation 32 describes the rate of the ratchet, as determined in simulations, extremely well. Figure 4B shows the same simulation data as Figure 4A, but this time rescaled by 8πσ2(λ)/Sλ*(0) as a function of NsSλ*(0). After this rescaling, we expect all data points to lie on the same curve given by exp[−NsSλ*(0)], as is indeed found for many different values of 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 NsSλ*(0) is correctly captured by Equation 32.

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 Sλ*(0) and plotted it in Figure 5B. We find that α(λ) 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λ*(0)eλ as shown in Figure 5B, we can recast Equation 32 into its traditional form and undo the scaling with s. In units of generations, the mean time between clicks is given byTclick≈2.5ζ(λ)α(λ)sNse−λeNsα(λ)e−λ, (33)where ζ(λ) 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.