A theory of resistance to multiplexed gene drive demonstrates the significant role of weakly deleterious natural genetic variation

Significance CRISPR-based gene drives have the potential for controlling natural populations of disease vectors, such as malaria-carrying mosquitoes in sub-Saharan Africa. If successful, they hold promise of significantly reducing the burden of disease and death from malaria and many other vector-borne diseases. A significant challenge to success is the evolution of resistance. Here, we develop a theory of resistance for multiplexed drive, which shows the importance of weakly deleterious naturally occurring genetic variation, whose effect is significantly amplified compared to de novo mutation. These results provide a fundamental basis to estimate how many guide RNAs are required to prevent resistance in the face of natural genetic variation.

where the mean fitness of males and females is given by the quadratic forms 8w f = x T t W f y t [S3] 9w m = y T t W mxt [S4] 10 and where the operation is the element by element (Hadamard) multiplication operation, i.e. z = x y is equivalent to 11 zi = xiyi. 12 The fitness matrix for males, W m is a 4 × 4 matrix, where each entry, (Wm)ij = 1. For females, we assume all functional females: [S5] Including mutations is straightforward using a mutation probability matrix M , whose elements Mij represent the probability Finally, we add the effect of drive by following (S1), where we assume drive acts during meiosis after viability selection.

31
Drive is only produced by heterozygotes W/D, whose frequency after selection and mutation is (1 − µ)w s 14 /ws(x1y4 + x4y1) for 32 females and males, where s = {f, m}, respectively, which gives the following allele frequency update equations [S10] where Ht = (x1)t(y4)t + (x4)t(y1)t is the drive-heterozygote frequency (before viability selection acts) and the vector κ W is a here Rm is the absolute mean number of offspring from females, and α is a parameter controlling density dependent growth, 48 which we typically tune to give a carrying capacity K = (Rm(w f +wm) − 1)α = (Rm − 1)α = N0, the initial population size we 49 set for the population in the absence of drive. In these simulations we take a typical value of Rm = 6.

50
Multiplexed drive. In general, notational complexity increases considerably when we consider 2-fold and higher orders of 51 multiplexed drive. We will only consider 2-fold and 3-fold multiplex in this paper, but we will set out a general framework for 52 any order of multiplex m. 53 Firstly, how many different haplotypes nH are there for each degree m of multiplex? For n = 3 alleles at each site (excluding 54 drive), for 2-fold multiplex, as we calculated above there are n(n + 1)/2 = 6 haplotypes, which is the number of unordered pairs, 55 corresponding to the number of upper diagonal elements for n × n matrix; including drive, we then need to track a total of 56 nH = 7 haplotype frequencies. For 3-fold multiplex, the number of haplotypes to be tracked is the number of unordered triplets, 57 which is not so trivially calculated, but amounts to calculating the analogue of the number of "upper diagonal" elements of a 58 n × n × n tensor. The answer is that there are n (m) /m! = 3 (3) /3! = 10 haplotypes, excluding drive, for m = 3 fold multiplex, 59 where n (m) = n(n + 1)(n + 2)...(n + m − 1) is the rising factorial -note this expression holds true for any m integer. So 60 including drive, for 3-fold multiplex, there are a total of nH = 11 haplotypes frequencies to be tracked. 61 We also need to decide on the ordering of the haplotypes in the frequency vectors x and y. The convention we use is column 62 ordering and its generalisation for higher dimensional tensors, which gives the mapping of indices to haplotypes as shown in 63 Table 1 and Table 2 64 The fitness matrices for 2-fold multiplex are 7 × 7, while for 3-fold they are 11 × 11. For males their fitness matrices are filled  for m = 2 the probability of mutation is where haplotype k1k2 maps to j and haplotype k 1 k 2 maps to i, using Table 1 For m = 3, equivalent considerations give the following mutation probability from the j th to i th haplotype where haplotype k1k2k3 maps to j and haplotype k 1 k 2 k 3 maps to i, using Table 2  This is formalised mathematically below and extended to any number of m of gNRAs.

