Interplay of population genetics and dynamics in the genetic control of mosquitoes

Some proposed genetics-based vector control methods aim to suppress or eliminate a mosquito population in a similar manner to the sterile insect technique. One approach under development in Anopheles mosquitoes uses homing endonuclease genes (HEGs)—selfish genetic elements (inherited at greater than Mendelian rate) that can spread rapidly through a population even if they reduce fitness. HEGs have potential to drive introduced traits through a population without large-scale sustained releases. The population genetics of HEG-based systems has been established using discrete-time mathematical models. However, several ecologically important aspects remain unexplored. We formulate a new continuous-time (overlapping generations) combined population dynamic and genetic model and apply it to a HEG that targets and knocks out a gene that is important for survival. We explore the effects of density dependence ranging from undercompensating to overcompensating larval competition, occurring before or after HEG fitness effects, and consider differences in competitive effect between genotypes (wild-type, heterozygotes and HEG homozygotes). We show that population outcomes—elimination, suppression or loss of the HEG—depend crucially on the interaction between these ecological aspects and genetics, and explain how the HEG fitness properties, the homing rate (drive) and the insect's life-history parameters influence those outcomes.


Introduction
Molecular biology tools are facilitating new genetics-based, species-specific methods of controlling insect pest populations, with the intention of improving efficacy over current methods and reducing adverse environmental consequences. In public health, the new technologies aim to reduce the burden of vector-borne diseases; key targets are the mosquito species Anopheles gambiae (major vector of human malaria, which causes 150-270 million cases and 0.6-1.6 million deaths annually [1,2]) and Aedes aegypti (principal vector of dengue virus, with ca 50-100 million infections and 18 000-19 000 deaths a year [3,4]). Transgenic approaches, involving the release of genetically modified (GM) mosquitoes into natural populations, include a broad class of population reduction methods aiming to suppress the numbers to a lower level or possibly local elimination, based on the principles of the sterile insect technique (SIT) [5]. Released GM mosquitoes mate with wild mosquitoes, potentially affecting the population genetics and population dynamics of the natural population, so ecological understanding is important if these vector control approaches are to improve human health [6,7].
Homing endonuclease genes (HEGs) are 'selfish' genes that can spread rapidly through populations even if they harm the host organism, because they are inherited at a higher than Mendelian rate [8]. HEGs exploit cellular repair mechanisms to copy themselves. A HEG encodes an endonuclease that recognizes and cuts a very short, specific DNA sequence (about 15-30 bp). The HEG sits inside its recognition sequence, so, in a heterozygous individual, the target sequence on the homologous chromosome is cut, and the cell's DNA repair machinery may use the HEG-bearing (HEG þ ) chromosome as a repair template, which results in conversion of a HEG heterozygote to a HEG homozygote. Engineered HEGs have potential to drive introduced traits, such as sterility or inability to transmit disease, through a mosquito population without needing large-scale sustained releases. Use of a meiosis-specific promoter to control the HEG allows heterozygous individuals to develop normally, but causes biased transmission of the HEG to their gametes. Two strategies are under development in Anopheles mosquitoes [9]. Releasing males carrying a HEG construct that targets a gene essential for female fertility could substantially reduce a target population of mosquitoes in similar manner to SIT. Males carrying an X-targeting HEG construct on their Y chromosome will produce a male-biased sex ratio, reducing both the number of mosquitoes (fewer females laying eggs) and disease transmission (fewer females biting). Some transgenic lines that were intended to target Xbearing sperm in this manner, actually created genetically sterile males ( protein transferred to the zygote with the sperm resulted in damage to the maternal X chromosome, causing embryo lethality) [10], and those mosquitoes are in early stage trials [11]. As part of a broader programme of research on the ecology and genetics of insect control, here we formulate a population ecological and genetic model and apply it to the simple case of a HEG, present in both sexes, that targets a gene essential for survival, thereby causing pre-adult lethality.
The population genetics of HEG-based systems has been established using discrete-time mathematical models [12,13]. However, this can only form part of our understanding of this approach as a vector control method, as several ecologically important aspects remain unexplored. The goal is fewer mosquitoes, so population dynamics must be combined with the population genetics. Here, we address three key population dynamic issues: overlapping generations, density-dependent larval competition and the timing of the HEG's effect on survival relative to density-dependent mortality.
Little is really known about intra-or interspecific competition among most species of vector mosquitoes, however, there is growing evidence that competition gives rise to density-dependent survival through the larval stages [14][15][16][17][18][19][20][21][22]. Although it is not thought to be common in natural populations, overcompensatory competition has been demonstrated in a variety of insect species [23][24][25][26], including Anopheles arabiensis [25] and one of several possible model fits to Ae. aegypti field data [14,15]. Recent work [27], coupling a HEG population genetics model with simple dynamics, has shown that HEG-based vector control can potentially suppress or eliminate a mosquito population in timescales that are practical for disease control; combination with simple malaria epidemiology suggests that there is little scope to eliminate disease transmission without also eliminating the mosquito population. Those results were generated using a density function that imposed extra larval mortality that increased monotonically with larval density, and did not allow for the possibility of overcompensating larval competition. Among other aims, here we extend this finding to explore a wider range of mechanisms of density-dependent competition. In view of the lack of empirical evidence, we consider a broad range of strengths of density dependence to examine the effects.
Mathematical modelling of other genetic strategies predicts that late-acting genetic lethality is better at controlling populations that are restricted by density-dependent larval mortality [27][28][29]. We explore fully and test the hypothesis that HEG mortality acting after that competition should be more effective for population control than early-acting HEG mortality.
Mosquitoes do not exhibit synchronized discrete generations, but overlapping generations occur with all life stages (eggs, various instar larvae, pupae and adults) coexisting in natural habitats. Here, we formulate and analyse acontinuoustime model to explore both mosquito population genetics and population dynamics. We adapt little-known work by Kostitzin [30], whose method applies Lotka-Volterra style competition equations (with linear competition based on the logistic equation) to different 'groups' (genotypes) rather than different species. We use a fixed time delay to represent density-dependent competition in immature stages, and extend Kostitzin's approach to incorporate a flexible twoparameter nonlinear intraspecific competition model which allows us to investigate the effects of under-or overcompensating density-dependent population regulation and to vary the competitive ability (as well as survival) between genotypes. We apply this set of delay differential equations to the release of mosquitoes bearing an early-acting or late-acting HEG that targets a gene important to survival, and investigate the interplay between population genetics and population dynamics. We provide a richer understanding of the population genetics and show that the potential efficacy of vector control depends crucially on the ecology ( particularly population dynamic regulation) and the genetics.

