Model I: Target site resistance
Cost-free resistance
In the absence of resistance, release of a driving Y in this model leads either to population suppression or population elimination, depending on the transmission rate of the driving Y (m) and the intrinsic rate of increase of the population (R
m
). Assuming that m = 0.95 (consistent with the results of Galizi et al. [11, 12]), and R
m
= 6 (consistent with the analyses of Deredec et al. [10]), then the model population tends to elimination, with time course calculated from (2) and shown in Fig. 1 (assuming an initial release of h
0 = 0.05 N
0, i.e. 5% of the pre-release equilibrium population size). As can be seen, the population size drops rapidly, such that it is only \(\cong 30\%\) of its original size after 10 generations, and \(\cong 0.2\%\) after 30 generations. The number of driving Y males initially increases due to the drive, and then decreases due to the reduction in total population size. The total number of driving Y males born over the time to elimination is \(\cong 4\;N_{0}\), and the number of daughters of these males is \(\cong 0.2\;N_{0}\) (see Additional file 1: Eq. A2.12). As it is assumed that mutations to target site resistance can only occur in the meiotic cells of driving Y males, and would only be transmitted to their daughters, these numbers give some indication of the opportunity for resistance to arise.
The stochastic model for the evolution of resistance requires one more parameter, the product of the initial population size (N
0) and mutation parameter u (i.e., the fraction of female progeny of a driving Y male H(t) that inherit a mutant resistant X chromosome (0 ≤ u ≤ 1). If one considers the baseline value of uN
0 = 1, then the probabilities of three possible outcomes can be calculated from the model: no resistant mutation arises and the population is eliminated (\(1 - P_{Mut} \cong 80\%\) for baseline parameter values); at least one mutation does arise, but none establish, and the population is eliminated \((P_{Mut} - P_{1} \cong 13\% )\); and at least one mutation arises and establishes, and the population persists indefinitely \((P_{1} \cong 7\% ).\) The overall probability the population is eliminated is 93%. For these parameter values, the population is usually eliminated by the driving Y before a single resistant mutation occurs, and even when one does arise, most of the time the mutation is lost and the population is still eliminated. In the Additional file 1: Figure A3.1 presents some exemplar runs from the Gillespie simulations (106 runs, N
0 = 106) for which at least one resistant mutation survives stochastic loss and the population recovers.
In those cases when resistance does evolve, it evolves quickly, and the number of females in the target population is suppressed for a relatively short time (e.g., it remains below 33 and 5% of its pre-release equilibrium value for 17.6 (median 17.3, interquartile range 16–19) and 12.3 (median 11.9, interquartile range 10.5–13.6) generations on average, respectively (full simulation model used, 105 runs). To put these numbers into context, An. gambiae mosquitoes may have 10–18 generations per year, depending on temperature [31, 32]. As an aside, it is noted for comparison that if, for example, m = 0.9, which is insufficient to eliminate the target population, resistance will always evolve (assuming u > 0), and the population will be suppressed below 33% of its initial density for an average of 90.4 (median 62.5, interquartile range 28–122) generations before it recovers. It never goes below 5% of its initial value (detailed results in Additional file 1: Figure A3.2, Section A3).
In this model, there are three parameters that affect the probability that resistance evolves before the population is eliminated: the product of the mutation rate and population size (uN
0), the intrinsic rate of increase of the population (R
m
), and the transmission rate of the driving Y (m). Now, the effect of varying each of these parameters individually on the probability of resistance evolving is investigated.
Varying uN
0
From (5) in the limit of low uN
0 ≪ 1, the probability that at least one mutation arises before elimination, P
Mut
, is proportional to uN
0:
$$P_{Mut} \cong 2\left( {1 - m} \right)\left( {uN_{0} } \right)\mathop \int \limits_{0}^{\infty } \frac{{\left( {R_{m} - \left( {R_{m} - 1} \right)N\left( \tau \right)} \right)F\left( \tau \right)H\left( \tau \right)}}{H\left( \tau \right) + M(\tau )}d\tau$$
Since the probability that a single mutation establishes does not depend on uN
0 in (3), the probability a mutation arises and establishes, allowing the population to persist (P
1), is also proportional to uN
0 for uN
0 ≪ 1:
$$P_{1} \cong 2\left( {1 - m} \right)\left( {uN_{0} } \right)\mathop \int \limits_{0}^{\infty } \frac{{p_{{est,F_{R} }} \left( \tau \right)\left( {R_{m} - \left( {R_{m} - 1} \right)N\left( \tau \right)} \right)F\left( \tau \right)H\left( \tau \right)}}{H\left( \tau \right) + M(\tau )}d\tau$$
This agrees with simulation results of Marshall et al. [19] for population-suppressing homing-based gene drive that show a linear relationship between \(1/N_{0}\) and the cleavage-resistant allele generation rate leading to a given elimination probability. For sufficiently small values of uN
0, the conditional probability therefore also does not vary with uN
0, as expected due to the proportionality of both P
1 and P
Mut
on uN
0, such that uN
0 cancels out in (7). Figure 2a confirms these dependencies for varying \(uN_{0} .\)
Varying R
m
Increasing R
m
while keeping everything else the same leads to an increase in the population size, as \(N_{0} = (R_{m} - 1)/\gamma\), and therefore in uN
0, and consequently has much the same effects as an increase in uN
0, as analysed above. Here it is asked, for populations of the same size but that differ in R
m
(and therefore also differ compensatingly in γ), how are they expected to differ in the probability of evolving resistance? Increasing R
m
in this way increases the time for a driving Y to eliminate the population, and therefore increases the opportunity for resistant mutations to arise before elimination (Fig. 2b, black lines). In addition, the probability that a mutation arising at a specified time \(t_{a}\) becomes established (in the absence of others), p
est
(t
a
), increases as \(R_{m}\) increases (Additional file 1: Section A2.1, Figure A2.2a). This is because at higher R
m
the recruitment (birth) rate of new resistant mutants is higher, reducing the probability of stochastic loss. Since at higher R
m
the probability of a mutation arising is higher, and its subsequent probability of surviving stochastic loss is also higher, the overall probability \(P_{1}\) that resistance evolves increases with increasing R
m
(red lines in Fig. 2b).
Varying m
Increasing the transmission rate of the driving Y (m) reduces the probability of at least one mutation arising before elimination (Fig. 2c, black lines). Again, there are several reasons for this. Firstly, there is a factor (1 − m) in the mutation rate because the proportion of mutant females born from driving Y males decreases according to the sex bias (1 − m) (Eq. 3). At the limit of m = 1, no mutations can arise at all because driving Y males only create other males. Secondly, with larger m there is less time for mutations to occur, since the driving Y eliminates the population more quickly. Close to m
crit
(for these parameters, m
crit
= 0.9167), the population is eliminated very slowly, providing more opportunity for a mutation to occur. The probability \(p_{est} \left( {t_{a} } \right)\) that a single resistant mutation (that arises in a female at time t
a
) gets established in a population is also affected by m (Additional file 1: Figure A2.2b). So, with stronger Y drive, the probability that a resistant mutation occurs before elimination is lower, and if one does occur, the probability it establishes is lower, and thus the overall probability of resistance evolving is lower (Fig. 2c, red line).
Costly resistance
Thus far, it has been assumed that resistance is cost-free, with no pleiotropic effects on other fitness components. Now the case of resistance having a cost is considered. This is modelled as a decrease in fertility of females with the resistant gene, and decreased participation in mating for males with the gene (Additional file 1: Section A1.1). There are four genotypes carrying one or more resistant alleles, and the system of equations (Additional file 1: Eq. A1.2) has separate parameters for all of them. However, for simplicity of analysis here, these are collapsed to a single parameter, assuming that heterozygous females have fitness w < 1, and homozygous females and hemizygous males have fitness w
2 (relative to fitness one for wild-types and driving Y males without the resistant gene).
To find the equilibria for a nonzero mutation rate, the time derivatives in the deterministic differential equations in the presence of mutation (Additional file 1: Eq. A1.2) are set to zero. Focussing initially on the deterministic model, Fig. 3 is a contour plot showing the fate of the resistant mutation and the resulting effect on the equilibrium number of females in the population as a function of w and R
m
, assuming m = 0.95 and u = 10−6 (see Additional file 1: Section A1.1b, for calculation of equilibria). The resistant allele fixes deterministically in the population for \(w_{1} < w \le 1\), and establishes at intermediate equilibrium with the wild-type for 0 ≤ w < w
1 (w
1, which is independent of R
m
, is given in Additional file 1: Section A1.1b). The total population goes extinct (shaded area) if fitness is below w = w
ex
(100% suppression of population, solid line), i.e. where net total population growth is not positive. Above w = w
ex
, the total population is nonzero and reduced, and the dotted contour lines in Fig. 3 show the percentage suppression of the total female population size, compared to its size when there is no fitness cost (N
0/2). Therefore, a resistant mutation might establish, and even spread to fixation, but still not rescue the population if it is too costly. As expected, populations with higher R
m
can be rescued by mutations with higher cost than populations with lower R
m
. Moreover, as the strength of Y drive \(m\) is increased, the minimum fitness required for fixation of the mutation, \(\,w_{1}\), decreases, as does w
ex
(above R
m
= 1/w
21
), since the higher the Y drive, the more of an advantage the resistant allele has over the sensitive allele in transmitting to the next generation and in restoring the 50:50 sex ratio, and thus a higher fitness cost can be tolerated (Additional file 1: Figure A1.2). It is also found that there is very weak dependence of these results on the mutation parameter u, assuming it is low (u ≪ 1). Although the focus is on high m (strength of drive) for population suppression and elimination, it is possible to compare results for allele frequencies to the model of Huang et al. [22] for population modification using Y-linked drive and X-linked resistance. For u = 0, using the analytical expressions for allele frequencies for general fitnesses, a region of low m where the resistant X chromosome mutation cannot spread is similarly found, \(m < 1 - \frac{{w_{{F_{R} }} w_{{H_{R} }} }}{{2(2 - w_{{F_{R} }} )}}\), which becomes the condition of Huang et al. (Equation (5), [22]) if their expressions for resistant heterozygote fitnesses are substituted.
Now the effect of fitness costs in the stochastic model is considered. Here, fitness costs do not affect the rate at which mutations arise (P
Mut
), which is constant, but do reduce the probability of surviving stochastic loss and ultimate establishment (P
1, Fig. 2d). For high costs, the branching process cannot be used (as explained in Additional file 1: Section A2.3), and in this region full simulations are used. The conditional probability that at least one mutation survives stochastic loss if one or more mutations arise is also less for higher cost. Thus, higher fitness cost of the mutation (i.e., lower density-dependent recruitment rate for the heterozygotes) results in lower probability that a mutant mosquito will survive early stochastic loss, because it will be less able to pass on the allele before dying. In summary, the costlier resistance is (i.e., the lower w), the lower the probability that it will evolve, and if it does evolve, the lower the impact on the driving Y intervention.
Seasonal cycles
It has been shown above that a key parameter affecting the probability of resistance evolving is the initial population size, through the combined parameter uN
0. In many locations, the number of mosquitoes shows dramatic fluctuations between wet and dry seasons [33, 34]. Previous theory has shown that the probability a beneficial mutation establishes and goes to fixation can be affected by such fluctuations, and can depend upon when in a seasonal cycle the mutation arises [35, 36]. Seasonality can also affect the time to elimination by a driving Y, depending on when in the cycle the releases are made [13]. To investigate the consequences of seasonal fluctuations in mosquito numbers for the evolution of resistance to a driving Y, periodic variation is incorporated into the model of mosquito demography via a sinusoidal time dependence in the parameter γ for density dependence in the mosquito recruitment rate, such that instead of being a constant, it varies seasonally:
$${\upgamma}\left( t \right) = {\upgamma}_{0} [a]\left( {1 + a\sin \left[ {\frac{2\pi t}{T}} \right]} \right)$$
(8)
Here, a is the amplitude of the oscillations in γ(t), and T is the seasonal period, in mosquito generations in (8) to be consistent with time t, which is normalized with generation time (1/μ). An increase in γ(t) (and thus reduction in recruitment rate) during the dry season could arise from a reduction in the number or productivity of breeding sites at this time. The effects of seasonality are investigated by varying the amplitude a, and as a is varied, γ0[a] is adjusted such that the mean population size over the cycle is kept constant (Additional file 1: Section A1.1a).
As an example, Fig. 4a shows the periodically-varying female population, with amplitude a chosen to give a peak/trough ratio of female mosquito numbers of 100:1, and a representative mosquito generation time of \(1/\mu = 20\) days. Figure 4a also shows the deterministic female wild-type population for an initial amount of synthetic driving Y males, \(h_{0 }\), introduced at different times of year. In Fig. 4a, the time of the first release (going into low season) is benchmarked as year one, with other sample releases at t
release
= 3 months (population trough), 6 months (rising population) or 9 months (coincides with population peak) after year one. The female population (versus time) after introduction of the driving Y follows different paths to extinction depending on the time of year that the driving Y is released, t
release
. For no seasonality, for these parameters, it takes roughly one year (i.e. T = 18.25 mosquito generations) for the driving Y to crash the population to \(\cong 1\%\) of its initial value for typical parameters.
Depending on the point in the yearly cycle that the driving Y is released, any mutation that then arises will encounter different changing conditions (i.e., population sizes and growth rates), which will in turn affect the probability that it will establish. In Fig. 4b, P
Mut
, P
1 and \(P_{con}\) are shown as functions of the times of synthetic driving Y release, t
release
, into the wild-type population during the year. The interplay between the yearly seasonal variation and dynamics of the driving Y establishment and extinction of the population strongly influences mutant creation and survival probabilities. The probability \(P_{Mut}\) that a mutation arises (and subsequently may or may not survive) is highest when the driving Y is released at \(t_{release} \cong 7.8\) months after each yearly benchmark in (Fig. 4b), which is well into the high season, so that subsequently the peak of the driving Y male coincides with the population peak and maximizes the rate of mutant creation. But mutants that are likely to have arisen during the population surge at high season subsequently experience periods of dropping recruitment rate, and thus decreased chance of the mutation establishing, so conversely \(P_{con}\) is at its lowest for driving Y release at \(t_{release} \cong 7.8\) months after each yearly benchmark.
The combination of these effects results in P
1, the probability that at least one mutant will arise and survive, and the population will persist, being less for all release times for high seasonal variation (peak/trough ratio = 100:1) than for no seasonal variation, and at two times of the year (just before the population peak and just before the trough), it is approximately tenfold less. For low-amplitude seasonal variation (not shown), there are some individual release times in the year for which P
1 is greater than in the equivalent non-seasonal model; however, P
1 averaged over all times of driving Y release is always less than for the model with no seasonality for all amplitudes of variation (Fig. 5). Thus, comparing populations with the same mean population size, seasonal variation decreases the overall chance of successful mutation and population rescue (when averaged equally over all possible release times).
Model II: Trans-acting suppressor
Cost-free suppression
Now the possibility is considered that a mutation arises on an autosome that suppresses the expression or activity of the driving Y. As before, to analyse the fate of new mutations, a combination of a branching process model and Gillespie stochastic simulations is used (averaged over \(10^{6}\) runs) (Additional file 1: Section A3). Here mutations can arise in all individuals, not just the daughters of driving Y males, though for simplicity it is still assumed that no mutations exist before release. v (0 ≤ v ≤ 1) is defined as the chance of the suppressor mutation arising on the relevant autosome (per individual per autosome, for all births, male or female).
With baseline parameters of R
m
= 6, and m = 0.95, \(\sim 7.3\)
N
0 individuals are born between release of the driving Y and elimination of the population (see explanation after Additional file 1: Eq. A2.17); and so there are more individuals in which suppressor mutations can arise than for target site resistance mutations. However, the probability of a new suppressor mutation surviving stochastic loss is less than for a resistant mutation at all times of arising [compare p
est
(t
a
) in Additional file 1: Figures A2.1b (Model I) and A2.4b (Model II)]. A contributing factor is that the suppressor is on an autosome rather than on the X-chromosome: firstly, new suppressor mutations arise in males and females rather than only in females, and secondly, mutant males H
S
and M
S
pass the suppressor mutation equally to males and females while \(H_{R}\) and \(M_{R}\) pass on the resistant X-chromosome at the same rate but only to females. Consequently, when the mutation is rare, the proportion of time spent in a male vs female mutant is higher for the suppressor mutation than for the X-chromosome mutation, and since mutant males have a lower probability of stochastic survival than females (see Additional file 1: Figure A2.4b), early suppressor mutant populations are overall less likely to survive than resistant ones. Despite the lower probability of each mutation establishing, the greater opportunity for mutations to arise means that for equal mutation rates (i.e., u = v), a suppressor is more likely to establish than a target site resistance allele. With \(v\varvec{N}_{0}\) = 1 and for baseline parameters, \(P_{Mut} \cong 1\) and \(P_{1} \cong 0.6.\)
The effect of varying the underlying parameters \(v\varvec{N}_{0}\), R
m
, and m is now considered. The effects are much like those for the target site resistance model, and for much the same reasons, though quantitative details differ. From (6) in the limit of low uN
0 ≪ 1, P
Mut
is proportional to \(v\varvec{N}_{{0\varvec{ }}}\), \(P_{Con}\) is largely independent of \(v\varvec{N}_{0}\), and therefore P
1 is proportional to \(v\varvec{N}_{0}\). These results are shown for baseline parameters in Fig. 6a. As R
m
increases (while keeping N
0 constant), the time to eliminate the population increases, giving more opportunity for suppressor mutations to arise, and if they do arise, then they have a higher probability of establishing because low-density recruitment rates are higher, resulting in a higher overall probability of a suppressor establishing (Fig. 6b). As the transmission rate of the driving Y (m) increases, the time to elimination decreases, and therefore the opportunity for a suppressor to arise decreases. A difference from the target site resistance model is that even for m = 1, mutations can arise in the suppressor model, while for the previous model, mutations only arise in females fathered by driving Y males, and therefore none can arise for a 100% male sex bias. The conditional probability of a mutation surviving also decreases with increasing m, and therefore the overall probability P
1 also decreases (Fig. 6c).
Costly suppression
The case of the suppressor allele having a cost is now considered. For simplicity, it is assumed that heterozygous females and males have fitness w < 1, and homozygous females and males have fitness w
2.
Focussing firstly on the deterministic model, Fig. 7 is a contour plot showing the fate of a suppressor and the population as a function of w and R
m
(see Additional file 1: Section A1.2a, for calculation of equilibria from the system of deterministic equations). With a fitness cost to the suppressor mutation, one difference to the resistant model is that the region where the allele with the suppressor mutation is fixed, \(1 > w \ge w_{1} = 1 - v\), is extremely narrow for v ≪ 1 and is not discernible on Fig. 7 with baseline parameters (m = 0.95, v = 10−7). Below \(w_{1}\) (all visible areas on the plot), the suppressor allele is at intermediate equilibrium with the wild-type allele. Also shown are curves of constant percentage suppression of the total female population (dotted lines). The population is eliminated in the shaded area below the 100% extinction line w = w
ex
. There is little dependence of the results on \(v\) for v ≪ 1, and an insignificant decrease in w
1 and w
ex
with increasing Y-drive m (not shown). Thus again, unless fitness cost to the mutation is too high and/or R
m
is low, a suppressor mutation may establish and spread, leading to population rescue.
Note that trans-acting suppressors are less tolerant of fitness costs than target site resistant alleles (e.g., the shaded area is larger in Fig. 7 than in Fig. 3). This is because the resistant X-chromosome has a direct transmission advantage compared to sensitive alleles in the presence of the driving Y [passed to half rather than (1 − m) of offspring from a driving Y male]. The autosomal allele has no such transmission advantage over the wild-type, although in both cases the mutation restores the 50:50 sex ratio. This difference also fits with the already noted difference that suppressors have a lower probability of establishing than target site resistant alleles. Finally, as expected, increasing the fitness cost of a suppressor reduces the probability of it surviving stochastic loss and establishing (Fig. 6d).