Model framework
We use the Mosquito Gene Drive Explorer (MGDrivE) modeling framework [26] to model the spread of translocations and UDMEL through spatially structured mosquito populations (Fig. 1). This is a genetic and spatial extension of the lumped age-class model of mosquito ecology [27] modified and applied by Deredec et al. [28] and Marshall et al. [29] to the spread of homing gene drive systems. The framework incorporates the egg, larval, pupal, and adult life stages, with egg genotypes being determined by maternal and paternal genotypes and the allelic inheritance pattern of the gene drive system. Spatial dynamics are accommodated through a metapopulation structure in which lumped age-class models run in parallel and migrants are exchanged between populations according to a zero-inflated exponential dispersal kernel with parameters defined in Additional file 1: Table S1 [22, 30,31,32,33,34,35,36,37,38]. Further details of the framework are described in the “Methods” section.
Applying the MGDrivE modeling framework to our research questions, we incorporate the inheritance patterns of reciprocal chromosomal translocations and UDMEL into the inheritance module of the model (Fig. 1a, b, Additional file 2: Fig. S1), the life cycle parameters of Aedes aegypti (Additional file 1: Table S1) into the life history module, and the distribution of households in Yorkeys Knob (923 households) and Trinity Park (1301 households) along with their expected mosquito population sizes and movement rates between them into the landscape module (Fig. 1c). The suburb of Trinity Park served as a control site for field releases of Wolbachia-infected mosquitoes, to quantify the extent to which the Wolbachia infection could spread from one community to another, and plays a similar role for our simulated releases of threshold-dependent gene drive systems.
The inheritance patterns that result from chromosomal translocations are depicted in Fig. 1a. Chromosomal translocations result from the mutual exchange between terminal segments of two non-homologous chromosomes. When translocation heterozygotes mate, several crosses result in unbalanced genotypes and hence unviable offspring, resulting in a heterozygote reproductive disadvantage. This results in bistable, threshold-dependent population dynamics, confirmed in laboratory drive experiments [14]. The inheritance patterns produced by the UDMEL system are depicted in Fig. 1b. This system consists of two unlinked constructs, each possessing a maternally expressed toxin active during oogenesis and a zygotically active antidote expressed by the opposite construct. Offspring are more likely to have the antidote(s) to the maternal toxin(s) when transgenes are present at moderate-to-high population frequency. This produces threshold-dependent dynamics since, above a critical frequency, the selective advantage of the antidotes outweighs the selective disadvantage of the toxins, and below the critical frequency, the selective disadvantage of the toxins dominates. Mathematical models predict a threshold frequency of ~ 24% in the absence of an additional fitness cost—a result consistent with laboratory drive experiments for a construct having a modest additional fitness cost [15].
Model validation
Using data from field trials of Wolbachia-infected Ae. aegypti mosquitoes in Yorkeys Knob and Gordonvale, Australia [21], we validated our modeling framework prior to application to other threshold-dependent drive systems. Wolbachia biases the offspring ratio in favor of those carrying Wolbachia through a mechanism known as cytoplasmic incompatibility, in which offspring of matings between infected males and uninfected females result in the death of some or all progeny, while matings involving infected females tend to produce infected offspring [39]. For the wMel strain of Wolbachia that was used in the Australian field trials, incompatible crosses produce no viable offspring, and Wolbachia is inherited by all offspring of infected females. There is also a fitness cost associated with Wolbachia infection, the value of which has been estimated between 0 and 20% [21, 40, 41].
We used Wolbachia surveillance data from Fig. 1 of Hoffmann et al. [21] to validate our model framework and parameter values. Based on the description of the field trials in Hoffmann et al. [21], we simulated weekly releases of 20 Wolbachia-infected mosquitoes (10 female and 10 male) per household at a coverage of 30% over 10 weeks in the spatially explicit landscapes of Yorkeys Knob and Gordonvale with the exception that in Gordonvale, the fifth release was postponed by a week due to a tropical cyclone. Model predictions were calculated for a variety of literature-based values of adult mortality rate and Wolbachia-associated fitness cost [21, 33, 40,41,42] and were compared to observed Wolbachia prevalence over time. We found that model predictions most closely matched field data for a baseline adult mortality rate of 0.090 per day [33] and that predictions matched field data quite well for both 5% and 10% fitness costs, with a 10% fitness cost being closer to that estimated elsewhere [21, 40, 41] (Additional file 3: Fig. S2). Model predictions in Fig. 2 use these parameter values, along with others listed in Additional file 1: Table S1, and their agreement with the observed field data provides good validation of our modeling framework.
Population replacement and remediation for translocations
The use of translocations for transforming pest populations was initially suggested by Serebrovskii [43] and later Curtis [44] for the introduction of disease-refractory genes into mosquito populations. A number of models have been proposed to describe their spread through randomly mating populations [14, 16, 45, 46]; however, with one recent exception addressing spatial structure [18], these have largely ignored insect life history and mating structure. Such models suggest that the translocation need only exceed a population frequency of 50%, in the absence of a fitness cost associated with the translocation, to spread to fixation in a population, which could conceivably be achieved through a single seeding release round. Here, we find that incorporating life history and population structure into mosquito population dynamic models significantly increases release requirements.
In Figs. 3 and 4, based on the precedent set by the 2011 Wolbcahia field trial [21], we consider weekly releases of 20 adult Ae. aegypti males homozygous for the translocation per household for given durations and coverage levels, where coverage level is the proportion of households that receive the releases. Releases are simulated in the community of Yorkeys Knob, in which prior releases of Wolbachia-infected mosquitoes suggested a local population of ~ 15 adult Ae. aegypti per household [22], and for mosquito movement rates inferred from previous studies [25, 37, 38] (Additional file 1: Table S1). For a coverage level of 100%, and in the absence of a fitness cost, four weekly releases of 20 Ae. aegypti males (~ 3:1 released to local males) are required for the translocation to spread to fixation in the community (Fig. 3), as opposed to the single release expected when ignoring life history and population structure [45]. As coverage is reduced to 50%, the required number of releases increases to 7, and for a coverage level of 25%, as seen for the World Mosquito Program in Yorkeys Knob, the required number of releases increases to 16 (Figs. 3 and 4). Although large, these releases are achievable, considering the much larger releases conducted for sterile insect programs [47].
To simulate remediation of a translocation, we consider weekly releases of 20 adult Ae. aegypti wild-type males in the community of Yorkeys Knob, whereby the translocation has already reached fixation in that community. In the absence of a fitness cost associated with the translocation, translocations are symmetrical in their threshold dynamics, and so, for a coverage level of 100%, four weekly releases are required for the translocation to be completely remediated from the community, and for a coverage of 25%, 16 weekly releases are required for the translocation to be completely remediated (Figs. 3 and 4). Encouraging features of these results are that (i) remediation can be achieved through releases of non-biting, non-disease-transmitting males; (ii) release sizes are achievable; and (iii) despite the spatial household structure, both replacement and remediation are complete within the community. The time to replacement is highly dependent on the coverage level and number of releases, but is reasonably quick given sufficient releases. At a coverage of 50%, 20 weekly releases led to the translocation spreading to a frequency > 99% within half a year of the final release (or within 300 days of the first release). For equivalent wild-type releases, this is the same as the time to > 99% elimination.
Population replacement and remediation for UDMEL
UDMEL was the first synthetic gene drive system to be engineered that displays threshold-dependent dynamics [15]. The system consists of two unlinked constructs, each possessing a maternally expressed toxin active during oogenesis and a zygotically active antidote expressed by the opposite construct (Fig. 1b). At low population frequencies, the maternal toxin confers a significant selective disadvantage, leading to elimination, while at high population frequencies, the zygotic antidote confers a selective benefit in the context of a prevalent toxin, leading to fixation. The dynamics of this system in randomly mating populations have been characterized by Akbari et al. [15], suggesting that the system need only exceed a population frequency of ~ 24%, in the absence of a fitness cost, to spread to fixation, while the wild-type must exceed a population frequency of ~ 76% to eliminate the construct. Both replacement and remediation should therefore be achievable with 1–2 releases of transgenic and wild-type organisms, respectively; however, as for translocations, we find that incorporating life history and population structure into our models increases release requirements in both cases.
In Figs. 3 and 4, we consider weekly releases of 20 adult Ae. aegypti males homozygous at both loci for the UDMEL system in the community of Yorkeys Knob. The lower threshold for UDMEL as compared to translocations means that replacement is much easier to achieve for UDMEL. For a coverage level of 50% or higher, and in the absence of a fitness cost, a single release of 20 Ae. aegypti males leads to the UDMEL system spreading to fixation throughout the community (Fig. 3). As coverage is reduced to 25%, the required number of releases to achieve fixation increases to two (Figs. 3 and 4). As for translocations, the time to replacement is highly dependent on the coverage level and number of releases. From Fig. 3, it is apparent that UDMEL reaches total allele fixation slowly, although the number of individuals having at least one copy of the transgene increases quickly. At a coverage of 50%, 10 weekly releases lead to wild-type individuals falling to a frequency < 2% within 2.6 years of the final release (or within 3 years of the first release).
Remediation, however, is more difficult to achieve for UDMEL compared to translocations due to the higher threshold that wild-type organisms must exceed to eliminate UDMEL. Additionally, wild-type females must be included in the releases to propagate the wild-type allele because, assuming continued functioning of UDMEL components, the maternal toxins of females having UDMEL at both loci kill all offspring that do not inherit UDMEL at both loci. To simulate remediation, we first consider weekly releases of 10 adult Ae. aegypti wild-type females and 10 adult males in the community of Yorkeys Knob. In the absence of a fitness cost associated with the UDMEL construct, and for a coverage level of 75%, nine weekly releases are required for a reduction in UDMEL allele frequency over the first year (Fig. 3); however, a closer inspection of the simulation results reveals that complete remediation of UDMEL from the community is not possible even with 20 releases, as the UDMEL allele frequency bounces back. Comparison of these results to those for a panmictic population with a population size equal to that of Yorkeys Knob reveals that complete remediation of UDMEL can occur at coverage levels as low as 25% (for 16 or more weekly releases) (Additional file 4: Fig. S3). Inspection of the spatially explicit simulation results suggests that the rebound in UDMEL allele frequency in the structured population is due to UDMEL remaining at super-threshold levels after the wild-type releases in a small number of households, and slowly recolonizing the landscape following that. Complete remediation of UDMEL is possible, however, for 15 or more releases at a coverage level of 100% (Fig. 4). These results make a strong case for translocations as preferred systems to introduce transgenes in a local and reversible way as (i) remediation of UDMEL requires releases of biting, vector-competent females and (ii) release requirements for these biting, vector-competent females are burdensomely high due to the high threshold that must be surpassed consistently throughout a spatially structured population.
Confinement of translocations and UDMEL to release site
Confinement of translocations and UDMEL to partially isolated populations has previously been modeled by Marshall and Hay [16] and Akbari et al. [15]. In both cases, two randomly mating populations were modeled that exchange migrants at given rates. Population structure was otherwise ignored, as was mosquito life history. Results from these analyses suggest that translocations would spread and remain confined to populations for which the migration rate is less than ~ 5.8% per mosquito per generation [16], and that UDMEL would remain confined to populations for which the migration rate is less than ~ 1.6% per mosquito per generation [15]. These migration rates are relatively low; however, this may be beneficial for the types of landscapes we are considering here, whereby the system may spread between neighboring households, but not from one suburb to another. Recently, Champer et al. [18] showed that translocations would remain confined to and persist in a population connected to another by a “migration corridor” under a range of parameter values.
For our landscape of interest—the suburbs of Yorkeys Knob and Trinity Park—it is very unlikely that Ae. aegypti mosquitoes will travel from one suburb to another by their own flight. Extrapolating the exponential dispersal kernel used in our simulations, fitted to data from mark-release-recapture experiments collated by Guerra et al., [38] suggests these events to be negligible, before accounting for the fact that the intervening vegetated area may serve as a barrier to Ae. aegypti flight [36]. Furthermore, rare migrant mosquitoes are unlikely to cause the threshold frequency for either drive system to be exceeded, thus making spatial spread due to such movements unlikely. In considering confineability to the release suburb, we therefore model “batch migration,” in which several mosquitoes are carried, perhaps by a vehicle, from one community to another at once. Batch migration events could be thought of as several adult mosquitoes being carried at once, or perhaps more likely, as a larval breeding site, such as a tyre, being carried from one household to another, with several adults emerging from the tyre following transport. We model batch migration events as occurring between randomly chosen households, and vary the number of daily migration events and the effective number of adults carried per event. For computational simplicity, we focus on migration events from Yorkeys Knob, in which either system has already reached fixation, to households in Trinity Park, which is initially fixed for wild-type mosquitoes.
In Fig. 4e, f, we see that both the number and size of daily batch migration events affect the chance of either system establishing itself in the neighboring suburb, Trinity Park. For translocations, ~ 16 daily migration events of batches of 5 adults are required for spread in Trinity Park. For batches of 10 adults, ~ 9 daily migration events are required, and for batches of 20 adults, ~ 5 daily migration events are required. For UDMEL, ~ 3 daily migration events of batches of 5 adults are required for spread in Trinity Park, and for batches of 10 adults, ~ 2 daily migration events are required.
These results continue to make a strong case for translocations as preferred systems to introduce transgenes in a local and reversible way as (i) many more batch migration events are required to lead to spread for translocations as opposed to UDMEL and (ii) the rate of migration events required for translocations to spread is higher than what would be expected between these communities. Specifically, Wolbachia releases in Yorkeys Knob in 2011 provide evidence for occasional batch migrations to the nearby suburb of Holloways Beach; however, the spatio-temporal pattern of Wolbachia spread, as inferred from monitored trap data, suggests only ~ 1–2 batch migration events over the course of a month, consisting of less than 5 adult females per event [21].
Sensitivity analysis
A theoretical study by Khamis et al. [48] on toxin-antidote-based underdominance gene drive systems, similar to UDMEL but for which the toxins are zygotic rather than maternal [20], found that the gene drive threshold frequency is highly sensitive to (i) the increase in adult mortality rate in organisms having the transgene, (ii) the duration of the larval life stage, and (iii) the parameters determining the character or strength of larval density dependence. In Fig. 5 and Additional file 5: Fig. S4, we explore the sensitivity of our model outcomes of replacement, remediation, and confinement for translocations and UDMEL as we vary (i) the duration of the larval life stage, (ii) the baseline adult mortality rate, (iii) the fitness cost associated with the gene drive system, and (iv) the mean adult dispersal distance. For translocations, we model a 10% fitness cost as a 10% reduction in mean adult lifespan for organisms homozygous for the translocation and a 5% reduction for organisms heterozygous for the translocation. For UDMEL, since its inheritance bias is induced through the action of maternal toxins, we model a 10% fitness cost as a 10% reduction in female fecundity for organisms homozygous for UDMEL at both loci, with 2.5% additive fitness costs contributed by each transgenic allele.
For translocations, the associated fitness cost had the greatest impact on the release scheme required for the system to be fixed or remediated from the population, given the life parameters considered (Fig. 5). A 10% fitness cost led to ~ 10 weekly releases at a coverage of 50% being required for the translocation to reach fixation (an increase of 3 releases), while a 20% fitness cost led to ~ 13 weekly releases being required (an increase of 6 releases). Small changes in the duration of the larval life stage had minor impacts on the release requirements, with an increase in larval lifespan of 2 days leading to one more weekly release being required for the translocation to reach fixation, and vice versa. A 2% change in the baseline adult mortality rate and 50% change in the mean migration distance had negligible impact on release requirements. Remediation, on the other hand, requires fewer wild-type releases when there is a fitness cost associated with the translocation. A 10% fitness cost led to 5 weekly releases at a coverage of 50% being sufficient to eliminate the translocation (a decrease of 2 releases), and a 20% fitness cost led to 4 weekly releases at a coverage of 50% being sufficient for elimination (a decrease of 3 releases). Small changes in the duration of the larval life stage had minor impacts on the wild-type release requirements for elimination, with an increase in larval lifespan of 2 days or a 2% decrease in the adult mortality rate leading to one fewer release being required.
The sensitivity of our predictions regarding confinement to the release site is of particular interest, as invasion of a neighboring community may be more likely under some parameter values than others. Fortunately, a fitness cost associated with the translocation leads to a higher threshold and hence more batch migration events required for invasion of a neighboring community. A 10% fitness cost led to ~ 2–3 additional daily migration events of 10 adults required for spread to Trinity Park, and a 20% fitness cost led to ~ 6–7 additional daily migration events required (Fig. 5). Also noteworthy, a 2% increase in the adult mortality rate led to ~ 2 fewer daily migration events required for spread to Trinity Park—i.e., ~ 7 migration events for batches of 10 adults and ~ 14 migration events for batches of 5 adults. While still above inferred batch migration rates, this highlights that there could exist parameter sets beyond those explored for which invasion is feasible.
UDMEL displays similar parameter sensitivities regarding fixation and batch migration outcomes as for translocations, with the exception that these outcomes are less sensitive to fitness costs (Additional file 5: Fig. S4), likely due to the fact that fitness is accommodated through a reduction in female fecundity rather than an increase in adult mortality. A 20% fitness cost led to ~ 1 additional weekly release being required for the system to spread to fixation, whether at a coverage of 25% or 50%. Similarly, for an invasion of Trinity Park, a 10% fitness cost required ~ 1 additional daily migration event of 5 adults, and a 20% fitness cost required ~ 2 additional daily migration events. Of note, a 2% increase in the adult mortality rate or a 2-day increase in the duration of the larval stage led to ~ 1 fewer daily migration event required for spread to Trinity Park, making this now very achievable—i.e., ~ 2–3 migration events for batches of 5 adults and ~ 1–2 migration events for batches of 10 adults.
Finally, we conducted an analysis of the sensitivity of our results to population structure, exploring the impact of (i) removing all population structure by treating Yorkeys Knob and Trinity Park as randomly mixing populations and (ii) incorporating heterogeneity in mosquito household population size. Results of the comparison to panmictic populations are depicted in Additional file 6: Fig. S5. Previously, we had seen that introducing population structure greatly increases the release requirements to eliminate UDMEL from a community (Additional file 4: Fig. S3). The trend of higher release requirements in structured populations is also seen for translocations, although to a lesser extent, with one additional release required for either replacement or remediation at a coverage of 75%, two additional releases required at a coverage of 50%, and 7–8 additional releases required at a coverage of 25%. Invasion of a neighboring population, on the other hand, requires ~ 1–2 fewer daily migration events in structured populations for batches of 10 adults having either the translocation or UDMEL.
Incorporating heterogeneity in household mosquito population size, we retain a mean of 15 adults per household, as inferred from Wolbachia field trial data in Yorkeys Knob [22], and distribute population sizes across households according to a zero-inflated, truncated exponential distribution with 55% of households having no mosquitoes and none having more than 45 adults. This distributional form, including zero inflation, is based on results of a large field survey conducted across a set of households in Kamphaeng Phet province, Thailand [49]. Including this source of heterogeneity substantially increases release requirements for replacement and remediation with translocations, with 3–4 additional releases required at a coverage of 75%, 5–6 additional releases required at a coverage of 50%, and 20 releases being insufficient at a coverage of 25% (Additional file 7: Fig. S6). Release requirements are marginally increased for UDMEL, with ~ 1–2 additional releases required at coverages of 25–100%. These results are likely due to the threshold frequency being more difficult to exceed in households with large numbers of mosquitoes, and this being less of an issue for UDMEL due to its lower threshold frequency. Fortunately, population size heterogeneity makes confinement more promising for both systems, increasing the required number of daily migration events for batches of 10 adults by 3–4 for translocations and by 1–2 for UDMEL.