Material and methods
We set out the state variables and parameters in table 1 and the  full mathematical model in table 2. This is a deterministic model, assuming a panmictic (random mating) population with no genetic mutation or genetic drift. A fuller, more detailed derivation of the model is given in the electronic supplementary material, section Methods.
N i denotes the number of adults of genotype i (i takes values 1: ww, 2: Hw, 3: HH, where H is HEG þ and w is HEG 2 wildtype). These variables are real numbers (not integers), in effect assuming a large population, and are expressed in units of adult mosquitoes per host. Adult mortality occurs at rate m, and intrinsic per capita growth rate r (net of density-independent mortality during immature stages) and m are identical for all types. Population growth terms rN i are replaced by rn i , where the progeny arising n i are functions that reflect mating crosses among the three genotypes [30], and incorporate the effects of homing. The homing rate g is the proportion of successful conversions of w (HEG 2 ) to H (HEG þ ) in heterozygous individuals during meiosis. Density-dependent competition among larvae is reflected in the equations, with some simplifying assumptions, using a fixed time delay t [36]. The composition and number of new adults emerging at time t depend on the number and type of adults that were mating and laying eggs at a previous time t 2 t.
We formulate two versions of this model (table 2): (1) 'early-acting' HEG fitness costs, where HEG mortality occurs after homing and before density-dependent competition effects-e.g. the knockout renders fertilized eggs non-viable; or (2) 'late-acting' HEG fitness costs, where HEG mortality occurs after density-dependent competition and before mating-e.g. the target gene is essential for pupation.
Relative fitness parameters c 1 ¼ 1, c 2 ¼ (12hs), c 3 ¼ (12s) incorporate HEG-induced mortality, such that homozygotes suffer a fitness penalty (s . 0) as a result of the target gene being knocked out and, unless that knockout is recessive (h ¼ 0), heterozygotes incur a partial fitness penalty (hs where 0 h , 1). For each genotype, the net population growth term is multiplied by the relative fitness c i . This applies to both early-and lateacting HEGs. However, early-and late-acting HEGs have a different impact on larval competition (see below), because of the timing of the induced mortality. Bellows [23] examined several density-dependent competition functions of different forms and assessed their ability to represent 30 datasets exhibiting a range of intraspecific population dynamic behaviours. We adapt the nonlinear competition function that Bellows judged most flexible, based on that of Maynard Smith & Slatkin [37]: All genotypes experience the same strength of density-dependent intraspecific competition (which can range from contest to scramble by varying the value of coefficient b, with b . 1 denoting overcompensating scramble-like competition). The density-dependent scale parameters are either a single coefficient a or are genotype-specific a i . This allows us to explore the effects where the HEG can alter the competitive performance of transgenic larvae and not just the proportion surviving to reproductive maturity; a different value of a i could be thought of as representing different resource needs (perhaps transgenic individuals are smaller or have different metabolism from wild-type), and therefore altering the scale at which density-dependent regulation acts (and the number of that genotype that the habitat is capable of supporting). Mortality from a late-acting HEG occurs after the larvae have participated in density-dependent competition, so it does not appear in the competition term, i.e. the denominator of the fraction. An early-acting HEG will already have reduced the numbers of larvae, so the survivors suffer less competitive conditions, with the relative fitnesses (c i ) reducing the effect of competition. The gametic frequency of the HEG at time, t, is determined after homing: we label the frequencies of H HEG þ (q) and w HEG 2 ( p) sperm and eggs produced by adults at time t, where p(t) þ q(t) ¼ 1. GM insects are released into the population at time t ¼ 0; the gametic frequency at that time is denoted by q 0 . Mortality of those released adults reduces the gametic frequency to q t by the time the first transgenic progeny emerge as adults. We numerically simulated a one-off release of a number of adults equal to 25% of the equilibrium of the natural population (0.25N*), all of them heterozygotes (Hw), giving an instantaneous HEG allele frequency of 0.1. (For any given homing rate g, the gametic frequencies q 0 and q t can then be calculated; for our default homing rate value, g ¼ 0.8, gametic frequencies are q 0 ¼ 0.18 and q t ¼ 0.0304.)

