Spatial structure undermines parasite suppression by gene drive cargo

Two kinds of engineered gene drives

The gene drives proposed and developed for genetic engineering fall into two classes. One class relies on homing, whereby the drive element cuts the genome at a specific site and inserts itself into that site (Burt, 2003; Gantz & Bier, 2015; Kyrou et al., 2018). A homing drive’s fitness advantage comes from a transmission bias in gametes of heterozygotes. CRISPR technology has greatly facilitated this type of engineering because CRISPR-Cas9 is a site-directed nuclease. The other class relies on biased offspring survival, many of which are known as ‘killer-rescue’ systems (Chen et al., 2007; Gould et al., 2008; Marshall & Hay, 2011; Legros et al., 2013; Akbari et al., 2013; Akbari et al., 2014; Buchman et al., 2018; Oberhofer, Ivy & Hay, 2019; Champer et al., 2019a). One of the major differences between these two classes of drive elements is the speed and ease with which they spread. A homing element spreads rapidly and can, in theory, be successfully introduced with a single individual. Killer-rescue systems spread more slowly and often must be introduced above a threshold density to spread, although that distinction is not absolute (Champer et al., 2019a).

Mass action dynamics

Gene drives have traditionally been modeled and understood in the context of well-mixed populations (e.g., Prout, 1953; Bruck, 1957; Hamilton, 1967; Burt, 2003; Marshall & Hay, 2011; Legros et al., 2013; Akbari et al., 2014; Unckless et al., 2015; Beaghton et al., 2017; Godfray, North & Burt, 2017). A homing gene drive gains its advantage from heterozygotes, the non-drive allele of a heterozygote being replaced with the drive allele during reproduction. Heterozygote frequency is enhanced with outcrossing (mass action), depressed with inbreeding and some other types of assortative mating. Even killer-rescue systems presumably benefit from mass action, to distribute the killing beyond those inheriting the rescue. With mass action, an efficient homing drive can spread from low frequencies to near-fixation in close to 10 generations (Burt, 2003; Godfray, North & Burt, 2017; Beaghton et al., 2017).

The evolution of a gene drive and its associated cargo can be divided into two phases. The first phase encompasses the short-lived spread of the drive. Although gene drives are potentially highly efficient, various types of fitness effects, imperfections in the drive mechanics and variation in the host population can limit the final coverage of the drive (Deredec, Burt & Godfray, 2008; Godfray, North & Burt, 2017; Beaghton et al., 2017; Champer et al., 2017). Once the drive has spread to its limit, phase two sets in, whereby evolution proceeds according to fitness effects on the host. Any fitness cost stemming from the drive allele or its genetic cargo now begins to select a population reversal toward loss of the cargo and/or drive, favoring alleles resistant to the drive or cargo-free drive states. The speed of this reversal depends heavily on fitness costs and on the initial frequencies of the different parties; it is typically much slower than the spread of the drive (Beaghton et al., 2017). In the long term, a genetic cargo with any fitness cost will be lost. The social benefits of the cargo must therefore be manifest in a time frame compatible with its expected duration.

Population spatial structure

Gene drives require reproduction. Their spread will thus follow the conduits of reproductive connections in the host population, which may well have a strong spatial component—as when individuals mate with neighbors (North et al., 2013; Beaghton, Beaghton & Burt, 2016; Tanaka, Stone & Nelson, 2017). Any genetic variation that arises in the gene drive or cargo, such as mutations that delete or down-regulate the cargo, will be propagated along those conduits and expand accordingly, leading to spatial structure in parasite inhibition (Beaghton et al., 2017). Even more simply, for purely dynamical reasons, the drive may fail to reach isolated pockets of the population (North et al., 2013). In turn, spatial structure of a genetically variable inhibitor will often mean that different locations of the parasite experience different levels of inhibition. With spatial structure, even small regions of reduced parasite suppression may enable parasite persistence which then facilitates parasite evolution of resistance to the inhibitor.

