Modelling Allee effects in a transgenic mosquito population during range expansion

ABSTRACT

Mosquitoes are vectors for many diseases that cause significant mortality and morbidity. As mosquito populations expand their range, they may undergo mate-finding Allee effects such that their ability to successfully reproduce becomes difficult at low population density. With new technology, creating target specific gene modification may be a viable method for mosquito population control. We develop a mathematical model to investigate the effects of releasing transgenic mosquitoes into newly established, low-density mosquito populations. Our model consists of two life stages (aquatic and adults), which are divided into three genetically distinct groups: heterogeneous and homogeneous transgenic that cause female infertility and a homogeneous wild type. We perform analytical and numerical analyses on the equilibria to determine the level of saturation needed to eliminate mosquitoes in a given area. This model demonstrates the potential for a gene drive system to reduce the spread of invading mosquito populations.

KEYWORDS:

• Mosquito dynamics
• Allee effect
• gene drive
• mathematical model

1. Introduction

Mosquitoes have been labelled the deadliest animal [19] as more than half a million people die each year from mosquito-borne diseases, including dengue and malaria [36,37]. Methods such as insecticide-treated mosquito nets, indoor residual spraying, space treatment, and larvicide application have led to significant headway in disease control by directly targeting the mosquito populations [18,34]. However, there are limitations to these methods. For example, extensive usage of common insecticides has led to mosquitoes rapidly acquiring resistance [20]. It is not a straightforward task to replace insecticides used in treated mosquito nets, which require low toxicity due to their close proximity to humans [20]. Effective control of mosquito populations may consequently require novel control measures to supplement existing tactics.

Beyond the rise of insecticide resistance, mosquito populations continue to expand into new territory and increase in density in other locations, further altering the potential for profound impacts on local transmission dynamics of vector-borne diseases. In the last 30 years, for example, Aedes albopictus has spread from Southeast Asia to Africa, Europe and the Americas due in part to an increase in global shipping [16]. In 2017, the Center for Disease Control showed that the spread of Ae. albopictus and Aedes aegypti increased in the United States [10]. Although Anopheles gambiae have not spread outside of Africa, their geographical range has shifted [31]. While mosquitoes are endemic to many locations, other places observe seasonal variation in their populations. As temperatures continue to increase globally, African highlands have observed increases in mosquito populations [23]. A study in Peru also revealed that weather is related to variation in Ae. aegypti populations and, importantly, has direct consequences on the transmission of dengue virus between humans and vectors [8]. Although eliminating mosquitoes that serve as disease vectors is insurmountable in many locations, minimizing their spatial spread to new areas may be a viable option.

One method that shows promise in facilitating control of mosquitoes and the infectious diseases they carry is the release of transgenic or otherwise genetically altered mosquitoes [3,17]. Manipulation of mosquito genetics has already been undertaken using the CRISPR-Cas9 system [14,17], such as on an Anopheles stephensi gene that exhibits anti-Plasmodium falciparum activity [14]. In An. gambiae, three recessive genes have been identified that produce female sterility. Using CRISPR-Cas9 constructs that function as a gene drive system, they observed transmission rates of these genes at over 90% [17]. Gene drive occurs when a gene is inherited with greater than random probability, in contrast to traditional Mendelian inheritance, where each gene has the same probability of being passed on to progeny. As the gene which creates infertility is more favourably inherited than the gene which does not, these results suggest that it is biologically feasible to use gene drive to express infertility in female Anopheles mosquitoes.

An additional feature of the population dynamics that may play a crucial role in the design and implementation of control via transgenic mosquitoes is the impact of low population density on the the ability of mosquitoes to successfully find a mate. Low population densities are most likely to occur as mosquitoes move into new territories and, consequently, this may generate a mate-finding Allee effect. More generally, an Allee effect is the loss of fitness or ability to reproduce due to decreased population size and arises from several mechanisms including mate-finding failure, foraging efficiency, and predation [6]. In particular, strong Allee effects exist if there is a population threshold below which per capita growth rates become negative [29]. Mosquitoes may have a naturally occurring Allee effect supported by evidence of skip oviposition in Ae. aegypti [39]. This preference implies that larvae survival may decrease with a low population density, suggesting the presence of an Allee effect.

We develop a mathematical model that describes the population dynamics of emerging mosquito populations following the introduction of transgenic mosquitoes and in the presence of a potential Allee effect. We conduct a stability analysis and perform numerical simulations of our model to determine conditions necessary for effective control. We find the ability of transgenic mosquitoes with fertility costs to reduce or eliminate the total population, with the magnitude of reduction dependent on the strength of the gene drive, the reduction in fertility, and the Allee effect.

2. Materials and methods

We develop two models of mosquito population dynamics. The first model describes the basic population dynamics of mosquitoes under the assumption that there is a pre-existing mate-finding Allee effect. The second model expands this framework by introducing a transgenic population of mosquitoes with genetic variation of two alleles at a single locus.

2.1. Basic model

We begin with a continuous-time model of the mosquito life cycle comprised of aquatic (J) and adult (A) classes. Egg, larvae, and pupae stages are included in the single aquatic class (J), in line with previous work [4,22].

We let r be the proportion of the adults that are males, and, thus, 1−r are females. Females have an associated fecundity of β eggs, and transition from the aquatic class to the adult class occurs at rate α. We assume there are constant, density-independent mortality rates ${\mu }_J$ and ${\mu }_A$ for the aquatic class and adult class, respectively. As there is evidence of density-dependent mortality in the first larval stage [13], we additionally assume that the aquatic class (J) has an associated density-dependent mortality rate ${\mu }_0$, which is in line with previous mosquito models [1,4,33]. Thus, in the absence of an Allee effect our basic model is given by (1) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =(1-r)\beta A-({\mu }_J+{\mu }_{0}J)J-\alpha J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\alpha J-{\mu }_{A}A.\end{array}$(1) We introduce a mate-finding Allee effect using the functional form $rA/(rA+\delta )$, where δ is the Allee constant, which dictates the strength of the Allee effect. A larger value indicates a stronger Allee effect. Thus, our basic model including an Allee effect is (2) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{r(1-r)\beta }{rA+\delta }{A}^2-({\mu }_J+{\mu }_{0}J)J-\alpha J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\alpha J-{\mu }_{A}A.\end{array}$(2) Notice that systems (1(1) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =(1-r)\beta A-({\mu }_J+{\mu }_{0}J)J-\alpha J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\alpha J-{\mu }_{A}A.\end{array}$(1) ) and (2(2) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{r(1-r)\beta }{rA+\delta }{A}^2-({\mu }_J+{\mu }_{0}J)J-\alpha J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\alpha J-{\mu }_{A}A.\end{array}$(2) ) are identical when $\delta =0$, that is, there is no Allee effect.

2.2. Gene drive

We extend our basic model to consider the effect of including and releasing transgenic mosquitoes into the population (Figure 1). Specifically, we consider transgenic mosquitoes that harbour a trait causing sterility that can be expressed by a single gene determined by two alleles [14]. We assume transmission of the gene does not follow Mendelian inheritance.

Figure 1. Schematic of the mosquito life-cycle model. Aquatic stages (J) and adult stages (A) are denoted by their genetic composition: wild type (W) or gene drive (D). All mating crosses are possible and the offspring are generated according to Table A1. Mortality occurs with at a density-dependent rate from the aquatic stages and a constant rate from the adults stages, both denoted by a diagonal arrow.

Let W denote the naturally occurring ‘wild type’ allele and D the new allele, which inhibits fertility and displays the gene drive phenomenon. Adults are categorised by their allele representation as denoted by subscripts. We assume that $A_{WW}$ individuals reproduce as before, but all homogeneous gene drive females ($A_{DD}$) have an associated fertility cost given by $1-f_d$ (where $f_d$ is defined as the fertility of the homogeneous females). The alleles are co-dominant so that $A_{DW}$ may have some fitness cost to fertility 1−f (where f is defined as the fertility of the heterogeneous females). We impose parameter bounds for fitness such that $0\le f_d,f\le 1$. The total adult population is now given by $A=A_{WW}+A_{DW}+A_{DD}$.