Population genetics
Our population genetic outcomes are consistent with previous work [13], albeit with one new element resulting from taking account of the biological time lag (q t substitutes for q 0 where gametic frequency is necessary to determine the outcome). Gametic frequencies q ¼ 0 (HEG not present or extinct) and q ¼ 1 (HEG reaches fixation) are always equilibrium solutions. An internal equilibrium q* (0 , q* , 1) exists if and only if a specified relationship between the genetic properties of the HEG (relative fitness costs of gene disruption s, dominance of fitness costs h and homing rate g) holds true. Where it exists q Ã ¼ gð1 À hsÞ À hs sð1 À 2hÞ : ð3:1Þ By boosting allele transmission rates, homing confers an advantage on the HEG. Homing occurs only in a heterozygous individual, during meiosis, where it converts (at rate g) a HEG 2 (w) gamete to HEG þ (H ). This results in reduced fitness (survival) of the resulting progeny. If the other parent contributes a w gamete, then the zygote that would have been wild-type (ww, no fitness penalty) becomes heterozygous (Hw, with lower relative fitness c 2 ¼ (1 2 hs)). This incurs a selection disadvantage relative to the fitness of the heterozygous parent in which homing occurred Similarly, if the other parent contributes an H gamete, then a zygote that would have been heterozygous (Hw, with fitness c 2 ¼ (1 2 hs)) is instead homozygous (HH, with lower fitness c 3 ¼ (1 2 s)) and the selection disadvantage relative to the fitness of the Hw parent in which homing occurred is In a synthetic HEG, intended to spread through a vector population, the engineered gene disruption effect would preferably be nearer recessive than dominant (h , 1/2). If so, then heterozygote fitness would be closer to that of wildtype than of HH homozygotes, and equation (3.2a) would be smaller than (3.2b). Otherwise, the converse is true. The population genetic outcomes, which are summarized in the upper part of figure 1, are a consequence of the relative Table 2. Model for genetic mosquito control strategy using a HEG affecting survival.
(1) early-acting HEG fitness cost (2) late-acting HEG fitness cost rsif.royalsocietypublishing.org J. R. Soc. Interface 11: 20131071 weight of these forces and so are dependent on the relationships between the homing rate g and the phenotypic (fitness) properties of the HEG (s, h). Figure 2 illustrates the four regions, in terms of s and h, where each of these population genetic outcomes occurs. If the drive advantage (g) is lower than both the selection disadvantages (equations (3.2a) and (3.2b)), then the HEG will go extinct, i.e. to gametic equilibrium q ¼ 0. In that case, q ¼ 1 is unstable, so any introduction of wild-type alleles into a HEG þ homozygote population should eventually remove the HEG. If the genetic drive (g) is high enough that it outweighs both causes of reduced fitness (equations (3.2a) and (3.2b)), then the HEG will spread to fixation; the gametic equilibrium is q ¼ 1 (q ¼ 0 is unstable, so the HEG will always invade even from very low levels). Assuming low dominance of fitness effects (h , 1/2), an intermediate drive (g) can favour Hw individuals (equation (3.2a)) but be insufficient to outweigh the selection against HH genotypes (3.2b), and these balancing forces drive the population to a stable internal equilibrium consisting of all viable genotypes. In that case, the equilibrium gametic frequency of the HEG is given by q* (equation (3.1)); either allele can invade from rare (q ¼ 1 and q ¼ 0 are unstable). If the HEG is strongly deleterious to heterozygotes (h . 1/2), then an intermediate genetic drive (g) produces opposing selective forces that are frequency dependent. If the H allele is predominant in the population (q t . q*), then the HEG will be driven to fixation, because the drive (g) outweighs the fitness effect of replacing Hw progeny with HH (equation (3.2b)), whereas if the HEG is relatively rare (q t , q*), then it will become extinct, because homing does not overcome the disadvantage of replacing ww with Hw (equation (3.2a)). In those circumstances, neither allele can invade from rare (q* is unstable, q ¼ 1 and q ¼ 0 are stable).
Note that the gametic equilibrium reached is independent of initial gametic frequencies q 0 and q t except for that last case, i.e. a HEG that has strong fitness effects (h . 1/2) and an intermediate homing rate (the relevant inequalities are satisfied with moderate-to-high fitness cost, s; figure 2 and the electronic supplementary material, section Population genetics outcomes). In comparison, in all cases, the initial gametic frequency does affect the time taken to reach the predicted equilibrium frequency.
One key feature of using a continuous-time framework rather than a discrete generation model is that it takes account of the biological time lag from egg to adult emergence. In the period from release of HEG-bearing insects (t ¼ 0) to the time when the first transgenic progeny emerge as adults (t ¼ t), the released GM insects reduce in number owing to mortality, but the wild-type adults are replenished by emergence of existing larvae. Consequently, in that period, the HEG gametic frequency is diluted from its instantaneous value at release to a frequency that may be significantly lower. With our default parameter values, gametic frequencies fall from q 0 ¼ 0.18 to HEG goes extinct equilibrium q = 0 population genetics outcome (depends on fitness properties s, h and homing rate g of HEG) initial (release) conditions q t < q* ?  Figure 1. Overview of outcomes. This summarizes the way in which mosquito population genetics and dynamics interact. The population genetic outcome depends on the fitness properties of the HEG (cost s, and its dominance h) and the homing rate (g). The population dynamic implications depend on the relationship between the life-history parameters (r ¼ r/m lifetime progeny net of density-independent juvenile mortality, and t egg-to-adult generation time) and the HEG parameters, particularly the HEG load (0 for ww, s in HH population, or L (equation (3.3)) at mixed genotype equilibrium). rsif.royalsocietypublishing.org J. R. Soc. Interface 11: 20131071 q t ¼ 0.0304. In circumstances where the gametic frequency determines the population genetic outcome, it is very likely that the HEG allele will be sufficiently rare (q t , q*) that the HEG will be lost from the population rather than driven to fixation. The quantity of insects released would have to be substantially higher to mitigate that effect. (See the electronic supplementary material, section Population genetics outcomes.) The 'HEG load' (analogous to genetic load in classical population genetics) has been defined as the relative reduction in the population growth rate of the population in the presence of the HEG [13]. The HEG load in a wildtype population is zero, in a pure HH population it is s, and at the mixed genotype equilibrium q* the HEG load L is which, alternatively, can be expressed as 1 À L ¼ ð1 À g 2 Þð1 À hsÞ 2 À ð1 À sÞ sð1 À 2hÞ : ð3:3bÞ At the mixed genotype equilibrium, the HEG load is the sum of fitness penalties weighted by their genotype frequency (equation (3.3a)) and is influenced by the homing rate as well as the relative fitnesses of HEG homozygotes (1 2 s) and heterozygotes (1 2 hs) (equation (3.3b)).