Pre-existing genetic variation in the host population may also affect gene drive efficacy, spread and cargo expression (e.g., Drury et al., 2017; Champer et al., 2019b). For example, some designs for a harmless homing drive have it target a genomic region that can be disrupted with little or no fitness effects; such a region may thus not be strongly selected to conserve sequences and may be variable across the host population, blocking gene drive spread in some regions. (One design avoids this problem by targeting an essential gene and carrying a cargo that replaces the targeted gene (Burt, 2003; Champer et al., 2019c)). Cargo gene expression may likewise be affected by the genome in which it resides, and geographic variation in genomic content may lead to geographic variation in cargo expression.

Our intent is to investigate the consequences of spatially structured inhibition of the parasite/pathogen. The details of structure will typically be implementation-specific, but an appreciation for the importance of spatial structure when it exists may be a requisite for successful application of a gene drive cargo.

Mathematical models

Our model is most easily applied to an asexual pest/parasite infecting a single host species. Although not specifically modeled, our problem may be extended in spirit to a parasite transmitted between two host species, as to Plasmodium transmitted by a mosquito to humans and back to mosquitoes; in this case, the gene drive is introduced into the mosquito to block Plasmodium reproduction and transmission. Our models merely omit the second host, but we conjecture that the effects they reveal apply to that case, subject to some conditions mentioned in the Discussion.

The social goal is to suppress parasite reproduction with a genetic cargo in the host. To keep the problem simple, we assume that a gene drive and its cargo have already swept through the host species. (The actual process of gene drive evolution is thus ignored, and indeed, non-gene-drive methods of cargo infusion may also be used to achieve this end (Okamoto et al., 2014)). In any one host individual, parasite inhibition by the cargo occurs in one of three states: (i) full inhibition, (ii) partial inhibition, or (iii) no inhibition. Partial inhibition would result from weak expression of the cargo in the host; no inhibition would result from loss of the cargo from the gene drive or resistance to the gene drive itself, such that the host individual lacks the gene drive and its cargo altogether. The formulation of the model is trivially extended to multiple states of partial inhibition.

The model counts numbers of parasites in each type of host. Host type merely translates into parasite reproduction. Notation is

• The relative frequency of hosts of type h

• The fecundity of a parasite genotype g in hosts of type h

• The current number of parasites of genotype g in hosts of type h

• The current number of progeny produced by genotype g across all patches

• The fidelity of parasite reproduction to hosts of the same type,

where parasite genotype g ∈ {0, 1, …, M} and host type h ∈ {0, 1, …, H}. Here, host type 0 indicates hosts with no cargo (hence no parasite suppression), and larger values of h correspond to increased levels of suppression, with H denoting the number of types of cargo-carrying hosts differing in some aspect of parasite suppression; parasite genotype 0 is the wild type (with no protection against the cargo), and larger values of g correspond to mutant strains with increased levels of resistance to the cargo.

To approximate the separation of phases between rapid gene drive spread and the subsequent effect of cargo on the parasite, we let the x_h and b_gh be constant in time, so the only changes are in parasite numbers. Time is discrete. For biological reasons, in the case of three host types, we also assume that the fecundity of parasite genotype 0 satisfies b₀₀ > 1 > b₀₁ > b₀₂ , which means that the parasite has negative growth in all host types except the one lacking cargo.

To set the stage for a structured host population, we suppose that hosts are clustered in patches of similar host types, patch types designated by subscript h (‘Host’ and ‘patch’ are used interchangeably below, but ‘patch’ helps convey structure). Patches could result, for example, from limited host migration and gene drive spread through the host population in a manner that follows host population structure. The clustering of hosts and the consequent movement of parasites between patches determines the extent to which structure is experienced by the parasites. To establish a mass-action baseline, adult parasites reproduce and release all progeny into a random pool, from which they settle into each of the H + 1 patch types at frequencies x₀, x₁, …, x_H.

Spatial structure is modeled indirectly by assuming that a fraction α of the progeny born in a patch type remains in the same patch type without entering the random pool; this ‘fidelity’ increases with the retention of progeny in their natal patch type. This process is fundamentally the same as migration in standard population genetics problems (Crow & Kimura, 1970). Our formulation is different in that α denotes a lack of movement from the natal site instead of movement between patches/populations; this formulation leads naturally to calculating the null case of mass action (α = 0), which is presumably the default expectation in gene drive applications. Note that fidelity to a patch type is imposed without inbreeding, so the numbers within a patch are assumed large enough that consanguinity can be ignored.