From the basic model, a proportion $rA/(rA+\delta )$ of females successfully mate where rA represents the total number of adult males in the population, which can be broken down by the type of male so that $\frac{rA}{rA+\delta }=\frac{r{A}_{WW}}{rA+\delta }+\frac{r{A}_{DW}}{rA+\delta }+\frac{r{A}_{DD}}{rA+\delta }.$ The fecundity for females of type $A_{WW}$ is β. Due to a fitness cost, heterogeneous females ($A_{DW}$) have a fecundity that is a fraction f of that for $A_{WW}$ females ($f\beta$). Under these assumptions, for example, the probability that a single female type DW mates with a type WW male is $\left(r{A}_{WW}/(rA+\delta )\right)$, and the female produces $f\beta (1-r)$ eggs.

Gene drive allows for alleles to be spread at a greater rate than Mendelian inheritance. We will assume that x is the strength of spread of the new allele where $0.5\le x\le 1$. When $x=0.5$, the gene exhibits normal Mendelian inheritance. The contribution of each adult class to the aquatic population is found by considering all possible matings as follows $\beta A_{female}\left(\frac{A_{male}}{A_{male}+\delta }\right)=\beta (1-r)\left(A_{WW}+f{A}_{DW}+f_d{A}_{DD}\right)\left(\frac{r\left(A_{WW}+A_{DW}+A_{DD}\right)}{rA+\delta }\right).$ We then split up the offspring of each male-female pairing into the appropriate aquatic group. For example, a heterogeneous gene drive ($A_{DW}$) male mating with a homogeneous wild type ($A_{WW}$) female, produces only two types of offspring: homogeneous wild ($J_{WW}$) or heterogeneous gene drive ($J_{DW}$). Based on the strength of the gene drive, a proportion of the offspring (1−x) are homogeneous wild type and the others (x) are heterogeneous gene drive. Details of each mating are provided in Table A1.

Our full model with gene drive is given by (3) $\begin{array}{rl}\frac{\mathrm{d}{J}_{WW}}{\mathrm{d}t}& =\frac{\beta (1-r)r}{rA+\delta }(A_{WW}^2+(1-x)(f+1)A_{DW}{A}_{WW}+(1-x{)}^{2}f{A}_{DW}^2)\\ & -({\mu }_J+{\mu }_0(J_{WW}+J_{DW}+J_{DD}))J_{WW}-\alpha J_{WW},\\ \frac{\mathrm{d}{J}_{DW}}{\mathrm{d}t}& =\frac{\beta (1-r)r}{rA+\delta }(x(1+f)A_{WW}{A}_{DW}+(1+f_d)A_{WW}{A}_{DD}+2x(1-x)f{A}_{DW}^2\\ & +(1-x)(f+f_d)A_{DW}{A}_{DD})\\ & -({\mu }_J+{\mu }_0(J_{WW}+J_{DW}+J_{DD}))J_{DW}-\alpha J_{DW},\\ \frac{\mathrm{d}{J}_{DD}}{\mathrm{d}t}& =\frac{\beta (1-r)r}{rA+\delta }(f{x}^2{A}_{DW}^2+x(f+f_d)A_{DD}{A}_{DW}+f_d{A}_{DD}^2)\\ & -({\mu }_J+{\mu }_0(J_{WW}+J_{DW}+J_{DD}))J_{DD}-\alpha J_{DD},\\ \frac{\mathrm{d}{A}_{WW}}{\mathrm{d}t}& =\alpha J_{WW}-{\mu }_A{A}_{WW},\\ \frac{\mathrm{d}{A}_{DW}}{\mathrm{d}t}& =\alpha J_{DW}-{\mu }_A{A}_{DW},\\ \frac{\mathrm{d}{A}_{DD}}{\mathrm{d}t}& =\alpha J_{DD}-{\mu }_A{A}_{DD}.\end{array}$(3)

3. Results and discussion

We begin with a general analysis of our models by determining the equilibria of the basic model of mosquito population dynamics and discussing their stability. We then find the equilibria of the model that includes transgenic mosquitoes in the context of fixed parameter sets. Finally, we present numerical simulations of the model under variable parameter sets.

3.1. Analysis of basic model of mosquito population dynamics

To ensure biological realism, we assume all parameter values are positive. In the absence of Allee effects ($\delta =0$), the basic model of mosquito population dynamics has two possible equilibria: the first corresponds to population extinction, $(J^{\ast },A^{\ast })=(0,0)$; and the second is a non-zero equilibrium, $(J^{\ast },A^{\ast })=\left(\right({\mu }_A/\alpha )A_b^{\ast },A_b^{\ast })$ where $A_b^{\ast }=\frac{(1-r){\alpha }^2\beta -{\mu }_J{\mu }_A\alpha -{\mu }_A{\alpha }^2}{{\mu }_A^2{\mu }_0}.$ The equilibria change stability at $r_1=\alpha (1-r)\beta /(\alpha {\mu }_A+{\mu }_J{\mu }_A)$. If $r_1<1$, then the extinction equilibrium is stable and $(J^{\ast },A^{\ast })=\left(\right({\mu }_A/\alpha )A_b^{\ast },A_b^{\ast })$ is unstable. If $r_1>1$, the opposite is true, that is, the extinction equilibrium is unstable and $(J^{\ast },A^{\ast })=\left(\right({\mu }_A/\alpha )A_b^{\ast },A_b^{\ast })$ is stable. Noticing the divergence of the vector field is always negative for biologically realistic parameter values and population sizes, applying Bendixson's criterion, there are no closed orbits. This implies when $r_1<1$ the extinction equilibrium is globally asymptotically stable and when $r_1>1$ the positive equilibrium is globally asymptotically stable.

