The extinction time under mutational meltdown driven by high mutation rates

Abstract Mutational meltdown describes an eco‐evolutionary process in which the accumulation of deleterious mutations causes a fitness decline that eventually leads to the extinction of a population. Possible applications of this concept include medical treatment of RNA virus infections based on mutagenic drugs that increase the mutation rate of the pathogen. To determine the usefulness and expected success of such an antiviral treatment, estimates of the expected time to mutational meltdown are necessary. Here, we compute the extinction time of a population under high mutation rates, using both analytical approaches and stochastic simulations. Extinction is the result of three consecutive processes: (a) initial accumulation of deleterious mutations due to the increased mutation pressure; (b) consecutive loss of the fittest haplotype due to Muller's ratchet; (c) rapid population decline toward extinction. We find accurate analytical results for the mean extinction time, which show that the deleterious mutation rate has the strongest effect on the extinction time. We confirm that intermediate‐sized deleterious selection coefficients minimize the extinction time. Finally, our simulations show that the variation in extinction time, given a set of parameters, is surprisingly small.

where k := min b := min x2N 8 < : x such that for ✓ := µ(1 s)/s and := ✓ k. Here N denotes the total population size, which Gessler assumed to be constant. In our model this translates to N = K, as our population remains at its carrying capacity during the ratchet phase.
For simplicity let us define and which appear in the definitions of the value k in (14) and b in (15). Since k, b 2 N, Gessler's speed of the ratchet, eq. (13), can be non-monotonous with respect to the parameters, as shown in Figure 2. For instance, both coefficients k and b are monotonous and non-decreasing functions of the mutation rate µ. However, the following situation might occur: i.e., with increasing µ, b may take on a larger value while k remains constant, which results in a sudden decrease of the ratchet speed (which should monotonously increase with the mutation rate), To avoid this behavior, we extend the domain of the functions f k , f b as defined in (16)-(17), to positive real numbers. We define where we use the Euler Gamma function which generalises the factorial operation over the real numbers (Abramowitz and Stegun 1972). In particular, f k (x) = f k (x) for all x 2 N. The extension of the function f b (x) in (17) to positive real numbers needs more care because x appears in the upper bound of summation of Defining n x,k := bx k 1c 2 N, i.e. n x,k  x k 1 < n x,k + 1, the function represents an interpolation of g(x k 1) := P x k 1 i=0 i i! on [0, 1). Indeed, when x k 1 2 N, equation (21) corresponds to the original sum, and g(n) = g(n) for all n 2 N. Therefore, the following function extends f b to positive real numbers: Since both f k and f b are increasing functions of x 2 N, the problem of finding the minimum integers k, b such that f k (k), f b (b) 0, as required in equations (14) and (15) Similar to the original formula by Gessler (1995), the above expression is well defined only for suitable choices of the parameters µ, s, N . Figure 2 in the main text compares the ratchet speed in its original discrete version by Gessler (1995) with our smooth monotonic extension.

Appendix S2 Description of the computational analysis
The simulation of the population dynamics and mutation accumulation was performed using the script 'simulate_mutational_meltdown.jl'. First, the parameters are initialised and a founder  (3) in the Population-size control step, if the population has grown over its carrying capacity, randomly selected individuals are removed until population size is reduced to K individuals. Possible outputs are the extinction time, the number or the distribution of mutations every generation, the beginning and end time points of the pre-ratchet, ratchet and meltdown phases, and the times until the fittest classes are lost. We ran the simulations for the parameter regime described in the manuscript using the script 'launcher_script.jl'. Here, we were interested in the times of pre-ratchet, ratchet and meltdown phase and the times until loss of the fittest class. The output is a data frame that contains the simulation output and data frame that contains the parameters used for this specific simulation. Each specific simulation can be identified and reproduced via a seed from the random number generator.
The analysis was performed in the Jupyter notebook 'Computational_analysis.ipynb', and the simulations were implemented in the programming language 'julia', version 1.6 (Bezanson et al. 2017;Besançon et al. 2021).

Power-law exponent B A
Selection coefficient K = 10 4 K = 10 3 K = 10 2 Figure S2: Numerical estimation of the power law relationship between the mean extinction time and the mutation rate T E ⇡ ↵µ . We used ordinary least squares regression to estimate (A) the scaling factor ↵ and (B) the power-law exponent and the respective standard errors depending on the selection coefficient s and for different carrying capacities K.

Appendix S4 Simulation data for carrying capacities K=100, 10000
In addition to the simulation data presented in Section 3.6 with a carrying capacity of K = 1000, we also performed simulations for carrying capacities of K = 100 and K = 10000, see Figures S3 and S4. They showed similar results and are consistent with our expectation on the dependence of the extinction time on the carrying capacity as described in Section 3.4, also shown in Figure  3B. The contour lines give parameter combinations with equal mean extinction time. Our analytical estimation is not applicable and returns an infinite extinction time if the mutation-selection balance is stable, which is indicated by the blue line in panel A. This is not observed in stochastic simulations shown in panel B. The predicted optimal selection coefficient for which the extinction time is minimal (green line panel A) is smaller than the one observed in stochastic simulations (green dots panel B). C The relative error of our analytical estimation is small in general but becomes large as the boundary of our parameter regime is approached. D The coefficient of variation of the extinction time in the simulations is comparably small, but positively correlated with the selection coefficient. (Other parameters: founder population size N 0 = 20, wild-type reproduction rate w 0 = 2 and n = 20 simulation runs.)