We first consider a simple system of linear dynamical equations. The overall progeny output of strain g across all environments is (1)${\displaystyle N_g=\sum _{h=0}^H{n}_{gh}{b}_{gh}.}$

Using primes to indicate one generation hence, (2)${\displaystyle n_{gh}^{\prime }=\alpha n_{gh}{b}_{gh}+\left(1-\alpha \right)N_g{x}_h.}$

The joint dynamics for genotype g across all host types can be written as a matrix projection recursion (Caswell, 2006). Dropping the genotype subscript g for ease of visualization, the recursion has form (3)${\displaystyle {\mathbf{n}}^{\prime }=\mathbf{M}\mathbf{n}}$ where n and n′ are (H + 1)-dimensional column vectors with hth element n_h and $n_h^{\prime }$, respectively. The (H + 1) × (H + 1) matrix M has diagonal elements M_hh = b_h[α + (1 − α)x_h] and off-diagonal elements M_hi = b_i(1 − α)x_h, i ≠ h. In words, matrix element M_hi describes the rate that individuals who originate in patch type i contribute to the abundance in patch type h at the next time step. The densities of the genotype g at any time t, n(t), can be computed as n(t) = M^tn(0) (Caswell, 2006).

This model allows for a biological anomaly: when a genotype i is initially assigned to a patch (j), a sufficiently large combination of fecundity and fidelity (b_ijα > 1) allows its numbers to persist even when the patch type is absent (even when x_j = 0). This effect occurs because, once a genotype exists in a patch (by initial conditions), a portion of its growth comes from offspring who stay in the patch to reproduce. This effect is independent of patch size. Because we interpret x_j as the fraction of hosts of type j, we require x_j > 0 for any host type that harbors parasites. This requirement is further imposed when enforcing a carrying capacity (see below).

These equations assume fixed fecundities and thus allow unlimited population growth. They should be adequate to decide whether parasites persist or die out, because they can be applied to deterministic dynamics at low population densities to describe the direction of population growth, when fecundities are intrinsic and unaffected by densities. For dynamics and evolution in high-density populations, we use a related model that imposes density regulation. For simplicity, we limit our description of density regulation to the case of two patches and two strains.

Parasite persistence is facilitated by spatial structure

The growth or suppression of the parasite depends on whether the magnitude of the leading eigenvalue of M in Eq. (3) exceeds 1. Density dependence can be ignored when addressing persistence, but an eigenvalue exceeding 1 merely indicates that parasites have positive growth somewhere in the environment, perhaps in only a tiny locale, with negative growth everywhere else. Although the characteristic equation is easily found, the leading eigenvalue (λ_max) of a genotype is tractable for arbitrary H only for the extremes of mass action (α = 0), and complete parasite isolation among host types (α = 1). For mass action, (6)${\displaystyle {\lambda }_{\max ,g}=\sum _h{x}_h{b}_{gh}.}$

As is well appreciated for mass action, parasite growth is just the weighted average fecundity across all host types. Small levels of weak suppression (i.e., low values of x₀ despite possibly high values of b₀₀) will not themselves enable parasite persistence except when parasite fecundity is extraordinarily high. When more than one genotype has an eigenvalue greater than 1, density dependence will determine which one prevails (see below).

With complete parasite separation across the host types (α = 1), (7)${\displaystyle {\lambda }_{\max ,g}=b_{g0}}$ (λ_max is associated with patch type 0 because cargo-free hosts are assumed to offer the highest fecundity of all patches, regardless of genotype). Here, the parasites inhabiting each host type have their own eigenvalue, and any host type with b_gh > 1 will allow parasites of genotype g to persist in that patch type (subject to competition among different parasite genotypes). In this extreme, the values of x_h no longer matter: even a small fraction of permissive hosts will allow the parasite to persist.

The question motivating our study is how sensitive parasite persistence is to fidelity α—reflecting spatial structure. As the eigenvalue for several patch/host types is unwieldy, we reduce the problem to just two patch types, h = 0 and 1 (fully permissive and fully blocking of wild-type); this reduction also simplifies patch type abundances: x₀ + x₁ = 1. For this case, the largest eigenvalue is