When an Allee effect exists ($\delta >0$) in the basic model, three equilibria are possible. For all biologically relevant parameter values, there exists a locally stable extinction equilibrium $(J^{\ast },A^{\ast })=(0,0)$. In addition, there may be two other real-valued equilibria in the system, arising from a saddle node bifurcation at $r_2=\frac{({\mu }_0{\mu }_A^2\delta +{\alpha }^2{\mu }_{A}r+\alpha {\mu }_J{\mu }_{A}r-{\alpha }^{2}r(1-r)\beta {)}^2}{4r\alpha \delta {\mu }_0{\mu }_A^3({\mu }_J+\alpha )}=1.$ When $|r_2|<1$, the non-zero equilibria are imaginary and the extinction equilibrium is globally asymptotically stable as there are no closed orbits (Bendixson's criterion). When $|r_2|>1$, there exist two real non-zero equilibria, both with the same sign as the term squared in the numerator, that is, $s={\mu }_0{\mu }_A^2\delta +{\alpha }^2{\mu }_{A}r+\alpha {\mu }_J{\mu }_{A}r-{\alpha }^{2}r(1-r)\beta$. The two positive equilibria in the basic model when s<0 are given by $\begin{array}{rl}{A}^{\ast }& =\frac{-s\pm \sqrt{s^2-4r\alpha \delta {\mu }_0{\mu }_A^3({\mu }_J+\alpha )}}{2{\mu }_0{\mu }_A^{2}r},\\ J^{\ast }& =\frac{{\mu }_A}{\alpha }{A}^{\ast }.\end{array}$ A bifurcation diagram of how the the positive equilibria depend on the Allee constant is given in Figure 2(a). Here, the smaller of the two positive equilibria corresponds to the Allee threshold (the population size below which the per capita growth rate is negative) and the larger equilibrium corresponds to the carrying capacity (the population size above which the growth rate becomes negative). Therefore, the Allee threshold is unstable and the carrying capacity is locally stable (as shown in Appendix A.1).

Figure 2. Equilibrium values and stability change with Allee effect and density-dependent mortality. (a) A bifurcation diagram as the Allee constant δ varies. The equilibrium value at the bifurcation, marked by a star in (a), is highly correlated with the Allee constant (b) and density-dependent mortality (c). The diagram in (a) is computed with the parameters from Table 1. The LHS parameter sets were used to calculate the bifurcation values for (b) and (c). Parameter sets are the same as those used in Figure 3.

Figure 3. Effect of varying gene drive on the equilibrium of the full model. (a) Without an Allee effect, gene drive stronger than 95% is necessary for extinction with most parameter sets. (b) With an Allee effect, population extinction is attained for a wider range of gene drive strengths. The colour represents the equilibrium population size, with brighter colours associated with larger population size (visualized on a log-scale). Parameter sets are ordered on the x-axis by their equilibrium population size when $\delta =0$. Parameters are found in Table 1.

To illustrate the dependence of the saddle bifurcation point (red cross in Figure 2(a)) on parameterization, we select 192 parameter sets (details below) and determine the Allee constant δ such that $A_1^{\ast }=A_2^{\ast }$, that is, at the saddle node bifurcation. We graph the computed bifurcation values against the Allee constant δ and the density-dependent mortality rate of larva ${\mu }_0$ in Figure 2(b,c), respectively. We observe a linear relationship in log-log space with a positive correlation ($R^2=0.9874$) between the Allee constant and the bifurcation value in contrast to a negative correlation ($R^2=0.9305$) between the density-dependent death rate and the bifurcation value.

3.2. Analysis of dynamics with transgenic mosquitoes

In the full model including gene drive, it is difficult to find a comprehensible analytical solution. We first consider the case where the equilibrium for $J_{DW}=0$, which implies that $A_{DW}=0$. At equilibrium, the differential equation for $J_{DW}$ becomes $\frac{\mathrm{d}{J}_{DW}}{\mathrm{d}t}=\frac{\beta (1-r)r}{rA+\delta }\left((1+f_d)A_{WW}{A}_{DD}\right)=0.$ Notice that the equilibrium for $\mathrm{d}{J}_{DW}/\mathrm{d}t=0$ only occurs if either $A_{WW}=0$ or $A_{DD}=0$. First, if $A_{DD}=0$, as we are under the assumption of $A_{DW}=0$, we return to the same results as the basic model. Thus, we consider the case where $A_{WW}=0$, in addition to $A_{DW}=0$. With the assistance of Mathematica, we obtain the following equilibria for $A_{DD}$ $\begin{array}{rl}{A}_{DD}& =\frac{-1}{2r{\mu }_0{\mu }_A^2}\left(\delta {\mu }_0{\mu }_A^2-f_{d}r{\alpha }^2\beta (1-r)+r\alpha {\mu }_A(\alpha +{\mu }_J)\phantom{\sqrt{\left(f_d\right(1-r)r{\alpha }^2\beta +{\mu }_A(r\alpha (\alpha +{\mu }_J)+\delta {\mu }_0{\mu }_A){)}^2-4r\alpha \delta {\mu }_0(\alpha +{\mu }_J){\mu }_A^3}}\right.\\ & \pm \left.\sqrt{\left(f_d\right(1-r)r{\alpha }^2\beta +{\mu }_A(r\alpha (\alpha +{\mu }_J)+\delta {\mu }_0{\mu }_A){)}^2-4r\alpha \delta {\mu }_0(\alpha +{\mu }_J){\mu }_A^3}\right).\end{array}$ The equilibria are not dependent on the fertility of the heterozygotes f, but the stability of these equilibria do change based on the parameter f. We examine the stability of these equilibria where $0\le f_d,f\le 1$ and $\delta =0$ or 150. Biologically if $A_{DD}<1$, only a fraction of a mosquito exists, so we will only discuss results where there is at least one mosquito. Further discussion of the stability of equilibria in the gene drive model can be found in Appendix A.2.

3.2.1. Absence of an Allee effect

Consider the case when there is no Allee effect $\delta =0$ and there can exist one positive $A_{WW}$ equilibrium and one positive $A_{DD}$ equilibrium in addition to the extinction equilibrium. As we are working under the assumption $A_{DW}=J_{DW}=0$ there are no positive equilibria where both $A_{WW}$ and $A_{DD}$ co-exist. The stability of the $A_{WW}$ equilibrium depends on the gene drive x and the fertility of the heterozygotes f. In the case of the parameters from Table 1, the wild type only equilibrium is locally stable when $x<1/(1+f)$. Notice that the stability is not affected by the fertility of the homozygotes $f_d$. In contrast, the stability of the $A_{DD}$ equilibrium additionally depends on the fertility of the homozygotes $f_d$ as well as the gene drive x and the fertility of the heterozygotes f. Furthermore, with the chosen parameter set (Table 1) when $x<1/(1+f)$, there exists multi-stability between the wild type only, gene drive only, and extinction equilibria. In other words, the initial conditions determine when the adult population will be composed entirely of $A_{WW}$ or $A_{DD}$ mosquitoes (or no mosquitoes).

Table 1. Parameter ranges and a specific parameter set for An. gambiae at ${26}^{\circ }$C.

			An. gambiae
Symbol	Description	Range${}^{\mathrm{a}}$	Value	Reference
β	Eggs laid per female	$[50/4,122/2]$	100/3	[15,24]
r	Proportion males	$[0.5,0.55]$	0.50	[32]
${\mu }_0$	Density-dependent mortality rate	$[3\cdot {10}^{-6},0.5]$	0.05	–
${\mu }_J$	Mortality rate for aquatic stage	$[0.05,0.8]$	0.2	–
${\mu }_A$	Mortality rate for adults	$[1/30,1/8]$	1/23	[25]
α	Maturation rate	$[1/30,1/8]$	1/9.8	[5]
δ	Allee constant	$[0,250]$	0	–
1−f	Fertility cost of $A_{DW}$	[0, 1]	0.5	–
$1-f_d$	Fertility cost of $A_{DD}$	[0, 1]	0	–
x	Gene drive strength	[0.5, 1]	0.5	–

Notes: The origin of the parameter ranges and ${}^{\mathrm{a}}$references are found in the text. The density-dependent and independent mortality for the aquatic stage of An. gambiae varied widely in the literature, as noted in the text. Below the horizontal line are parameters for which Anopheles specific parameterization does not exist so standard values used are listed.

When $x>1/(1+f)$, the only possible equilibria either are composed entirely of gene drive homozygous mosquitoes or no mosquitoes at all, that is, extinction. There is a clear division of this space into where the gene drive equilibrium is stable or unstable, but the analytical formula is complicated and rather uninformative. Numerically, with the parameterization in Table 1, the $A_{DD}$ equilibrium is unstable when $x>\frac{1}{f+f_d}\left(-3.86\cdot {10}^{-3}+f+0.5f_d-0.5\sqrt{5.97\cdot {10}^{-5}-1.55\cdot {10}^{-2}{f}_d+f_d^2}\right).$ Note, these values were commuted in Mathematica and rounded.

3.2.2. Presence of an Allee effect

When an Allee effect is included, $\delta \ne 0$, more equilibria are possible, even with the simplifying assumption $A_{DW}=J_{DW}=0$. For example, consider $\delta =150$, with the parameterization in Table 1. There are five possible equilibria: the extinction equilibria, two positive equilibria with only wild type mosquitoes ($A_{WW}>0$) and two positive equilibria with only gene drive homozygotes ($A_{DD}>0$). The latter two equilibria are only real and positive when $f_d>\frac{r(1-r){\alpha }^2\beta {\mu }_A\left(r\alpha (\alpha +{\mu }_J)+\delta {\mu }_0{\mu }_A\right)+2\sqrt{r^3(1-r{)}^2{\alpha }^5{\beta }^2\delta {\mu }_0(\alpha +{\mu }_J){\mu }_A^3}}{r^2(1-r{)}^2{\alpha }^4{\beta }^2}$ or numerically with the parameterization from Table 1, when $f_d⪆0.243$. The larger equilibria in magnitude for $A_{WW}$ and $A_{DD}$ demonstrate qualitatively similar local stability patterns to those observed for the non-zero equilibrium in the basic model. The smaller $A_{WW}$ equilibrium is always unstable. The smaller $A_{DD}$ equilibrium has more complicated stability changes, which we only determine numerically.

3.2.3. Coexistence equilibria

It is possible for equilibria where $A_{WW},A_{DW},A_{DD}>0$ to exist in our system. However, these are analytically difficult to determine. For particular parameter sets we find coexistence equilibria numerically, for example in the lower right in Figure 6. In addition, there are multiple equilibria which are not biologically relevant, that is, $A_i<0$ with $i\in \{WW,DW,DD\}$.

3.3. Numerical results

3.3.1. Parameter ranges

Parameter values are highly species specific and dependent on many factors including temperature, diet, and mosquito size. We consider mosquitoes of genus Anopheles and Aedes, and specifically focus on An. gambiae, as these mosquitoes are primary vectors for the diseases that cause significant global public health impact [9,26,35–38]. We set upper and lower bounds for each parameter based on published experimental data, and vary the parameters within these bounds (Table 1). Based on temperature-specific data, we also fix a specific parameter set for An. gambiae at approximately 26${}^{\circ }$C as a baseline (Table 1).

Parameters are ascertained from the life history of mosquitoes. The fecundity β is determined by both the number of eggs per gonotrophic cyclic (approximately 50 to 122 eggs per gonotrophic cycle) and the number of days in a gonotrophic cycle (laying every 2 to 4 days in favourable environmental conditions) [7,11,12,15]. The proportion of males r is similar or slightly higher than females [7,12]. The transition from eggs to adults α takes between 8 and 30 days [5,12,27,32]. Adult mosquito life expectancy, related to adult mortality rate ${\mu }_A$, varies between 8 and 30 days [2,25,28]. Survival of larvae to adulthood, related to aquatic mortality rate ${\mu }_J$, is extremely temperature dependent and can vary from 3% to 93% [5,12,27,32]. Determining the density-dependent mortality rate from experimental evidence is more challenging. Previous estimates for this parameter vary significantly [13,21], so we allow ${\mu }_0$ to vary over a broad range of values.

The value of the Allee constant δ is not known for mosquitoes, although there exist estimates for other insect species. For example, in a study on gypsy months in different states, two different Allee threshold values were determined: 0.78% and 3.1% of the estimated carrying capacity [30]. As a result, we assume the mate-finding Allee effect to be a small percentage of the carrying capacity, if it exists. We choose δ such that the Allee effect is at most 15%. We run all simulations with $\delta =0$ and $\delta \in [50,250]$.

3.3.2. Parameter variation

We use Latin Hypercube Sampling (LHS) and filter our parameterization based on biologically derived criteria in the absence of gene drive. Specifically, we use n=225 LHS samples and a two-step filter process on the data. The first filter step chooses parameter sets such that there exist only non-zero equilibria that are real and positive. This amounts to only allowing parameter sets which fall to the left of the bifurcation value (star in Figure 2(a)). Approximately 85% of the total LHS samples are retained. Biologically if there exists a mate-finding Allee effect, we expect the number of mosquitoes necessary to find a mate to be relatively small compared to the carrying capacity. The relative Allee effect can be calculated as $A_1^{\ast }/A_2^{\ast }$. We therefore use a second step in our filter process to refine our parameter sets to those with relative Allee effects of less than $15\mathrm{\%}$. This retains ∼99% of the parameter sets from the first stage of the filter process. In total, we are left with 190 parameter sets that fulfill both criteria.

3.3.3. Simulation results

We carry forward the parameterization from the basic model for the gene drive model. We vary all the additional parameters found in the gene drive model: x, f, and $f_d$. We assume greater-than-Mendelian inheritance bias of the gene drive allele, that is, x>0.5. We consider the effect of the gene drive allele on the fertility of mosquitoes, such that $0\le f,f_d\le 1$.

To begin, we fixed the fertility of heterozygotes to be f=0.5 and homozygous gene drive to be $f_d=0$, and varied the gene drive strength for our parameter sets. We numerically determined the equilibria of the gene drive model with an initial population size set just above the Allee threshold $A_1^{\ast }$ found from the basic model. In Figure 3, the parameter sets are ordered by the non-zero equilibrium from the basic model without an Allee effect, although the relative Allee effect ($A_1^{\ast }/A_2^{\ast }$) gives a similar ordering (not pictured). From Figure 3, we find that gene drive must be stronger than 95% to drive the population to extinction in most parameter sets. At lower values of the equilibrium of the basic model increases, particularly when an Allee effect is present (further to the left in Figure 3(b)), a lower level of gene drive is sufficient to achieve elimination of the mosquito population. Notice that parameter sets which produce nearly identical equilibrium values, and are thus ordered sequentially along the x-axis, may require very different levels of gene drive to reach extinction (Figure 3). The reason for this variation in gene drive requirements is due to the values of individual parameters that differ between the parameter sets.

Next we relaxed our assumption on the fertility of heterozygotes, while the fertility of gene drive homozygotes remain fixed at zero ($f_d=0$), and varied gene drive strength. All other parameters were fixed to the Anopheles specific parameters found in Table 1. Elimination was observed for a wide range of f when gene drive was above 95% (Figure 4), regardless of the inclusion of an Allee effect. In the absence of an Allee effect, we find that based on the initial conditions, the adult population becomes either all wild type mosquitoes $A_{WW}$ or all homozygous gene drive type $A_{DD}$ (Figure 4(a)). In the presence of an Allee effect, we find that the population equilibrium decreases in all instances where gene drive mosquitoes previously dominated (Figure 4(b)). Notice, that there is a minimum level of fertility of heterozygotes for persistence of the gene drive population. Below this level mosquitoes with gene drive alleles lack sufficient progeny, so the allele cannot be passed on to the next generation (Figure 4, lower left corner).

Figure 4. Equilibrium for the full model when gene drive strength and heterozygous fertility are varied. (a) The equilibrium size when there is no Allee effect and (b) the fraction reduction in the population size when the Allee constant is $\delta =150$. Darker colours indicate smaller population size in (a) or smaller change in the presence of an Allee effect (b). Solid lines separate the space into regions with qualitatively different equilibria. Dashed lines show the regions on the other plot, for ease of comparison. Parameters from Table 1.

To consider variation in the fertility of both gene drive homozygotes and heterozygotes, we assume presence of an Allee effect and examine three scenarios of gene drive strength: 50%, 80% and 95%. Under Mendelian inheritance (x=0.5), the wild type only equilibrium is locally stable for all levels of homozygous and heterozygous fertility (Figure 5(a)). Depending on the initial conditions, however, it is also possible to reach a gene drive only equilibrium, as a region of multi-stability is present in the system. Furthermore, the region where only the wild type population exists is divided into two distinct regions: where no gene drive equilibrium is possible ($f_d\backslash lap0.243$) and where the equilibrium is unstable.

Figure 5. Bifurcation diagrams of the full model when varying both heterozygous and homozygous fertility but assuming the absence of heterozygotes. In the presence of an Allee effect ($\delta =150$), the regions of existence and stability change with increasing gene drive strength, x: (a) 50%, (b) 80%, and (c) 95%. Below the horizontal line $f\backslash lap0.243$ a non-zero equilibrium does not exist for the gene drive population. To the left of the vertical line in (b) and (c) the wild type only equilibrium is not stable. All other parameters from Table 1. Note: we assume $A_{DW}=J_{DW}=0$ for the determination of equilibria and stability.

Figure 6. Adult population size after three years when varying both heterozygous and homozygous fertility. The plots differ only in their initial levels of adult homozygous mosquitoes: (a) larger initial size of the wild type population relative to the gene drive population, $A_{WW}=500$, $A_{DD}=10$; (b) a larger initial size of the gene drive population $A_{WW}=10$, $A_{DD}=500$. Here a gene drive strength of 80% and an Allee effect with Allee constant $\delta =150$ are assumed. All other parameters from Table 1. The bifurcation lines in the absence of a heterozygous population (from Figure 5) are overlaid.

In the presence of gene drive (x>0.5) the region where the wild type only equilibrium is locally stable shrinks to low values of the heterozygous fertility, f (Figure 5(b,c)). For larger heterozygous fertility $f⪆30\mathrm{\%}$ and low homozygous fertility $f_d\backslash lap24\mathrm{\%}$ extinction is predicted assuming no heterozygous equilibria are possible, that is, $A_{DW}=J_{DW}=0$ (Figure 5(b,c), lower right corner). Our biological expectation is a monotonic increase of the fertility trait from heterozygotes to homozygotes. If we restrict ourselves to that region, below the y=x line in Figure 5, we see that it is a region that is dominated by extinction and gene drive only equilibria.

A detailed numerical examination of our model with an Allee effect supports our analytical results under our simplifying assumption of no heterozygotes. We discuss in detail when gene drive strength is 80%, but similar results occur for all gene drive considered (not pictured). In regions of known bi- or multi-stability our initial levels of gene drive or wild type mosquitoes matter (Figure 6(a,b), upper left corner). When the wild type population excludes the gene drive population, an identical wild type equilibrium is reached regardless of the fertility of the homozygous gene drive population. In contrast, when the gene drive population excludes the wild type population, the equilibrium population size is dependent on the fertility of both the heterozygotes and the homozygous gene drive population. When fertility of the gene drive population is very low ($f,f_d\backslash lap0.2$) then the population may get pushed below the Allee threshold causing it to crash (Figure 6(b)). Finally, in our numerical examination we also observe the persistence of population where we predicted the extinction of the population. This is due to the presence of co-existence equilibrium (Figure 6, positive values in the lower right). These co-existence equilibria require the persistence of the heterozygous population, that is, $A_{DW},J_{DW}>0$.

4. Conclusion

In this paper, we examined the impact of a novel control method on an invading mosquito population. Releasing sterile mosquitoes has been modelled previously; however one of the biggest issues with sterile release is the success of the released mosquitoes in competing with the native populations. For this to be successful, most models show a need for continual releases of sterile mosquitoes. A more efficient method would be if the gene of sterility could spread through the population. To investigate this, we developed a mathematical model of mosquito populations and then introduced gene drive. Our results show that it can be effective, but the drive through the population needs to be at a high level in order to guarantee conditions which will stop invasion. Overall, if gene drive is over 99% then for all parameter values considered it will prove effective. This high level of gene drive is unlikely in reality (although rates of greater than 90% gene drive have been observed experimentally) as resistance of the gene in the host may occur [17]. However, our model shows that even with Mendelian inheritance (gene drive of 50%), it is possible to eliminate the invading mosquito population under certain conditions. In particular, if a strong Allee effect exists in mosquitoes, the use of gene drive is much more effective.

Preventing the invasion of mosquitoes into new territories will impact their ability to spread disease to those areas. In most infectious disease modelling, spread of the disease is quantified by the reproduction number $R_0$. When calculating $R_0$ for mosquito-borne diseases, the population of mosquitoes is an important component. A strong Allee effect decreases the population number, which would decrease $R_0$. Coupled with the effects of introducing gene drive into an invading mosquito population, our results suggest gene drive can offer a viable method of control in the efforts to reduce mosquito-borne diseases.

Disclosure statement

The authors declare no conflicts of interest.

ORCID

Lauren M. Childs  http://orcid.org/0000-0003-3904-3895

Additional information

Funding

LMC acknowledges the support of the Simons Foundation: Collaboration Grant for Mathematicians Award 524390.

Appendix 1. Gene drive model derivation

We provide a detailed description of all fecundity terms, that is, the result of all possible crosses in Table A1. These terms are included in System 3(3) $\begin{array}{rl}\frac{\mathrm{d}{J}_{WW}}{\mathrm{d}t}& =\frac{\beta (1-r)r}{rA+\delta }(A_{WW}^2+(1-x)(f+1)A_{DW}{A}_{WW}+(1-x{)}^{2}f{A}_{DW}^2)\\ & -({\mu }_J+{\mu }_0(J_{WW}+J_{DW}+J_{DD}))J_{WW}-\alpha J_{WW},\\ \frac{\mathrm{d}{J}_{DW}}{\mathrm{d}t}& =\frac{\beta (1-r)r}{rA+\delta }(x(1+f)A_{WW}{A}_{DW}+(1+f_d)A_{WW}{A}_{DD}+2x(1-x)f{A}_{DW}^2\\ & +(1-x)(f+f_d)A_{DW}{A}_{DD})\\ & -({\mu }_J+{\mu }_0(J_{WW}+J_{DW}+J_{DD}))J_{DW}-\alpha J_{DW},\\ \frac{\mathrm{d}{J}_{DD}}{\mathrm{d}t}& =\frac{\beta (1-r)r}{rA+\delta }(f{x}^2{A}_{DW}^2+x(f+f_d)A_{DD}{A}_{DW}+f_d{A}_{DD}^2)\\ & -({\mu }_J+{\mu }_0(J_{WW}+J_{DW}+J_{DD}))J_{DD}-\alpha J_{DD},\\ \frac{\mathrm{d}{A}_{WW}}{\mathrm{d}t}& =\alpha J_{WW}-{\mu }_A{A}_{WW},\\ \frac{\mathrm{d}{A}_{DW}}{\mathrm{d}t}& =\alpha J_{DW}-{\mu }_A{A}_{DW},\\ \frac{\mathrm{d}{A}_{DD}}{\mathrm{d}t}& =\alpha J_{DD}-{\mu }_A{A}_{DD}.\end{array}$(3) .

Table A1. Origin of aquatic populations from crosses of wild type and gene drive adults.

		Fraction progeny
Female	Male	WW	DW	DD	Coming from adults
WW	WW	1	0	0	$\beta (1-r)A_{WW}\left({\displaystyle \frac{r{A}_{WW}}{rA+\delta }}\right)$
WW	DW	$(1-x)$	x	0	$\beta (1-r)A_{WW}\left({\displaystyle \frac{r{A}_{DW}}{rA+\delta }}\right)$
WW	DD	0	1	0	$\beta (1-r)A_{WW}\left({\displaystyle \frac{r{A}_{DD}}{rA+\delta }}\right)$
DW	WW	$(1-x)$	x	0	$f\beta (1-r)A_{DW}\left({\displaystyle \frac{r{A}_{WW}}{rA+\delta }}\right)$
DW	DW	$(1-x{)}^2$	$2x(1-x)$	$x^2$	$f\beta (1-r)A_{DW}\left({\displaystyle \frac{r{A}_{DW}}{rA+\delta }}\right)$
DW	DD	0	$(1-x)$	x	$f\beta (1-r)A_{DW}\left({\displaystyle \frac{r{A}_{DD}}{rA+\delta }}\right)$
DD	WW	0	1	0	$f_d\beta (1-r)A_{DD}\left({\displaystyle \frac{r{A}_{WW}}{rA+\delta }}\right)$
DD	DW	0	$(1-x)$	x	$f_d\beta (1-r)A_{DD}\left({\displaystyle \frac{r{A}_{DW}}{rA+\delta }}\right)$
DD	DD	0	0	1	$f_d\beta (1-r)A_{DD}\left({\displaystyle \frac{r{A}_{DD}}{rA+\delta }}\right)$

Appendix 2. Derivation of equilibria

A.1. Basic model

In the main text we find the equilibria $(J^{\ast },A^{\ast })$ of the basic system. Here we discuss the stability of these equilibria.

We determine the stability of the equilibria in the presence of an Allee effect, when $\delta >0$, by examining the eigenvalues of the Jacobian matrix (A1) $\mathcal{J}=\left|\begin{array}{cc}-\alpha -{\mu }_J-\frac{2{\mu }_0{\mu }_A}{\alpha }{A}^{\ast }& \frac{r(1-r)\beta A^{\ast }\left(r{A}^{\ast }+2\delta \right)}{(r{A}^{\ast }+\delta {)}^2}\\ \alpha & -{\mu }_A\end{array}\right|$(A1) evaluated at each equilibrium.

Considering the trivial equilibrium $(J^{\ast },A^{\ast })=(0,0)$ and evaluating the Jacobian, we find that our eigenvalues are ${\lambda }_1=-\alpha -{\mu }_J$ and ${\lambda }_2=-{\mu }_A$. As all parameters are non-negative, the eigenvalues are always negative, indicating the extinction equilibrium is always locally stable.

We now show that the non-zero equilibrium smaller in magnitude (i.e. closer to zero), $A_2^{\ast }$, is unstable and the one larger in magnitude, $A_1^{\ast }$, is stable. It is sufficient to test if the trace($\mathcal{J})<0$ and determinant($\mathcal{J})>0$ for local stability. Considering (A1(A1) $\mathcal{J}=\left|\begin{array}{cc}-\alpha -{\mu }_J-\frac{2{\mu }_0{\mu }_A}{\alpha }{A}^{\ast }& \frac{r(1-r)\beta A^{\ast }\left(r{A}^{\ast }+2\delta \right)}{(r{A}^{\ast }+\delta {)}^2}\\ \alpha & -{\mu }_A\end{array}\right|$(A1) ), the trace is always negative for biologically realistic parameters and equilibria, so if the determinant is positive, $\alpha {\mu }_A+{\mu }_J{\mu }_A+\frac{2{\mu }_0{\mu }_A^2}{\alpha }{A}^{\ast }-\frac{\alpha r(1-r)\beta A^{\ast }\left(r{A}^{\ast }+2\delta \right)}{(r{A}^{\ast }+\delta {)}^2}>0,$ the equilibrium will be stable. Let us assume that we have positive equilibriums, guaranteed when $s:={\mu }_0{\mu }_A^2\delta +{\alpha }^2{\mu }_{A}r+\alpha {\mu }_J{\mu }_{A}r-{\alpha }^{2}r(1-r)\beta <0$ as noted in the main text.

Then, we can rearrange determinant($\mathcal{J}$) as follows $\begin{array}{rl}& \alpha {\mu }_A+{\mu }_J{\mu }_A+\frac{2{\mu }_0{\mu }_A^2}{\alpha }{A}^{\ast }-\frac{\alpha r(1-r)\beta A^{\ast }\left(r{A}^{\ast }+2\delta \right)}{(r{A}^{\ast }+\delta {)}^2}\\ & =\alpha {\mu }_A+{\mu }_J{\mu }_A+\frac{2{\mu }_0{\mu }_A^2}{\alpha }{A}^{\ast }-\frac{\alpha r(1-r)\beta A^{\ast }}{r{A}^{\ast }+\delta }\frac{r{A}^{\ast }+2\delta }{r{A}^{\ast }+\delta }\\ \text{From equivalences of our non-zero equilibrium}& & =\alpha {\mu }_A+{\mu }_J{\mu }_A+\frac{2{\mu }_0{\mu }_A^2}{\alpha }{A}^{\ast }-\left({\mu }_J{\mu }_A+\frac{{\mu }_0{\mu }_A^2}{\alpha }{A}^{\ast }+{\mu }_A\alpha \right)\frac{r{A}^{\ast }+2\delta }{r{A}^{\ast }+\delta }& =\frac{{\mu }_0{\mu }_A^2}{\alpha }{A}^{\ast }+\left(\alpha {\mu }_A+{\mu }_J{\mu }_A+\frac{{\mu }_0{\mu }_A^2}{\alpha }{A}^{\ast }\right)\left(1-\frac{r{A}^{\ast }+2\delta }{r{A}^{\ast }+\delta }\right)\\ & =\frac{{\mu }_0{\mu }_A^2}{\alpha }{A}^{\ast }+\left(\alpha {\mu }_A+{\mu }_J{\mu }_A+\frac{{\mu }_0{\mu }_A^2}{\alpha }{A}^{\ast }\right)\left(\frac{-\delta }{r{A}^{\ast }+\delta }\right)\\ & & =\left(1-\frac{\delta }{r{A}^{\ast }+\delta }\right)\frac{{\mu }_0{\mu }_A^2}{\alpha }{A}^{\ast }-\frac{\delta \left(\alpha {\mu }_A+{\mu }_J{\mu }_A\right)}{r{A}^{\ast }+\delta }& =\frac{r{\mu }_0{\mu }_A^2(A^{\ast }{)}^2}{\alpha (r{A}^{\ast }+\delta )}-\frac{\delta {\mu }_A(\alpha +{\mu }_J)}{r{A}^{\ast }+\delta }\\ & =\frac{r{\mu }_0{\mu }_A^2(A^{\ast }{)}^2-\alpha \delta {\mu }_A(\alpha +{\mu }_J)}{\alpha (r{A}^{\ast }+\delta )}\\ & =\frac{r{\mu }_0{\mu }_A^2}{\alpha (r{A}^{\ast }+\delta )}\left((A^{\ast }{)}^2-\frac{\alpha \delta (\alpha +{\mu }_J)}{r{\mu }_0{\mu }_A}\right)\\ & =\frac{r{\mu }_0{\mu }_A^2}{\alpha (r{A}^{\ast }+\delta )}\left((A^{\ast }{)}^2-\frac{4r{\mu }_0{\mu }_A^3\alpha \delta (\alpha +{\mu }_J)}{(2r{\mu }_0{\mu }_A^2{)}^2}\right).\end{array}$ So we see that if (A2) $(A^{\ast }{)}^2-\frac{4r{\mu }_0{\mu }_A^3\alpha \delta (\alpha +{\mu }_J)}{(2r{\mu }_0{\mu }_A^2{)}^2}>0,$(A2) then the determinant($\mathcal{J}$) is positive and we have a stable positive equilibrium. In order to evaluate the when this occurs, let (A3) $A^{\ast }=\frac{-Q-W+Z\pm \sqrt{(Q+W-Z{)}^2-4QW}}{2{\mu }_0{\mu }_A^{2}r},$(A3) where $\begin{array}{rl}Q& ={\mu }_0{\mu }_A^2\delta ,\\ W& =\alpha {\mu }_{A}r(\alpha +{\mu }_J),\\ Z& ={\alpha }^{2}r(1-r)\beta .\end{array}$ Assuming there exists a positive equilibrium, then we require $\begin{array}{rl}(\ast )& (Q+W-Z{)}^2-4QW\ge 0,\\ (\ast \ast )& Z>Q+W.\end{array}$ Now if we plug $A^{\ast }$ into the condition found from (A2(A2) $(A^{\ast }{)}^2-\frac{4r{\mu }_0{\mu }_A^3\alpha \delta (\alpha +{\mu }_J)}{(2r{\mu }_0{\mu }_A^2{)}^2}>0,$(A2) ), we have the following (A4) $\begin{array}{rl}& {\left(\frac{-Q-W+Z\pm \sqrt{(Q+W-Z{)}^2-4QW}}{2{\mu }_0{\mu }_A^{2}r}\right)}^2-\frac{4QW}{(2r{\mu }_0{\mu }_A^2{)}^2}\\ & =\frac{{\left(-Q-W+Z\pm \sqrt{(Q+W-Z{)}^2-4QW}\right)}^2-4QW}{(2{\mu }_0{\mu }_A^{2}r{)}^2}\\ & =\frac{(-Q-W+Z{)}^2-4QW+(Q+W-Z{)}^2-4QW}{(2{\mu }_0{\mu }_A^{2}r{)}^2}\end{array}$(A4) (A5) $+\frac{\pm 2(-Q-W+Z)\sqrt{(Q+W-Z{)}^2-4QW}}{(2{\mu }_0{\mu }_A^{2}r{)}^2}$(A5) (A6) $=\frac{2\left((-Q-W+Z{)}^2-4QW\pm (-Q-W+Z)\sqrt{(Q+W-Z{)}^2-4QW}\right)}{(2{\mu }_0{\mu }_A^{2}r{)}^2}.$(A6) From (*), we know that the term on line (A4(A4) $\begin{array}{rl}& {\left(\frac{-Q-W+Z\pm \sqrt{(Q+W-Z{)}^2-4QW}}{2{\mu }_0{\mu }_A^{2}r}\right)}^2-\frac{4QW}{(2r{\mu }_0{\mu }_A^2{)}^2}\\ & =\frac{{\left(-Q-W+Z\pm \sqrt{(Q+W-Z{)}^2-4QW}\right)}^2-4QW}{(2{\mu }_0{\mu }_A^{2}r{)}^2}\\ & =\frac{(-Q-W+Z{)}^2-4QW+(Q+W-Z{)}^2-4QW}{(2{\mu }_0{\mu }_A^{2}r{)}^2}\end{array}$(A4) ) is positive. Considering only the larger equilibrium value $A_1^{\ast }$, we can also say that that the term on line (A5(A5) $+\frac{\pm 2(-Q-W+Z)\sqrt{(Q+W-Z{)}^2-4QW}}{(2{\mu }_0{\mu }_A^{2}r{)}^2}$(A5) ) is positive. This implies (A2(A2) $(A^{\ast }{)}^2-\frac{4r{\mu }_0{\mu }_A^3\alpha \delta (\alpha +{\mu }_J)}{(2r{\mu }_0{\mu }_A^2{)}^2}>0,$(A2) ) will be positive. This means the determinant($\mathcal{J}$) is positive and $A_1^{\ast }$ is stable.

We know want to determine the stability of $A_2^{\ast }$. Since all parameters will be positive, then Q,W,Z>0. Thus, $\begin{array}{rl}-4QW& <0\\ (-Q-W+Z{)}^2-4QW& <(-Q-W+Z{)}^2\\ \sqrt{(-Q-W+Z{)}^2-4QW}& <-Q-W+Z\\ \sqrt{(-Q-W+Z{)}^2-4QW}\sqrt{(-Q-W+Z{)}^2-4QW}& <(-Q-W+Z)\\ & \sqrt{(-Q-W+Z{)}^2-4QW}\\ (-Q-W+Z{)}^2-4QW& <(-Q-W+Z)\\ & \sqrt{(-Q-W+Z{)}^2-4QW}\end{array}$ (A7) $(-Q-W+Z{)}^2-4QW-(-Q-W+Z)\sqrt{(-Q-W+Z{)}^2-4QW}<0.$(A7) Now if we consider only the smaller equilibrium value $A_2^{\ast }$, we see that (A6(A6) $=\frac{2\left((-Q-W+Z{)}^2-4QW\pm (-Q-W+Z)\sqrt{(Q+W-Z{)}^2-4QW}\right)}{(2{\mu }_0{\mu }_A^{2}r{)}^2}.$(A6) ) is negative due to (A7(A7) $(-Q-W+Z{)}^2-4QW-(-Q-W+Z)\sqrt{(-Q-W+Z{)}^2-4QW}<0.$(A7) ). Thus, given a positive $A_2^{\ast }$, (A2(A2) $(A^{\ast }{)}^2-\frac{4r{\mu }_0{\mu }_A^3\alpha \delta (\alpha +{\mu }_J)}{(2r{\mu }_0{\mu }_A^2{)}^2}>0,$(A2) ) will be negative, demonstrating determinant($\mathcal{J}$) is negative and $A_2^{\ast }$ is unstable.

A.2. Gene drive model

The equilibria for the gene drive model are more complicated. The extinction equilibrium ($A_{WW}=A_{DW}=A_{DD}=J_{WW}=J_{DW}=J_{DD}=0$) always exists. Here we consider four other equilibria that occur when $A_{DW}=J_{DW}=0$. When $A_{DD}=J_{DD}=0$, we find the same two equilibria as in the basic model, such that only wild type mosquitoes persist, $A_{WW},J_{WW}>0$. Maintaining the assumption $A_{DW}=J_{DW}=0$, we also find two equilibria with only gene drive mosquitoes, $A_{DD},J_{DD}>0$. The stability of these equilibria differ from in basic model.

A.2.1. Stability of the extinction equilibrium

To begin, we examine the extinction equilibrium and find that there are only two eigenvalues, $-\alpha -{\mu }_J$ and $-{\mu }_A$, each with multiplicity of three. Thus, under our assumption of positive parameter values, the extinction equilibria is always stable.

A.2.2. Stability of the wild type only equilibrium

Consider $A_{WW}^{\ast }>0$ and $J_{WW}^{\ast }>0$ with all other populations zero. We find with the assistance of Mathematica the following eigenvalues of the Jacobian, written in terms of $A_{WW}^{\ast }$, as $\begin{array}{rl}{\lambda }_1& =-{\mu }_A,\\ {\lambda }_2& =-\left(\frac{{\alpha }^2+\alpha {\mu }_J+{\mu }_0{\mu }_A{A}_{WW}^{\ast }}{\alpha }\right),\\ {\lambda }_{3,4}& =-\left(\frac{Q_1+r{\mu }_0{\mu }_A(A_{WW}^{\ast }{)}^2+(\delta {\mu }_0{\mu }_A+r\alpha (\alpha +{\mu }_J+{\mu }_A))A_{WW}^{\ast }\pm \sqrt{Q_3}}{Q_2}\right),\\ {\lambda }_{5,6}& =-\left(\frac{Q_1+2r{\mu }_0{\mu }_A(A_{WW}^{\ast }{)}^2+(2\delta {\mu }_0{\mu }_A+r\alpha (\alpha +{\mu }_J+{\mu }_A))A_{WW}^{\ast }\pm \sqrt{Q_4}}{Q_2}\right),\end{array}$ where $\begin{array}{rl}{Q}_1& ={\alpha }^2\delta +\alpha \delta {\mu }_J+\alpha \delta {\mu }_A,\\ Q_2& =2\alpha (r{A}_{WW}^{\ast }+\delta ),\\ Q_3& =(r{A}_{WW}^{\ast }+\delta )\left({\alpha }^2\delta \right(\alpha +{\mu }_J-{\mu }_A{)}^2+r{\mu }_0^2{\mu }_A^2(A_{WW}^{\ast }{)}^3+{\mu }_0{\mu }_A(2r\alpha (\alpha +{\mu }_J-{\mu }_A)+\delta {\mu }_0{\mu }_A)\\ & (A_{WW}^{\ast }{)}^2+\alpha (-4(1+f)r^{2}x{\alpha }^2\beta +r\alpha ({\alpha }^2+({\mu }_J-{\mu }_A{)}^2+2\alpha \left(2\right(1+f)x\beta +{\mu }_J-{\mu }_A))\\ & +2\delta {\mu }_0(\alpha +{\mu }_J-{\mu }_A){\mu }_A\left)A_{WW}^{\ast }\right),\\ Q_4& ={\alpha }^2{\delta }^2(\alpha +{\mu }_J-{\mu }_A{)}^2+4r^2{\mu }_0^2{\mu }_A^2(A_{WW}^{\ast }{)}^4+4r{\mu }_0{\mu }_A\left(r\alpha \right(\alpha +{\mu }_J-{\mu }_A)+2\delta {\mu }_0{\mu }_A)(A_{WW}^{\ast }{)}^3\\ & +2A_{WW}^{\ast }\alpha \delta (-4r^2{\alpha }^2\beta +r\alpha ({\alpha }^2+({\mu }_J-{\mu }_A{)}^2+2\alpha (2\beta +{\mu }_J-{\mu }_A\left)\right)\\ & +2\delta {\mu }_0(\alpha +{\mu }_J-{\mu }_A){\mu }_A)+(-4r^3{\alpha }^3\beta +r^2{\alpha }^2({\alpha }^2+({\mu }_J-{\mu }_A{)}^2\\ & +2\alpha (2\beta +{\mu }_J-{\mu }_A)\left)\right)(A_{WW}^{\ast }{)}^2+(8r\alpha \delta {\mu }_0(\alpha +{\mu }_J-{\mu }_A){\mu }_A+4{\delta }^2{\mu }_0^2{\mu }_A^2\left)\right(A_{WW}^{\ast }{)}^2.\end{array}$ Given positive parameters, as well as ${\mu }_J>{\mu }_A$ and $A_{WW}>0$, then four of the eigenvalues are always negative, ${\lambda }_{1,2,3,5}<0$. It is possible for ${\lambda }_4$ and ${\lambda }_6$ to be positive. When either of these cases occur, the equilibrium would be unstable. Both ${\lambda }_4$ and ${\lambda }_6$ are dependent on the Allee constant δ. In addition, ${\lambda }_4$ is also dependent on the level of gene drive x and the heterozygous fertility f. Based on parameter values (Table 1), we can determine numerically where stability changes.

A.2.3. Stability of the gene drive only equilibrium

Consider $A_{DD}>0$ and $J_{DD}>0$ with all other populations zero. We find with the assistance of Mathematica the following eigenvalues of the Jacobian, written in terms of $A_{DD}^{\ast }$, as $\begin{array}{rl}{\lambda }_1& =-{\mu }_A,\\ {\lambda }_2& =-\left(\frac{{\alpha }^2+\alpha {\mu }_J+{\mu }_0{\mu }_A{A}_{DD}^{\ast }}{\alpha }\right),\\ {\lambda }_{3,4}& =-\left(\frac{Q_1+r{\mu }_0{\mu }_A(A_{DD}^{\ast }{)}^2+(\delta {\mu }_0{\mu }_A+r\alpha (\alpha +{\mu }_J+{\mu }_A))A_{DD}^{\ast }\pm \sqrt{Q_6}}{Q_5}\right),\\ {\lambda }_{5,6}& =-\left(\frac{Q_1+4r{\mu }_0{\mu }_A(A_{DD}^{\ast }{)}^2+2(\delta {\mu }_0{\mu }_A+r\alpha (\alpha +{\mu }_J+{\mu }_A))A_{DD}^{\ast }\pm \sqrt{Q_7}}{Q_5}\right),\end{array}$ where $\begin{array}{rl}{Q}_5& =2\alpha (r{A}_{DD}^{\ast }+\delta )\\ Q_6& =(r{A}_{DD}^{\ast }+\delta )\left({\alpha }^2\delta \right(\alpha +{\mu }_J-{\mu }_A{)}^2+(A_{DD}^{\ast }{)}^{3}r{\mu }_0^2{\mu }_A^2+{\mu }_0{\mu }_A(2r\alpha (\alpha +{\mu }_J-{\mu }_A)\\ & +\delta {\mu }_0{\mu }_A\left)\right(A_{DD}^{\ast }{)}^2+\alpha \left(4\right(f+f_d\left)r^2\right(-1+x){\alpha }^2\beta +2\delta {\mu }_0(\alpha +{\mu }_J-{\mu }_A){\mu }_A\\ & +r\alpha ({\alpha }^2+({\mu }_J-{\mu }_A{)}^2)A_{DD}^{\ast }-2{\alpha }^2(2f(-1+x)\beta +2f_d(-1+x)\beta -{\mu }_J+{\mu }_A\left)A_{DD}^{\ast }\right),\\ Q_7& ={\alpha }^2{\delta }^2(\alpha +{\mu }_J-{\mu }_A{)}^2+16r^2{\mu }_0^2{\mu }_A^2(A_{DD}^{\ast }{)}^4+16r{\mu }_0{\mu }_A\left(r\alpha \right(\alpha +{\mu }_J-{\mu }_A)+\delta {\mu }_0{\mu }_A)(A_{DD}^{\ast }{)}^3\\ & +4\alpha \delta (-2f_d{r}^2{\alpha }^2\beta +r\alpha ({\alpha }^2+({\mu }_J-{\mu }_A{)}^2+2\alpha (f_d\beta +{\mu }_J-{\mu }_A\left)\right)\\ & +\delta {\mu }_0(\alpha +{\mu }_J-{\mu }_A){\mu }_A)A_{DD}^{\ast }+4(-2f_d{r}^3{\alpha }^3\beta +r^2{\alpha }^2({\alpha }^2+({\mu }_J-{\mu }_A{)}^2\\ & +2\alpha (f_d\beta +{\mu }_J-{\mu }_A))+4r\alpha \delta {\mu }_0(\alpha +{\mu }_J-{\mu }_A\left){\mu }_A\right)(A_{DD}^{\ast }{)}^2+4{\delta }^2{\mu }_0^2{\mu }_A^2(A_{DD}^{\ast }{)}^2.\end{array}$ Similar to Section A.2.2, we find that given positive parameters, as well as ${\mu }_J>{\mu }_A$, and $A_{DD}^{\ast }>0$, then all of the eigenvalues are always negative, except ${\lambda }_4$ and ${\lambda }_6$. Both ${\lambda }_4$ and ${\lambda }_6$ are dependent on the gene drive strength x, Allee constant δ, and the fertility of the homozygous gene drive females $f_d$. Additionally, ${\lambda }_4$ is dependent on the heterozygous fertility f. Based on parameter values (Table 1), we can determine numerically where stability changes.

Appendix 3. Gene drive model with migration

In our models presented in the main text, we ignored the continued invasion of mosquitoes. The assumption of a single event populating a new population is biologically possible if a batch of eggs are laid in a new area and hatch at approximately the same time. To consider the continue invasion of mosquitoes, we adapted our basic model to include a time-decaying adult migration term: $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{r(1-r)\beta }{rA+\delta }{A}^2-({\mu }_J+{\mu }_{0}J)J-\alpha J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\alpha J-{\mu }_{A}A+M,\\ M& =m_0{\mathrm{e}}^{-ct},\end{array}$ where M is the migration rate of adult mosquitoes; c is the deceleration of migrating mosquitoes into the area, that is, the rate mosquitoes lose interest in migrating; and $m_0$ is the initial size of the migrating mosquito population. As under these assumptions the migrating population will eventually disappear, the addition of migration in this form does not effect the asymptotic dynamics. However, the decaying migration does effect the transient dynamics of the mosquito population.

We also make a similar adaptation to the gene drive model, that is, the addition of an extra equation for the migrating adult population, resulting in the system $\begin{array}{rl}\frac{\mathrm{d}{J}_{WW}}{\mathrm{d}t}& =\frac{\beta (1-r)r}{rA+\delta }\left(A_{WW}^2+(1-x)(f+1)A_{DW}{A}_{WW}+(1-x{)}^{2}f{A}_{DW}^2\right)\\ & -({\mu }_J+{\mu }_{0}J)J_{WW}-\alpha J_{WW},\\ \frac{\mathrm{d}{J}_{DW}}{\mathrm{d}t}& =\frac{\beta (1-r)r}{rA+\delta }\left[x\right(1+f)A_{WW}{A}_{DW}+(1+f_d)A_{WW}{A}_{DD}+2x(1-x)f{A}_{DW}^2\\ & +(1-x)(f+f_d)A_{DW}{A}_{DD}]-({\mu }_J+{\mu }_{0}J)J_{DW}-\alpha J_{DW},\\ \frac{\mathrm{d}{J}_{DD}}{\mathrm{d}t}& =\frac{\beta (1-r)r}{rA+\delta }\left(f{x}^2{A}_{DW}^2+x(f+f_d)A_{DD}{A}_{DW}+f_d{A}_{DD}^2\right)\\ & -({\mu }_J+{\mu }_{0}J)J_{DD}-\alpha J_{DD},\\ \frac{\mathrm{d}{A}_{WW}}{\mathrm{d}t}& =\alpha J_{WW}-{\mu }_A{A}_{WW}+M,\\ \frac{\mathrm{d}{A}_{DW}}{\mathrm{d}t}& =\alpha J_{DW}-{\mu }_A{A}_{DW},\\ \frac{\mathrm{d}{A}_{DD}}{\mathrm{d}t}& =\alpha J_{DD}-{\mu }_A{A}_{DD},\\ M& =m_0{\mathrm{e}}^{-ct},\end{array}$ where $J=J_{WW}+J_{DW}+J_{DD}$, $A=A_{WW}+A_{DW}+A_{DD}$, and all other parameters are defined as in the gene drive model. All migrating mosquitoes are assumed to be wild type.

All equilibria remain unchanged, but in the regions of multiple stability, which of the stable equilibria is reached may change due to the altered initial conditions. In other words, the migrating population contributes to the wild type population and acts like a larger initial influx of wild type mosquitoes. In addition, the inclusion of migrating mosquitoes can produce large changes in the transient behaviour (Figure A1).

Figure A1. Transient dynamics of the model with and without migration. (a) The inclusion of immigrating mosquitoes ($m_0=10$, c=0.5) boosts the transient mosquito population. (b) The rate of decay of migration c affects the size of the transient peak. Parameters from Table 1 with f=0.5, $f_d=0$, and x=0.9.