Population dynamics
The natural wild-type population (time-delayed version of equation (2.1)) has the following equilibrium We assume that r . m (otherwise, the wild-type population would naturally die out). We fix this natural equilibrium at 60 adult mosquitoes per host (table 1)   is larger than the hypothetical steady state; the stronger the density dependence, the greater the amount by which N exceeds N*.
The effect of releasing HEG-bearing insects into the population is to introduce a HEG load. If the HEG spreads, then this reduces the reproductive ability of the population (in effect, a scaling factor is applied to r). Where the population is suppressed but not eliminated, this effect is more likely to give stable dynamics at the post-release equilibrium, as the stability boundaries correspond to higher values of b, the strength of density-dependent competition (the electronic supplementary material, section Population dynamics: stability). If the HEG goes extinct, then the population reverts to its natural equilibrium (equation (3.4)) or oscillatory behaviour. To predict these outcomes and, in particular, to ascertain the properties that are expected to lead to local elimination of the vector population and thereby eliminate disease transmission, it is necessary to explore the combined effects of population genetics and population dynamics.

Population dynamics combined with genetics
If the HEG's genetic properties allow it to spread through the vector population to some extent, then it disrupts the target gene and reduces survival. If a sufficiently high HEG load can be imposed, then the population can be locally eliminated. A less effective HEG load will suppress the population to a lower abundance which may (or may not) be significant. These effects can be achieved in timescales that are potentially useful for disease control. (See figure 3, illustrating a recessive lethal knockout.) As expected, the late-acting version has a lower post-release equilibrium than early-acting (figure 3 compare bottom row with top). An early-acting HEG can actually increase the population (figure 3c,e). By killing eggs or early instars, the early-acting HEG system creates a lower-density environment in which unaffected individuals and survivors will compete, and with overcompensating (scramble-like) competition that population rebounds to a higher abundance. This effect is similar to that already known for classical and genetics-based SIT [39,40].
For an overview of results, see figure 1. Population elimination ðN Ã 1 ,N Ã 2 ,N Ã 3 Þ ¼ ð0,0,0Þ is an equilibrium for both early-and late-acting HEG constructs. It is achieved when the net population growth rate r multiplied by 1 -HEG load is below unity, otherwise the population will persist. Where the HEG properties are such that the HEG will go to fixation, the equilibrium is 0; 0; 1 a 3 ð1 À sÞ rð1 À sÞ m À 1 1=b ! early-acting or ð3:5aÞ where r(1 2 s) . 1; otherwise, a population of HEG homozygotes would die out. Where a mixed genotype population occurs, and r(1 2 L) . 1, the equilibria satisfy the following (q* as per equation (3.1), p* ¼ 1 2 q*, and L per equation (3.3)).
N Ã ¼ ½rð1 À LÞ À 1 1=b a 1 p Ã2 þ 2a 2 p Ã q Ã ð1 À hsÞ þ a 3 q Ã2 ð1 À sÞ early-acting, or ð3:6aÞ N Ã ¼ ½rð1 À LÞ À 1 1=b a 1 p Ã2 þ 2a 2 p Ã q Ã þ a 3 q Ã2 late-acting ð3:6bÞ Equations (3.5a,b) and (3.6a,b) confirm that any persistent late-acting HEG achieves a lower population equilibrium than the equivalent early-acting HEG (the electronic supplementary material, section Population dynamic outcomes). If density-dependent competition is very strong (high b), then the population will oscillate around the relevant equilibrium, with average value ( N) higher than that steady-state value (N*). A late-acting HEG achieves a lower average abundance than the equivalent early-acting HEG. The stability boundaries in a HEG-bearing population occur at higher values of b than for the natural population (the electronic supplementary material, section Population dynamics: stability). A population suppressed by HEG vector control is very likely to exhibit stable dynamics (e.g. figure 3 and the electronic supplementary material, table S2). The greater the HEG load imposed, the more stable the system is, and the more resilient it is to small perturbations.
In the mixed genotype case, increasing the homing rate g increases the HEG load (inspect equation (3.3b)) and so generally decreases the post-release population equilibrium (N*, equations (3.6)) and, if unstable (very high values of b), decreases the average density ( N). The exception is an early-acting HEG in a population with overcompensating density dependence (b . 1), where the HEG 'load' benefits the larval population, increasing the average adult abundance. These findings are illustrated for recessive lethal knockouts (s ¼ 1, h ¼ 0; figure 4). A recessive lethal HEG is always driven to an intermediate equilibrium frequency q* ¼ g (equation (3.1)), whatever the homing rate (equation (3.2a) equals 0 and equation (3.2b) equals 1), and HEG load L ¼ g 2 (equation (3.3b)). In this case, the homing rate above which r(1 2 L) , 1 and the population is eliminated is With our parameter values (table 1), this threshold value for a recessive lethal knockout is g ¼ 0.942.
The HEG fitness properties, (s, h), affect the equilibrium (N* where stable; figure 5) or average ( N where unstable) population density, through the four underlying population genetics cases (figures 1 and 2). With homing rate g ¼ 0.8, for most combinations of s, h, the HEG spreads to fixation, and the equilibrium density (equations (3.5)) decreases as the HEG-induced mortality s increases ( figure 5a,b,d). With an early-acting HEG, a population under strong densitydependent competition rebounds to higher than natural densities (figure 5c). Where fitness cost (s) is high and nearer recessive than dominant (h , 1/2), the population reaches a mixed genotype equilibrium giving intermediate suppression ( figure 5, lower right of each panel a,b,d). Where s and h are both high ( figure 5, upper right of each panel a-d), the HEG is driven extinct, and the population returns to its natural equilibrium. For much of that region (except a triangular area at the rsif.royalsocietypublishing.org J. R. Soc. Interface 11: 20131071 high s end; figure 2), the outcome might be altered by releasing far more transgenic insects, so that the HEG is sufficiently common to tip the balance to HEG fixation instead of extinction. Even with this high homing rate (0.8) and a fairly large release (augmenting the original population by one-quarter), there are few combinations of s and h that suppress the population to a level at or close to local elimination. These combinations occur when s , 2g/(1 þ gÞ and h , 1/2, close to the point ðs; hÞ ¼ ð2g=ð1 þ gÞ; 1/2Þ; where the four regions of population genetic outcomes intersect (figure 2); the HEG is at fixation. The effect of differential competition among genotypes is explored by varying u (. or ,0) to alter the relative competition effect of HH genotypes and f, the dominance of that impact (table 1; a 3 ¼ (1 þ u)a 1 and a 2 ¼ (1 þ fu)a 1 ). This changes the equilibrium N* (and the average N, if unstable) in a suppressed population. Relative to the result with equal competition, the effect of unequal competitive effects is to multiply the equilibria by: ; HEG fixed ð3:8aÞ 1 À L 1 À L þ 2wup Ã q Ã ð1 À hsÞ þ uq Ã2 ð1 À sÞ; mixed genotypes, early-acting ð3:8bÞ and 1 À L 1 À L þ 2fu p Ã q Ã þ uq Ã2 ; mixed genotypes, late-acting.

ð3:8cÞ
If the effect is dominant (f ¼ 1), then the expressions (3.8b) and (3.8c) collapse to ð1 À LÞ/½1 À L þ uð1 À L À p Ã2 Þ and 1/½1 þ uð1 À p Ã2 Þ, respectively. If transgenic insects require greater resources and therefore exert competitive pressure at lower scales (u . 0), then the suppressed population will settle to a lower density than if all genotypes had competed equally (and the opposite for u , 0). The stronger the dominance of this alteration to competitive effect (higher f ), the larger the impact on population size (illustrated in figure 6). We explored a comprehensive range of values of s, h, g, b and a i . Our findings are also robust to changes in t, r and m (the electronic supplementary material, section Population dynamics outcomes).

Discussion
We have devised a novel mathematical modelling method for combining population genetics and population dynamics in a continuous-time framework with density-dependent competition. We applied it to a genetic vector control strategy using HEGs to drive a gene knockout through a population, inspired by a system that is in development and progressing towards implementation. We have shown that across a range of possible parameter sets, the outcome is generally more likely to be population suppression than to be local elimination. There have been few ecological studies of the population dynamics of A. gambiae mosquitoes, of which very few were conducted in semi-natural conditions [17,41]. Given our poor understanding of population regulation in malaria vectors, our results highlight the variation in possible outcomes of HEG vector control under different assumptions about the strength of density-dependent larval competition.
The influences on population genetics behaviour can be separated out between the HEG's fitness properties (the penalty to survival, s, and the dominance of that fitness cost h) and the homing rate g. Engineering a synthetic HEG targeting a gene that is important for survival and creating transgenic insect lines with that HEG positioned within its recognition site is challenging, and fitness properties, homing rates and regulation of gene expression cannot be precisely controlled [42]. A homing rate g . 0.942 is needed for a recessive lethal knockout to locally eliminate a population (equation (3.7)). Given the homing rates achieved in experimental strains, 0.56 to approximately 0.9 [10,35], such a high rate presents a challenging target. It has been proposed that to achieve a high enough HEG load in practice with realistic homing rates, multiple HEGs could be used, each targeting a different gene, with independent fitness effects across loci [9]. More sophisticated HEG systems are being developed for field use, targeting female fertility or imposing a male bias, and it is intended to combine multiple copies in a released strain to increase the likelihood of control success.
If the HEG persists within the population (at fixation or some intermediate frequency), then the relationship between the genetic parameters (s, h and g within HEG load L, or simply the HEG load s at fixation) and the life-history traits (r) determine whether the mosquito population is locally  eliminated or merely suppressed. The higher the reproductive potential of the population (r), the greater the HEG load needed to overcome it. For example, the conditions for population elimination are 1 2 s . 1/r for a population in which the HEG goes to fixation, or g . ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 À 1=r p for a recessive lethal knockout (equation (3.7)). Even in these simpler  rsif.royalsocietypublishing.org J. R. Soc. Interface 11: 20131071 cases, it is necessary to understand the population dynamics and not just the population genetics, to predict whether local elimination is feasible. Close to thresholds, elimination can take a long time to achieve. Our model is restricted to consider only a single value of b; all genotypes experience the same strength of densitydependent competition. This is a mathematical consequence of using nonlinear population dynamic models of this kind (electronic supplementary material, section Methods). We have modelled a wide range of strengths from very weak (undercompensating competition with low b) to very strong (overcompensating with high b), to explore the range of possible effects. Our illustrations present results generated with some parameter values for b that may be regarded as implausibly high or low, not to suggest that such extremes are to be expected, but for ease of visibly discerning the nature of the effects identified across the spectrum of strength of density dependence.
One novel consideration in our study is altering the competitive impact of larvae bearing one or two copies of the HEG transgene. The density scale parameter in this model is manifested as a competitive effect (a change affects larvae of all genotypes) rather than as a competitive response. If the parameter is higher for HEG homozygotes than wild-type insects (a 3 . a 1 ), then HH larvae contribute greater competitive pressure on all genotypes. This might plausibly occur, for example, if HH individuals process food less efficiently and need to consume more nutrients during development, reducing the nutritional resources available to larvae of all types. As the transgenic insects have greater resource utilization needs, the environment would be able support fewer of them. By contrast, an example of different competitive response might be less ability to use an alternative food resource to ameliorate some of the pressure. More sophisticated models would be required to tease out such effect/response distinctions. In our model, the density scale parameters do not affect whether a population is eliminated or persists following HEG release, nor the equilibrium mix of genotypes. They do affect the population size, by changing the scale at which density-dependent effects have significant impact, and so would affect the disease transmission capability of the vector population following release.
Where the population can be suppressed but not eliminated, as is the case across most of the parameter space for properties of the HEG even with our generous homing rate 0.8 ( figure 5), a crucial question is whether the new equilibrium is below the entomological threshold necessary to sustain disease transmission. As estimates of the basic reproduction number of malaria (R 0 , or Z 0 ) across Africa vary from near 1 to over 3000 [43], broad general assertions are inappropriate, but some patterns can be observed. If the wild population is under strong density-dependent competition with oscillatory dynamics, then HEG suppression of the population would have a dual effect on population size, by both inducing stable dynamics (the steady-state N* is lower than the average N) and reducing the equilibrium size, potentially increasing the chances of a successful epidemiological outcome. A HEG that affects both survival and competitive ability (increased density scale parameter) can be more effective at population suppression than a HEG that impacts on survival alone, which also enhances the potential for disease reduction. Depending on the extent of suppression, stochastic effects (which are not included in our deterministic framework) could become relevant, as the mosquito population could be eliminated at low numbers. In principle, a deterministically persistent HEG might become extinct in a small stochastic population owing to genetic drift.
As we and others have noted [15,16], there is a lack of data about density-dependent competition in the field for mosquito disease vector species. The few experiments published have almost all studied container-dwelling Aedes mosquitoes, where overcompensating competition is thought to be possible, and almost nothing is known about Anopheles mosquitoes (an exception being ref. [17]), which will be critical for applying genetic control tools to malaria vectors. Methods have been demonstrated for studying natural container-breeding populations in natural settings (without artificial addition of water or food materials), varying mosquito densities across natural ranges [21,22]; more of such manipulative ecological studies to determine how density dependence acts and to test for compensatory effects of added mortality would have long-term value, especially if they can quantify the form and shape of density dependence functions for modelling genetic vector control. This is one of several challenges in mosquito ecology, another is estimating the intrinsic rate of population increase. Addressing these ecological challenges could make a real contribution to developing, regulating and implementing feasible vector control methods. For example, regulators might be concerned at the theoretical prospect of early-acting (but not late-acting) HEG releases possibly increasing a vector population that is regulated by overcompensating competition, although there is no evidence of overcompensating density-dependent larval competition in A. gambiae in the field and the potential risk might be unfounded. Incorporating stochasticity and spatial and temporal heterogeneity into models, might also dampen the propensity for strong oscillations in population dynamics.
Although we applied our method to a particular application of HEGs to mosquito control, this approach has much broader relevance. It can be extended to other HEG strategies (for example, a HEG targeting female fertility can be modelled by expanding the system with separate equations for males and females) and, in principle, to other genetic vector control strategies. Our framework could be expanded to include a second insect species, with intraspecific and interspecific competition following the same functional form (with different resource utilization) [44]. Further aspects of mosquito ecology could be incorporated, such as seasonal fluctuations in population size or extension to structured populations with limited interbreeding. Entomological models could be coupled with epidemiological models to explore the potential for alleviating human disease. This would be necessary for assessing strategies that aim to propagate a genetic refractoriness to disease transmission through the vector population. Extension to include a health economic analysis could be used to assess the consequent reduction in disease burden [45].
Previous models of selfish genetic elements used discretetime population genetic frameworks to deduce mathematical conditions for invasion of the novel element primarily in terms of biased gene conversion (transmission) and reduced survival (fitness) [9,12,13,46]. Those papers that explicitly considered vector population control applications explored outcomes in terms of allele or gametic frequencies, either with no direct modelling of the effect on population size [9,12,13], or with fairly simple population dynamics [27]. Our study investigated the effects of the opposing forces of drive and natural selection in a richer population dynamic framework that is able to represent a wide range of densitydependent larval competition dynamics, the relative timing of key biological processes, and fitness effects that alter competitive ability as well as reduce survival. Our key novel findings include that: the HEG gametic frequency is significantly diluted between the release of insects and the adult emergence of their first transgenic progeny, during which time further wild-type larvae emerge as adults (with practical implications for how many insects to release); a persistent HEG has a stabilizing effect on population dynamics, which could further suppress the mean vector abundance (for more effective disease control); an early-acting HEG performs worse than a late-acting HEG and would be counterproductive in the presence of overcompensating density-dependent larval competition; a HEG altering competitive effect as well as survival can give more effective vector control. We have shown that mosquito control outcomes depend on the interplay between genetics and ecology (here focusing on overlapping generations, larval competition and the maturation time delay from egg to adult) and argue that both aspects need to be appropriately considered when assessing the potential effectiveness of genetic vector control methods.