(8)${\displaystyle {\lambda }_{\max ,g}=\frac{1}{2}\left\{\alpha \left[b_{g1}{x}_0+b_{g0}\left(1-x_0\right)\right]+b_{g1}\left(1-x_0\right)+b_{g0}{x}_0\right\}+}$ ${\displaystyle \frac{1}{2}\sqrt{{\left[b_{g1}\left(\left(1-\alpha \right)\left(1-x_0\right)+\alpha \right)+b_{g0}\left(\left(1-\alpha \right)x_0+\alpha \right)\right]}^2-4\alpha b_{g1}{b}_{g0}\left[\left(1-\alpha \right)\left(1-x_0\right)+\left(1-\alpha \right)x_0+\alpha \right]}.}$

Values of the eigenvalues for different fecundities in the two patch types are shown in Fig. 1, Figs. 1A–1D differing in fidelity (α) and cargo-free patch size (x₀) values. The graphs show four contour lines, but λ_max = 1—the boundary between persistence and extinction—is thicker than the others. It is clear that persistence is enhanced by spatial structure, though typically a large cargo-free fecundity is required if the cargo is effective (b₀₀ must be well above 1 when α and b₀₁ are small). But for small cargo-free patch sizes (small x₀), parasite persistence becomes possible despite low fecundity values in cargo-free hosts when the cargo becomes moderately ineffective (when b₀₁ exceeds 1). In addition to the effect of fecundity, there are also effects of α and x₀; one interesting effect is that the isoclines are visually step-like except in Fig. 1A. For those cases, fecundity in cargo-bearing hosts has little effect on the parasite growth rate until b₀₁ exceeds λ_max.

Most empirical interest is likely to be in the extreme case that the cargo completely suppresses wild-type parasite reproduction (b₀₁ = 0), as that would be the goal of the engineer. Indeed, any gene drive release could be avoided until such an appropriate inhibitor was found. This case corresponds to the the sliver defined by the vertical axes in Fig. 1. Furthermore, this case is highly tractable: (9)${\displaystyle {\lambda }_{\max ,0}{\mid }_{b_{01}=0}\equiv {\lambda }_0=b_{00}\left[x_0+\alpha \left(1-x_0\right)\right].}$

Persistence in this case requires ${\lambda }_0=b_{00}\left[x_0+\alpha \left(1-x_0\right)\right]\ge 1$. This implies, for example, that persistence of the parasite is assured even when completely inhibited by the cargo as long as fecundity in cargo-free hosts (b₀₀), exceeds ${\left[x_0+\alpha \left(1-x_0\right)\right]}^{-1}$. Consideration of this minimum cargo-free fecundity (Fig. 2) shows that spatial structuring with fidelity (α) well below 1 (e.g., 0.5) enables parasite persistence for even rare cargo-free hosts (small x₀), as long as the parasite can grow there moderately well (e.g., b₀₀ > 2). The results are largely insensitive to patch size x₀ when fidelity reaches 0.6.

It is also easy to see that, when cargo-free patch size x₀ is small, as it should be if the engineering worked as expected, the effect of increasing fidelity α is approximately the same as increasing the frequency of permissive hosts (x₀)—both are mostly linear effects and of the same magnitude. Whereas x₀ is somewhat under the control of the engineer, α is not and has the potential to thwart parasite eradication. Thus, a very small x₀ could easily cause parasite population collapse in the mass action case, but b₀₀α > 1 is sufficient for persistence no matter how small x₀. Persistence could be achieved by a cargo-free patch merely large enough that parasite progeny often did not disperse beyond the patch edges.

Ease of resistance evolution depends on ecology

At the population level, the evolution of resistance to a genetic cargo has many parallels with the evolution of resistance to antibiotics and pesticides. The latter problems are thoroughly studied, and it is well appreciated that resistance is especially prone to evolve under intermediate levels of drug/chemical application (e.g., Gould & MacKenzie, 2002; Andersson & Hughes, 2012; Tabashnik & Gould, 2012; Neve et al., 2014; Gould, Brown & Kuzma, 2018). Inhibition by a gene drive cargo is different from chemicals in that the levels of inhibition are established at fixed, semi-permanent levels in the near term. They are also largely unchangeable, at least in the short term, should it be discovered that they are inadequate.