97
To consider how multiplexed drive affects haplotype frequencies, we need to consider now that there is more than just a 98 single genotype that is affected by drive. All genotypes with drive on one chromosome and at least a single W allele on the 99 other can be converted to homozygous drive; for m = 2 there are three which are enumerated in Table 3 and for m = 3 there 100 are six enumerated in Table 4. site. wild type is the background colour, a white marker indicates the cleavage state, a long red bar represents drive inserted at the target site, a grey marker is an NHEJ event, which are functional (green) or non-functional (red). [S17] Note that the first column of κ is κ W , but with κ W 1 → κ W 1 + 1.

106
With this matrix the fraction of gametes that have the i th haplotype from drive heterozygote j can be calculated for any 107 value of m. For haplotypes i < nH (i.e. not including drive) this is given by where n W is the number of W alleles in the drive heterozygote j, and n W , n R , n N are the number of W, R, N alleles in haplotype 110 i, and as before the haplotypes k1k2 and k 1 k 2 map to indices using [S19]

D R A F T
Together these equations allow us to write down expressions for the update equations for haplotype frequencies between generations for m−fold multiplex [S20] where Hm is the set of indices of haplotypes that have at least a single W allele, which is a function of m the degree of RR and RRR, respectively for m = 1, 2, &3 gRNAs), none of these rise to sufficiently high frequency to achieve population 131 rescue as they are deleterious as homozygotes, and as drive heterozygotes, haplotypes containing W can still be converted to D.

132
The resistance haplotypes that have m copies of R do increase to high frequencies giving resistance and population rescue. The results below indicate that in general 1) we expect resistance mutations to arrive concurrently for NHEJ vs sequentially 143 for de novo SNPs and 2) for the same individual rate of producing functional mutants at all m sites, the critical population size 144 for de novo generation of mutants is smaller, and it is easier for resistance to arise de novo compared to NHEJ generation of 145 mutants and the disparity increases as the number of gRNAs m increases. This second effect, however, will be largely masked 146 by the generally much higher rate of generating NHEJ mutants compared to de novo SNPs. probability that none of these mutants establishes:

D R A F T
where N * n = (4s b νβQ(τ )) −1 , which gives us the result that we expect resistance to arise significantly, when N > Nn or the 169 rate of generation of NHEJ mutants is 2N νβ ∼ (2s b Q(τ )) −1 , then the probability of resistance is large.

170
For m = 2 gRNAs, a single resistance mutation (WR or RW) is not sufficient to prevent the copying of drive and so only 171 haplotypes with two functional resistance mutations (RR) have a selective advantage in the presence of drive. Assuming the 172 frequency of all mutants is small, and so assuming the fact that single mutants are selected against can be ignored (due to 173 drift), the difference equations describing this scenario are where r1 is the frequency of mutants with a single functional resistance allele at either of the two target sites and r2 is the 177 frequency of mutants that have a functional resistance allele at both target sites. If we make the same assumption that initially 178 the frequency of drive will be small, such that q(1 − q) ≈ q, then the solution for the frequency of single-fold functional 179 resistance mutants is r1(t) = 2(1 − ) νβQ(t), and so we have a difference equation for only r2: [S28] 190 and so comparing to Eqn. S22, we have the same difference equation but νβ → ( νβ) 2 and so would expect the probability of where we assume that initially q 1. The solution is simply r(t) = ξµt. Using the same simple heuristic as for NHEJ, this then 199 leads to the probability of fixation of mutants generated up to some time τ , related to the timescale for fixation of drive, to be: For m = 2 gRNAs, the difference equations for de novo generation of single and double functional resistance mutants is where again we assume the frequency of drive is initially small and we ignore any fitness effects for mutants at small frequency 209 in calculating the frequency of mutants. As above the frequency of single mutants will be r1(t) ≈ 2ξµt, where now as there 210 are two paths from the wild type to generate a single mutant there is a factor of 2, and so the difference equation for double 211 functional mutants will be 212 r2(t + 1) ≈ r2(t) + 2(ξµ) 2 t + (ξµ) 2 . [S33]

D R A F T
then be p ≈ 1 − e −N/N * d with N d = (4(ξµ) 2 s b τ 2 ) −1 . Again the arguments can be extended to m = 3 gRNAs and we find 218 N d = (4(ξµ) 2 s b τ 3 ) −1 , to leading order in τ . 219 We can fit the simulation data of the de novo probability of resistance vs N as ξ is varied, using this heuristic theory with 220 a single fitting constant γm and N * d = (4(ξµ) m γm) −1 , for each value of m and we find γ1 = 0.76 ± 0.004, γ2 = 6 ± 0.04 and 221 γ3 = 55 ± 0.3, as shown in Fig. S5; using the relation γm = s b τ m , we find these data are consistent with τ ≈ 8.5 generations 222 and a beneficial selection coefficient of s b ≈ 0.09. which has shape parameter θ = 4N ξµ and rate α. The mean of this distribution is r = θ/α = ξµ/σ, which is the classic 235 mutation-selection balance frequency. Now the probability of fixation given an initial frequency r0 and beneficial selection where the last approximation assumes α b = 4N s b 1. The average probability of fixation assuming a distribution ρ(r0) is 239 then simply calculated by averaging over the initial frequency r0, which gives where Z is the normalisation constant of the gamma distribution with population scaled fitness f : where extending the integral to ∞ will have small penalty in the strong selection limit, i.e. as long as f 1. Evaluating this 243 and re-arranging we find the result of Hermisson & Pennings in the strong selection limit: And so the probability of resistance p = 1 − e −N/N * s , where N * s = 4ξµ ln 1 + s b σ −1 , which is of the same form as for NHEJ 246 and de novo SNP mutations. There is only a single unknown parameter s b , and so fitting this form to the simulation data 247 D R A F T for the p vs N , for varying ξ we find s b ≈ 0.5 as shown in Fig. S6 (m = 1 panel), which is larger than found from de novo 248 simulations, and likely reflects some mechanistic aspect not reflected in these simple heuristic considerations.

249
For m = 2, we could try to calculate exactly the joint steady-state distribution of the frequency of all mutants, under the 250 sequential mutation scheme WW → {WR, RW} → RR, but this is in general an unsolved problem in population genetics, since 251 there is a net flux or current of probability through the states of the system (detailed balance is not obeyed) and so calculating 252 this joint distribution is not straightforward. Here, instead we appeal to a simple heuristic, which reproduces the scaling 253 with σ seen in the simulation data; we assume the distribution of the double-mutant frequency r2 follows the same gamma 254 distribution, but that mutation to this allele is solely from single mutants whose frequency is assumed fixed at its equilibrium 255 frequency r1 = 2ξµ/σ, where the factor 2 arises from the multiplicity 2 of paths to single mutants from the wild type; this is 256 reasonable given that sequential mutation rather than concurrent is dominant, as demonstrated for de novo mutation. Given 257 this and that the fitness cost of double mutants, before drive is introduced, is −2σ, the distribution of the double functional 258 resistance mutants is Note that there is an implicit assumption in all these calculations, that establishment must arise before the population is 270 removed, which will be true if ∼ 1/s b is much smaller than the timescale over which the population is removed.

271
Equilibrium frequency of multi-site deterministic mutation selection balance 272 If before the introduction of drive we assume the wild type is almost fixed and the k-fold resistance haplotype has frequency r k , 273 that σ µ and x k−1 x k then the ODE for the k th resistance haplotype is where we have ignored combinatorial factors in the rate of mutation between haplotypes with k resistance alleles and the extra 276 fitness cost as k increases. We find the fixed point for this set of couple equations recursively, starting with r * 2 = µ/σ (given the In the stochastic analysis of the previous section this becomes the mean frequency of the k-fold resistance haplotype: 282 r k = r * k .