Of the many factors to consider, an important one is the mutational spectrum of resistance: a cargo for which simple, single mutations can allow parasite persistence seems doomed to fail, and intuition suffices for preliminary understanding, at least deterministically. Our interest instead lies in gradual evolution and the selection of weak resistance mutations. It might be hoped that cargo-based inhibitors can be found for which resistance mutations are impossible, but a more realistic hope is that inhibitors could be found for which resistance can evolve only gradually.

The evolution of resistance can be considered in two contexts. One is known as ‘evolutionary rescue’, whereby the population is in decline and a resistance mutation potentially reverses the decline (Martin et al., 2013; Uecker, Otto & Hermisson, 2014; Hufbauer et al., 2015; Gomulkiewicz, Krone & Remien, 2017). The other context, the one addressed here, is resistance evolution in a persisting population. For our application, we imagine that the parasite is persisting because of spatial structure and would go extinct under mass action—the large majority of hosts inhibit parasite reproduction because of the cargo. We address how selection acts on a weakly resistant mutation, a mutation that is not necessarily sufficient to provide positive parasite growth from inhibitory hosts alone.

It might seem valid to evaluate long term resistance evolution from a comparison of eigenvalues of wild-type and mutant growth, in which case the preceding figures could be used to infer evolution of alternative genotypes. However, we are considering resistance evolution in established populations at which density dependence is operating. When density dependence operates locally, as assumed here, it will have a different effect on parasite growth in cargo-free patches than in inhibitory patches. If the wild-type parasite cannot grow in inhibitory patches but the mutant can, the fecundities of both genotypes will be suppressed by density dependence in cargo-free patches but the mutant’s fecundity in the inhibitory environment will be unaffected (at least while rare). Eigenvalue calculations do not include these fecundity modifications. Our model of resistance evolution uses the system of equations in (3) but with the density-dependent carrying capacity enforced as in Eqs. (4) and (5).

This model was evaluated numerically for different combinations of fidelity (α), cargo-free patch size (x₀), and genotype fecundities (b_ij). We focused on the case of a wild-type (starting genotype) that was unable to grow in the inhibitory environment of cargo-bearing hosts (b₀₁ = 0), and for which the mutant could grow in the inhibitory environment at some reduction in its ability to grow in the cargo-free environment (i.e., a trade-off was imposed between growth in the two patch types). The patterns across trials are qualitatively similar and easily comprehended (Fig. 3). Resistance could invariably evolve if it did not incur too much of a cost to growth in the cargo-free hosts.

Some ‘ecological’ patterns are evident. (1) A large effect of α—host patch fidelity—on evolution exists in some parameter ranges. Thus, to displace wild-type, the mutant is more sensitive to a reduction in cargo-free fecundity at larger α—until fecundity in cargo-bearing hosts is high enough that the mutant could sustain itself in those hosts alone (until b₁₁ exceeded 1). This effect is evident when comparing lower portions of Figs. 3B and 3D and was seen in other trials (not shown). We interpret this effect as that higher values of α increasingly partition growth in the two patch types, so that any fitness loss in the permissive patch—the one sustaining the population—is increasingly penalized, but only to the point that the mutant parasite can maintain itself in cargo-bearing hosts.

(2) Patch size x₀ also appears to have a large effect, but in the opposite direction as that of α: as x₀ increases, resistance evolves more easily (i.e., tolerates larger reductions in cargo-free fecundity, b₀₀ − b₁₀). This effect is seen by comparing Figs. 3A–3C which vary only x₀. Although not shown, the trend does not continue at high values of x₀, and the effect reverses. Indeed, mutants whose b₁₁ > 1 are not favored at high x₀ values if they suffer too much cost in b₁₀. We do not pretend to grasp these patterns intuitively. One challenge in understanding these outcomes is that x₀ has effects at different steps of the life cycle: with higher x₀, the cargo-free patch has an increased carrying capacity and thus produces more progeny; some of these progeny return to patch 0 and the others go to the random pool. But the increase in x₀ reduces the number of random-pool progeny that land in patch 1 and thereby reap the benefits of resistance.