Skip to main content

Fine-scale estimation of key life-history parameters of malaria vectors: implications for next-generation vector control technologies

Abstract

Background

Mosquito control has the potential to significantly reduce malaria burden on a region, but to influence public health policy must also show cost-effectiveness. Gaps in our knowledge of mosquito population dynamics mean that mathematical modelling of vector control interventions have typically made simplifying assumptions about key aspects of mosquito ecology. Often, these assumptions can distort the predicted efficacy of vector control, particularly next-generation tools such as gene drive, which are highly sensitive to local mosquito dynamics.

Methods

We developed a discrete-time stochastic mathematical model of mosquito population dynamics to explore the fine-scale behaviour of egg-laying and larval density dependence on parameter estimation. The model was fitted to longitudinal mosquito population count data using particle Markov chain Monte Carlo methods.

Results

By modelling fine-scale behaviour of egg-laying under varying density dependence scenarios we refine our life history parameter estimates, and in particular we see how model assumptions affect population growth rate (Rm), a crucial determinate of vector control efficacy.

Conclusions

Subsequent application of these new parameter estimates to gene drive models show how the understanding and implementation of fine-scale processes, when deriving parameter estimates, may have a profound influence on successful vector control. The consequences of this may be of crucial interest when devising future public health policy.

Graphic abstract

Background

Despite considerable progress, the global burden of disease caused by malaria remains high, and the demand for new, highly effective interventions is urgent [1]. To help reduce transmission, next-generation methods of control utilising genetic modification are being developed, which aim to reduce either the density or competence of vector populations [2]. With the release of genetically engineered organisms being strictly regulated, mathematical models remain a key resource for understanding their potential. However, gaps in our knowledge lead to an oversimplification of fine-scale mosquito dynamics, which can have a significant impact on our predictions of vector control efficacy.

The malaria parasite is transmitted principally by mosquitoes within the Anopheles gambiae complex, which comprises eight morphologically indistinguishable species [3]. In Africa, where 90% of malaria cases occur [1], Anopheles gambiae sensu stricto and Anopheles funestus are the dominant vectors in the majority of regions [4, 5], with Anopheles arabiensis and Anopheles coluzzii also being of significant importance [6]. Although each species exhibits some variation in behaviour and environmental tolerance [7,8,9,10], larval requirements for aquatic habitats mean that mosquito densities tend to peak at the times of year associated with high rainfall [11,12,13]. Interventions to reduce vector populations during these periods are an established method for reducing transmission—including larvicides [14], insecticide-impregnated bed nets [15] and pyrethrum spraying or trapping of adults [16]. These methods, however, have limitations, including the intensive use of chemicals, short half-lives of insecticides within the environment and the costs of continued redeployment. Next-generation “gene drive”-based interventions have the potential to overcome many of the problems inherent in traditional vector control [17]. They insert novel genes into a vector species which code for specific traits likely to cause either population suppression or reductions in vector competence, in a manner that allows those genes to propagate through the vector population. A promising technology in this field is the use of homing endonuclease genes (HEGs) [18], a class of nuclease genes found in simple single-celled organisms which can copy themselves from one chromosome to another. HEGs can be placed at specific sites on a chromosome and in a heterozygous individual produce an enzyme which cuts the DNA on the homolog of the HEG-bearing chromosome. When this site is repaired, the HEG-bearing chromosome is used as a template [19]. Such interventions have the potential to significantly reduce malaria transmission [2, 20,21,22], but are also susceptible to the emergence of resistance [23] or reinvasion of wild-type into depopulated landscapes [24].

Modelling shows that the spatial invasion dynamics of gene drive systems are highly sensitive to the fine-scale dynamics of mosquito populations [20]. Uncertainty around key aspects of these dynamics makes reliable predictions of the impact and spread of gene drive technologies challenging. In previous studies, the reproductive potential of a mosquito population (typically denoted as Rm) has been found to be a key driver of the efficacy of gene drive constructs [20, 21, 24, 25]. Rm is the maximum rate at which a population can grow in the absence of density-dependent constraints (i.e. in the infinitesimal population density limit). However, since real populations are discrete and clumped, and subject to differing regimes of resource competition, the theoretical limit represented by Rm may never be attained. Understanding how a more realistic representation of the mosquito life cycle affects estimates of Rm is therefore important for refining assessments of the likely impact of gene drive technologies.

Although not exclusively, Rm in mosquito populations is often limited principally at the aquatic stages by density-dependent larval mortality, characterised by two parameters: the baseline density-independent mortality rate and the carrying capacity of the local environment for mosquito larvae. Carrying capacity describes the maximum number of mosquito larvae that can be sustained by the resources available within the environment, whilst density dependence affects the mortality (and developmental delays) caused by intraspecific competition for these resources [26, 27]. The two are intrinsically linked—for example, decreasing carrying capacity would lead to an increase in resource competition and hence density-dependent effects. How mathematical models approach larval density dependence is currently limited by the few experimental studies that exist, and analysis of experimental data has shown both linear and quadratic relationships between density and mortality [26,27,28]. Recent modelling studies have accounted for a lack of empirical data by implementing flexible density dependence functions in models to simulate a range of scenarios [29, 30].

Although a number of approaches have been used to simulate density dependence, most models tend to assume that in the low-density limit (e.g. one mated female entering an environment otherwise empty of mosquitoes), density-dependent effects will not apply, and the population will grow at its theoretical maximum (i.e. Rm). However, this is likely to overestimate Rm, as it is inconsistent with how a mosquito lays eggs in the natural environment, where habitat heterogeneity, local environmental conditions and behavioural traits will cause females to lay eggs in temporally and spatially clumped batches [31]. This behaviour has been studied in both Aedes and Anopheles mosquitoes, in laboratory conditions [32,33,34] and in the wild [35,36,37].

In this paper, using a discrete-time stochastic mathematical model of mosquito population dynamics, we aim to identify the effects of egg “clumping” on the estimation of Rm. Additionally, as any conferred density effects from clumping are intrinsically linked with the functional form of density dependence, we explore multiple density dependence regimes. From this we aim to identify whether clumped egg-laying and specific forms of density dependence provide a better model fit to observed data, and the implications associated Rm estimates have on gene drive technologies.

Methods

Mosquito population model

We developed a discrete-time stochastic model that simulates the life cycle of a mosquito population. A schematic of the model is outlined in Fig. 1 and is described below.

Fig. 1
figure1

Model schematic showing transitions between stages, points of mortality and density-dependent mortality, immigration and egg-laying. Each arrow represents a Poisson draw, binomial or series of binomial events based on fitted probability parameters

The life cycle of a mosquito consists of four main stages: egg, larvae, pupae and adult. The first three stages are aquatic and, although there is variation in the wild, last approximately 5–14 days [38]. Larvae emerge from eggs and feed on algae, bacteria and other small particulate matter whilst developing through four separate larval instars, before metamorphosing into pupae and then finally adults. After maturation into the adult stage and having successfully mated, a female begins its first gonotrophic cycle. During the gonotrophic cycle, which lasts approximately 3 days, the mosquito takes a blood meal from a human or animal host, after which the blood is digested during a resting phase and used to develop the eggs in the ovaries. Once the eggs have matured, the mosquito is considered gravid and will search for a suitable site for oviposition. This process is then repeated until the eventual death of the adult.

As in previous models [28], to reduce complexity but maintain biological realism, we group the first two larval stages into a single class, early instar larval stage E, and the third and fourth stages into a late larval stage class L. Pupae P and adults M remain as single stages. Transition between early- and late-stage larvae and associated density-dependent mortality can be described as a binomial event occurring within discrete time \(\delta t\). The probability of mortality or transition can be defined as \(\delta t\) multiplied by the sum of the two respective rates, development to late-stage larvae (\(d_{E} )\) and density-dependent mortality \(\mu_{E} \left( {E,L} \right)\). A proportion \(\frac{{d_{E} }}{{\mu_{E} \left( {E,L} \right) + d_{E} }}\) of all the larvae transitioning out of E are expected to enter the state L, again simulated by sampling from the corresponding binomial distribution. Similarly, transitions between late-stage larvae and pupae can be described with a probability \(d_{E}\) and density-dependent mortality of \(\mu_{E} \left( {E,L} \right)\). Pupae then develop into adults with a probability of \(d_{P}\) and a density-independent probability of mortality \(\mu_{p}\). Half of all emerging adults are assumed to be female, and male adults are not explicitly modelled (we assume that there are enough males to successfully mate with all females). For simplicity and model tractability we assume that adults die at an age-independent rate \(\mu_{{\text{m}}}\) and have an approximate average life expectancy of 9–12 days in the wild [39]. Regular immigration of adult mosquitos from external populations is known to occur in most wild populations [40]; therefore, migration is modelled as a Poisson process with rate \(I\). The input of new larvae into the system comes from eggs laid at a constant rate, \(M\beta \delta t\), where \(\beta\) is a fitted parameter for the average daily egg-laying rate of adult females.

The dynamics of a mosquito population can thus be described as a series of stochastic equations where \({\text{Bin}}\left( {p,n} \right)\) denotes a binomial draw with \(p\) probability and \(n\) trials, and \({\text{Poisson}}\left( \lambda \right)\) a draw from a Poisson distribution with given lambda.

$$E\left( {t + \delta t} \right) = E\left( t \right) - B_{E} + M\beta \delta t,$$
$$L\left( {t + \delta t} \right) = L\left( t \right) - B_{L} + {\text{Bin}}\left( {\frac{{d_{E} }}{{\mu_{E} \left( {E,L} \right) + d_{E} }},B_{E} } \right),$$
$$P\left( {t + \delta t} \right) = P\left( t \right) - B_{P} + {\text{Bin}}\left( {\frac{{d_{L} }}{{\mu_{L} \left( {E,L} \right) + d_{L} }},B_{L} } \right),$$
$$M\left( {t + \delta t} \right) = M\left( t \right) + \frac{1}{2}{\text{Bin}}\left( {\frac{{d_{P} }}{{d_{P} + \mu_{P} }},B_{P} } \right) - {\text{Bin}}\left( {\mu_{M} \delta t,M} \right) + {\text{Pois}}\left( { I} \right),$$
(1)

where

$$B_{E} = {\text{Bin}}\left( {\left( {d_{E} + \mu_{E} \left( {E,L} \right)} \right)\delta t,E} \right),$$
$$B_{L} = {\text{Bin}}\left( {\left( {d_{L} + \mu_{L} \left( {E,L} \right)} \right)\delta t,L} \right),$$
$$B_{P} = {\text{Bin}}\left( {(d_{P} + \mu_{P} )\delta t,P} \right).$$

Density-dependent mortality

Density dependence in mosquito larval stages results from the competition for resources during the early and late instars; however, pupae are not pressured by resources, as they do not feed. For simplicity, although adults may have additional density constraints, it is assumed in this model that adults are not constrained by density effects [41].

Density dependence for early-stage larvae E(t) and late-stage larvae L(t) at time \(t\) can therefore be described as \(\mu_{E} \left( {E\left( t \right),L\left( t \right)} \right)\) and \(\mu_{L} \left( {E\left( t \right),L\left( t \right)} \right)\).

To quantify the effect of egg-laying under different potential systems of density dependence, three functional forms were considered within the model:

A linear relationship between larval density and density-dependent mortality, as used by White et al. [28]:

$$\mu_{E} \left( {E\left( t \right),L\left( t \right)} \right) = \mu_{E}^{0} \left( {1 + \frac{E\left( t \right) + L\left( t \right)}{{K\left( t \right)}}} \right)^{ } ,$$
(2)
$$\mu_{L} \left( {E\left( t \right),L\left( t \right)} \right) = \mu_{L}^{0} \left( {1 + \gamma \frac{E\left( t \right) + L\left( t \right)}{{K\left( t \right)}}} \right).$$

Density-dependent mortality, which grows as some (fitted) power (\(\Omega\)) of larval density, resulting in \(\Omega\) determining the intensity of density dependence, where \(\Omega\) < 1 intensity grows at a decreasing rate, whilst \(\Omega\) > 1 intensity grows at an increasing rate.

$$\mu_{E} \left( {E\left( t \right),L\left( t \right)} \right) = \mu_{E}^{0} \left( {1 + \frac{E\left( t \right) + L\left( t \right)}{{K\left( t \right)}}} \right)^{ \Omega } ,$$
(3)
$$\mu_{L} \left( {E\left( t \right),L\left( t \right)} \right) = \mu_{L}^{0} \left( {1 + \gamma \frac{E\left( t \right) + L\left( t \right)}{{K\left( t \right)}}} \right)^{\Omega } .$$

Exponential density dependence, where mortality increases exponentially with larval population density:

$$\mu_{E} \left( {E\left( t \right),L\left( t \right)} \right) = \mu_{E}^{0} e^{{\left( {\frac{E\left( t \right) + L\left( t \right)}{{K\left( t \right)}}} \right)}} ,$$
(4)
$$\mu_{L} \left( {E\left( t \right),L\left( t \right)} \right) = \mu_{L}^{0} e^{{\left( {\gamma \frac{E\left( t \right) + L\left( t \right)}{{K\left( t \right)}}} \right)}} .$$

Here, \(\mu_{E}^{0}\) and \(\mu_{L}^{0}\) are the death rates for early and late instars at low densities, \(K\left( t \right)\) is the environmental carrying capacity at time \(t\) and \(\gamma\) is the difference in effect of density dependence on late-stage instars compared with early-stage instars. Larval carrying capacity (K) was assumed to be proportional to rainfall from the previous \(\tau\) days weighted by an exponential distribution with mean \(2\tau\), where \({\text{rain}}\left( t \right)\) in millimeters is the daily rainfall and \(\lambda\) is a village-specific scaling factor as used by White et al. [28]:

$$K\left( t \right) = \lambda \frac{1}{{\tau \left( {1 - e^{{ - {\raise0.7ex\hbox{$t$} \!\mathord{\left/ {\vphantom {t \tau }}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$\tau $}}}} } \right)}}\mathop \int \limits_{0}^{t} e^{{ - \left( {t - t^{\prime}} \right)/\tau }} {\text{rain}}\left( {t^{\prime}} \right){\text{d}}t^{\prime}.$$
(5)

As there were inconsistent records of temperature in the Garki project for the study sites, we assumed a single fixed value. Similarly to White et al. [28], we used the mean air temperature of 24 °C to obtain an approximate mean value for the water temperature of 28 °C, and based our priors on this value.

Clumped egg-laying

To simulate clumped egg-laying, the input of eggs into the system (\(M\beta \delta t)\) in Eq. (1) is replaced by a stochastic “clumping” process. Here, the probability of each adult female laying eggs during any one time step \(\delta t\) is \(\frac{\delta t}{S}\), where \(S\) is the duration of the gonotrophic cycle. The total number of females laying eggs, \(E_{x}\), is thus sampled from a binomial distribution \(E_{x} = {\text{Bin}}\left( {\frac{\delta t}{S},M} \right)\). The number of eggs laid in that time step is then sampled from a Poisson with mean \(E_{x} n\), where \(n\) is the mean clutch size \({\text{Pois}}\left( {E_{x} n} \right)\). As \(n\), unlike \(\beta\), explicitly considers the gonotrophic cycle \(S\), the relationship between the two values is dependent on the average number of gonotrophic cycles within a mosquito's lifespan (\(\mu_{{\text{m}}}\)); therefore, approximately, \(n = \frac{{ \beta \left( {e^{{S\mu_{{\text{m}}} }} - 1} \right)}}{{\mu_{{\text{m}}} }}.\)

Data

The model was fitted to rainfall and mosquito catch data obtained from the Garki project [39]. The project was a World Health Organization-funded study into the effects of malaria and vector interventions during the 1970s and presents one of the most detailed longitudinal data sets for combined mosquito and meteorological data. The data consist of environmental, demographic and mosquito catch data from 24 villages over an approximately 5-year period. Time series of Anopheles gambiae caught by pyrethrum spray in houses were extracted for eight villages and the data pooled to the village level. From these villages, two were from control areas where no interventions were carried out after the first rainy season, allowing the data to be split into two successive seasons, giving a total of ten population time series.

Parameter estimation

We used an adaptive particle Markov chain Monte Carlo algorithm (pMCMC) to directly fit our stochastic model of mosquito dynamics to the data collected during the Garki project. pMCMC is a form of MCMC which uses a particle filter component to calculate a marginal likelihood for parameter acceptance/rejection within a Metropolis-Hastings algorithm [42, 43]. To account for the overdispersion and fractional sampling of the population, a beta-binomial likelihood function was used, with the parameter accounting for overdispersion \(r\) and the fraction of the population sampled \(p\) being fitted (see Additional file 1 for details of the pMCMC algorithm and likelihood methods).

The number of particles and chain length in a pMCMC has no set rules, and there is a trade-off between minimising the noise in the model and maintaining a realistic compute time. For each scenario of density dependence and egg-laying, we looked at the decrease in variation for set parameters when increasing the particle number and estimated the time to convergence in the chains. Therefore, the pMCMC chain was run with 150 particles for 2,000,000 iterations, with a 100,000-iteration burn-in. Attempts to increase the particle number and achieve convergence when fitting with an exponential density dependence were not successful, so full chains were not run. Priors were based on lab and field studies detailed in White et al. [28]. Approximate initial conditions for each of the villages were obtained by assuming deterministic equilibrium, allowing the value of each state variable to be derived from the value of a single estimate of the initial value of E, L, P, or M (see Additional file 1). This single parameter was denoted as z and estimated for each village. Most parameters were assumed to be common to all villages and were estimated by fitting the model to the data from all villages simultaneously, except village-/season-specific scaling factors \(\lambda_{{\left( {1 \ldots 10} \right)}}\) and initial condition parameters \(z_{{\left( {1 \ldots 10} \right)}}\), which were assumed to be village- and season-specific. Median values of the posterior distributions from the pMCMC chains for each parameter were then used to simulate the models for plotting, and 95% credible intervals were obtained by repeatedly sampling the joint posterior distribution.

Mosquito population growth rate

Rm in the context of this paper is considered as the population growth rate in the low-density limit, as quantified by the number of adult female offspring produced per adult female. To calculate this, two values were derived: a raw value in the absence of density dependence and a second value which included the effects of density dependence. To obtain a raw analytical value we followed White et al. [28], where a female laying \(n\) viable eggs during an oviposition cycle can expect to oviposit approximately \(n\left( {e^{{ - S\mu_{{\text{m}}} }} + e^{{ - 2S\mu_{{\text{m}}} }} + \cdots } \right) = \frac{n}{{e^{{S\mu_{{\text{m}}} }} - 1}}\) eggs, where \(e^{{ - S\mu_{{\text{m}}} }}\), \(e^{{ - 2S\mu_{{\text{m}}} }}\), etc., is the proportion surviving one and two gonotrophic cycles, respectively. The fraction surviving to adulthood can then be defined respectively for clumped and non-clumped egg-laying as:

$$R_{{\text{m}}} = \frac{1}{2}\left( {\frac{n}{{e^{{S\mu_{{\text{m}}} }} - 1}}} \right)\left( {\frac{1}{{1 + \mu_{{\frac{E}{{d_{E} }}}}^{0} }}} \right)\left( {\frac{1}{{1 + \frac{{\mu_{L}^{0} }}{{d_{L} }}}}} \right)\left( {\frac{1}{{1 + \frac{\mu P}{{d_{L} }}}}} \right),$$
(6)
$$R_{{\text{m}}} = \frac{1}{2}\left( {\frac{ \beta }{{\mu_{{\text{m}}} }}} \right)\left( {\frac{1}{{1 + \mu_{{\frac{E}{{d_{E} }}}}^{0} }}} \right)\left( {\frac{1}{{1 + \frac{{\mu_{L}^{0} }}{{d_{L} }}}}} \right)\left( {\frac{1}{{1 + \frac{\mu P}{{d_{L} }}}}} \right).$$
(7)

However, clumped egg-laying means larval populations always experience some density-dependent competition; thus the raw value for Rm never predicts actual population growth. Accounting for minimal density dependence is analytically challenging, however. We therefore evaluated Rm numerically by simulating the model for a single female mosquito until its death in a system empty of other adults and recorded the number of female offspring surviving until adulthood. Simulations were repeated using multiple random samples of parameter sets sampled from the pMCMC posterior to estimate credible intervals. To account for stochastic variation, the mean was taken from 50 simulations of each parameter set.

Effect on gene drive

To illustrate the effect on vector control, we looked at a test case of the impact on gene drive from changes in both density dependence and the addition of clumped egg-laying, by reconstructing the gene drive model formulated by Deredec et al. [20]. This model considers HEGs, a class of selfish genetic elements which are one of the most credible gene drive methods currently being developed. HEG efficacy is determined by the number of HEGs being used and their homing rate (the fraction of potential recipient heterozygous chromosomes that require the HEG). We ran the model using parameter estimates from our results for each of our scenarios and established the impact on the required HEGs and homing rate for elimination.

Results

Density-dependent mortality

The pMCMC chains gave good convergence, and the model visually fitted the village-level time series well, both when using a linear (Eq. 2) and fitted power (Eq. 3) density dependence (see Fig. 2 and Additional file 1). Final parameter estimates were taken as the median of the posterior estimates with 95% credible intervals obtained by repeatedly sampling parameter sets from the joint posterior, Table 1. It was not possible to fit an exponential density dependence (Eq. 4), with the pMCMC algorithm unable to obtain reliable posterior estimates for any parameters, suggesting that exponential density dependence is highly unlikely. Deviance information criterion (DIC) values show that fitting a power provided the best model, 1104.149, compared to 1115.039 for linear density.

Fig. 2
figure2

Model fits to Garki Project data for non-clumped egg-laying when fitting a power to density dependence with 95% credible intervals. Red points show counts of adult female mosquitoes (M) aggregated over individual villages for the first recorded rainy season in the data; for villages 4 and 5, a second rainy season denoted by S2 is also fitted to. Parameters for simulations were obtained from the median posteriors estimated by pMCMC fitting and 95% credible intervals from repeated samples of the joint posterior estimate

Table 1 Model parameters with priors, posterior estimates and 95% credible intervals for linear and fitted power density dependence, clumped and non-clumped egg-laying

Egg clumping and clutch size estimates

For both clumped and non-clumped egg-laying, when using a linear (Eq. 2) and fitted power density (Eq. 3) dependence, the model fitted to the data well. Again, however, it was not possible to fit an exponential density dependence (Eq. 4). DIC values show that fitting a power provided the best model, 1107.822, compared to 1122.032 for linear density. These values also show that non-clumped egg-laying provides a better fit to the data under both density dependence regimes, although this difference is marginal, particularly when fitting a power. Simulations of mosquito population dynamics using median posterior parameter estimates are shown in Figs. 2 and 3 for non-clumped and clumped egg-laying, respectively, of the fitted power density dependence (linear density dependence fits can be seen in Additional file 1).

Fig. 3
figure3

Model fits to Garki Project data for clumped egg-laying when fitting a power to density dependence with 95% credible intervals. Red points show counts of adult female mosquitoes (M) aggregated over individual villages for the first recorded rainy season in the data; for villages 4 and 5, a second rainy season denoted by S2 is also fitted to. Parameters for simulations were obtained from the median posteriors estimated by pMCMC fitting and 95% credible intervals from repeated samples of the joint posterior estimate

Clutch size \(n\) for fitted power (Eq. 3) and linear density dependence (Eq. 2) was estimated at 12.049 (95% CI 2.945–28.051) and 3.19 (95% CI 1.563–5.925), respectively. When considering non-clumped egg-laying, eggs laid per day \(\beta\) was estimated as 11.502 (95% CI 2.544–26.585) for fitted power and 1.305 (95% CI 1.00–2.372) for linear. Note that, as previously described, for the continuous (non-clumped) egg-laying model variant, when estimating \(\beta\), the gonotrophic cycle \(S\) is not explicitly considered (unlike in the clumped egg-laying equation); therefore, a comparative value to \(n\) (i.e. number of eggs laid in the same period) would be \(\frac{{ \beta \left( {e^{{S\mu_{{\text{m}}} }} - 1} \right)}}{{\mu_{{\text{m}}} }}\), giving 39.566 (95% CI 8.751–91.449) for the fitted power and 4.489 (95% CI 3.439–8.157) for linear density dependence.

High levels of correlation between \(\Omega\) and egg-laying also occurred during model fitting, where lower input of eggs into the system compensates for an increase in density-dependent mortality (Fig. 4). This is particularly pronounced under clumped egg-laying scenarios, showing that clutch size and the form of density dependence are intrinsically linked.

Fig. 4
figure4

Cross-correlation plots between \(n\) and \(\Omega ,\) showing the corresponding two dimensions of parameter space explored by the pMCMC algorithm; the hexagon colour and count value represent the number of accepted parameter proposals

Mosquito population growth rate

We found that simulated and raw analytical estimates of Rm were marginally lower for clumped egg-laying than continuous oviposition (Table 2), reflecting the finite lower bound on density-dependent mortality in the clumped model. However, for both clumped and non-clumped egg-laying, Rm was significantly lower when considering a linear density dependence over a fitted power.

Table 2 Rm estimates with 95% credible intervals for clumped and non-clumped egg-laying under linear and fitted power density dependence

Effect on gene drive

Both clumped egg-laying and fitting a power to the density dependence affected the required number and homing rate of HEGs for successful population elimination (Fig. 5). As expected, this is in line with the changes Deredec et al. [20] found when they modified values for Rm. Changes in density dependence had the greatest effect, whilst clumped egg-laying showed a relatively minor effect that cannot be concluded to be significant due to the overlap of credible intervals.

Fig. 5
figure5

Repeat of model by Deredec et al. [20] using parameters estimated by our analysis, Rm was derived both analytically and numerically for all density and egg-laying scenarios. The model estimates the number of HEGs needed in relation to their homing rate (a measure of efficacy) for successful elimination of a mosquito population. The bold centre line is the parameter estimates from the median of the posterior; the shaded bands represent the 95% credible intervals

Discussion

Here we have quantified the effects of differing density dependence regimes and mosquitoes laying their eggs in a clumped, non-homogeneous fashion. From this we can see how both egg-laying and density dependence have a clear impact on estimates of Rm and subsequent vector control efficacy.

Much recent development work on genetically modified mosquitoes has aimed at increasing homing and transmission rates of constructs [44,45,46,47]. In this context, our analysis suggests that the challenges to elimination posed by high values of Rm could vary substantially, especially if density dependence and egg-laying regimes are not universal. The consequences of this finding could have a significant impact on both the efficacy required from gene drive constructs and the release strategies necessary for successful deployment both positive and negative. Numerical and analytical estimates of Rm were significantly lower under a linear density dependence and clumped egg-laying scenario than under continuous egg-laying and a fitted power density dependence (Table 2). This level of reduction in Rm is substantial when considering predictions of the required homing rate and number of HEGs from our runs of the model by Deredec et al. [20].

Whilst we have looked only at HEG efficacy in this study, the variation in Rm will also be further compounded by wild-type resistance. For example, in work by Beaghton et al. [45] looking specifically at mutation rates in driving-Y interventions, Rm was key to the probability of mutations arising. This occurs because, with increasing Rm, the time to elimination is increased, and thus the probability that resistant mutations arise before population elimination is achieved. In addition, the probability of a resistant mutation becoming established also increases, as at a higher Rm the probability of stochastic loss decreases.

The fitting of a value for \(\Omega\) < 1 for both clumped and non-clumped egg-laying scenarios suggests a negative, non-linear density effect, i.e. as the larval population increases towards the environmental carrying capacity, the rise in severity of density effects decreases. It could be suggested from these results that whilst the overall death rate in the population is increasing, at the individual scale the mortality rate is decreased. This could be explained by there being a protective effect or efficiency gain in foraging/energy expenditure within larger groups of larvae, as often seen in other organisms [48, 49]. Additionally, the effects of density dependence may be mitigated to an extent by deceased or smaller larvae acting as an alternate food source.

The difference in Rm estimates between density dependence regimes is quite stark; however, it is difficult to identify the most probable functional form, particularly as DIC values are relatively close. In the wild, there is likely variation between sites and locations which can impact on larval development and adult fitness. These may be dependent on multiple factors, e.g. highly localised microclimates or differences in faunal communities. What this study perhaps more clearly shows is the importance of understanding this component of life history and the need for empirical data to help predict density dependence effects. This is something which was also highlighted in Khamis et al. [29], who similarly identified how density dependence could have a significant impact on vector control and cost-effectiveness. Additionally, our Rm estimates of non-clumped laying are lower than those from the comparable deterministic model by White et al. [28], likely a result of the more realistic, stochastic nature of our model and the fitting techniques used.

Estimated values for n and \(\beta\) in our study are relatively low, and lower than reported clutch sizes of Anopheles. However, the class E in our model is not a value for eggs, but early instar larvae, and thus will be lower than actual egg numbers when egg survival/hatching is considered. In addition, as the model does not include space or attempt to model individual breeding sites, the representation of density is not precise.

Generally, our model fitted to the data well, despite over- or underestimating some of the more extreme data points, particularly the higher peaks. However, it is important to consider how much overdispersion is inherent in such population counts. Further, as our model does not simulate individual breeding sites, there is likely variation between and within villages we do not capture. For example, we assume a single mosquito species and do not consider the interactions and competition with other species in larval pools. Whilst the environmental niches of mosquitoes both within and outside the Anopheles gambiae complex vary (some species such as An. arabiensis prefer larger or more permanent pools, whilst An. gambiae s.s. prefer smaller temporary sites [50]), there is still the possibility of interspecies competition for resources, and transient environmental conditions may lead to temporary niche overlap. As rainfall leads to higher carrying capacity, it also affects the type of habitat available, which may vary between villages/sites, and this could lead to a shift in species-specific density dependence effects. Lastly, in focusing solely on density dependence in larval mortality, our model does not account for density or other detrimental effects in the adult population, such as the size and availability of blood meals and infection status [51, 52], or the adult fitness costs resulting from density-dependent impeded larval development [53,54,55].

Conclusions

In this study, we have identified how general assumptions over egg-laying and the functional form of density dependence in mathematical models can have a significant impact on the estimates of life-history parameters. These in turn fundamentally affect the likely impact of vector control interventions such as gene drive. Our results suggest that these aspects of mosquito ecology should not be ignored in future studies, and understanding of local population dynamics are vital in our predictions of vector control efficacy. To build on this work, future data collection should consider detailed temporally regular, and preferably spatially stratified, population counts of adults and larvae in highly seasonal settings to better characterise these effects. Further, it is important we understand how these effects may change under differing regimes of mosquito and aquatic diversity, predator abundance and infection status, necessitating the need for comprehensive longitudinal studies conducted in realistic environments.

Availability of data and materials

Data is publicly available at http://garkiproject.nd.edu.

Code availability

Code is available at: https://github.com/aaronm70/mosModel.

References

  1. 1.

    WHO. World malaria report 2017. Geneva: World Health Organization; 2017.

    Google Scholar 

  2. 2.

    Burt A, Coulibaly M, Crisanti A, Diabate A, Kayondo JK. Gene drive to reduce malaria transmission in sub-Saharan Africa. J Responsible Innov. 2018;5(sup1):66–80.

    Article  Google Scholar 

  3. 3.

    Sinka ME, Bangs MJ, Manguin S, Coetzee M, Mbogo CM, Hemingway J, Patil AP, Temperley WH, Gething PW, Kabaria CW. The dominant Anopheles vectors of human malaria in Africa, Europe and the Middle East: occurrence data, distribution maps and bionomic précis. Parasites Vectors. 2010;3(1):117.

    PubMed  PubMed Central  Article  Google Scholar 

  4. 4.

    Guerra CA, Gikandi PW, Tatem AJ, Noor AM, Smith DL, Hay SI, Snow RW. The limits and intensity of Plasmodium falciparum transmission: implications for malaria control and elimination worldwide. PLoS Med. 2008;5(2):e38.

    PubMed  PubMed Central  Article  Google Scholar 

  5. 5.

    Wiebe A, Longbottom J, Gleave K, Shearer FM, Sinka ME, Massey NC, Cameron E, Bhatt S, Gething PW, Hemingway J. Geographical distributions of African malaria vector sibling species and evidence for insecticide resistance. Malar J. 2017;16(1):85.

    PubMed  PubMed Central  Article  Google Scholar 

  6. 6.

    Coetzee M. Distribution of the African malaria vectors of the Anopheles gambiae complex. Am J Trop Med Hyg. 2004;70(2):103–4.

    PubMed  Article  Google Scholar 

  7. 7.

    White GB. Anopheles bwambae sp. n., a malaria vector in the Semliki Valley, Uganda, and its relationships with other sibling species of the An. gambiae complex (Diptera: Culicidae). Syst Entomol. 1985;10(4):501–22.

    Article  Google Scholar 

  8. 8.

    Gillies MT, Coetzee M. A supplement to the Anophelinae of Africa south of the Sahara (Afrotropical region). Publ S Afr Inst Med Res. 1987;55:1–43.

    Google Scholar 

  9. 9.

    Coluzzi M. Heterogeneities of the malaria vectorial system in tropical Africa and their significance in malaria epidemiology and control. Bull World Health Organ. 1984;62(Suppl):107.

    PubMed  PubMed Central  Google Scholar 

  10. 10.

    Harbach RE. The classification of genus Anopheles (Diptera: Culicidae): a working hypothesis of phylogenetic relationships. Bull Entomol Res. 2004;94(6):537–53.

    CAS  PubMed  Article  Google Scholar 

  11. 11.

    Koenraadt CJM, Githeko AK, Takken W. The effects of rainfall and evapotranspiration on the temporal dynamics of Anopheles gambiae ss and Anopheles arabiensis in a Kenyan village. Acta Trop. 2004;90(2):141–53.

    CAS  PubMed  Article  Google Scholar 

  12. 12.

    Pascual M, Cazelles B, Bouma MJ, Chaves LF, Koelle K. Shifting patterns: malaria dynamics and rainfall variability in an African highland. Proc R Soc Lond B Biol Sci. 2008;275(1631):123–32.

    CAS  Google Scholar 

  13. 13.

    Kristan M, Abeku TA, Beard J, Okia M, Rapuoda B, Sang J, Cox J. Variations in entomological indices in relation to weather patterns and malaria incidence in East African highlands: implications for epidemic prevention and control. Malar J. 2008;7(1):231.

    PubMed  PubMed Central  Article  Google Scholar 

  14. 14.

    Floore TG. Mosquito larval control practices: past and present. J Am Mosq Control Assoc. 2006;22(3):527–33.

    CAS  PubMed  Article  Google Scholar 

  15. 15.

    Lengeler C. Insecticide-treated bed nets and curtains for preventing malaria. Cochrane Database Syst Rev. 2004. https://doi.org/10.1002/14651858.CD000363.pub2.

    Article  PubMed  Google Scholar 

  16. 16.

    Mabaso MLH, Sharp B, Lengeler C. Historical review of malarial control in southern African with emphasis on the use of indoor residual house-spraying. Trop Med Int Health. 2004;9(8):846–56.

    PubMed  Article  Google Scholar 

  17. 17.

    Kyrou K, Hammond AM, Galizi R, Kranjc N, Burt A, Beaghton AK, Nolan T, Crisanti A. A CRISPR–Cas9 gene drive targeting doublesex causes complete population suppression in caged Anopheles gambiae mosquitoes. Nat Biotechnol. 2018;36(11):1062–6.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  18. 18.

    Burt A. Site-specific selfish genes as tools for the control and genetic engineering of natural populations. Proc R Soc Lond B Biol Sci. 2003;270(1518):921–8.

    CAS  Article  Google Scholar 

  19. 19.

    Burt A, Koufopanou V. Homing endonuclease genes: the rise and fall and rise again of a selfish element. Curr Opin Genet Dev. 2004;14(6):609–15.

    CAS  PubMed  Article  Google Scholar 

  20. 20.

    Deredec A, Godfray HCJ, Burt A. Requirements for effective malaria control with homing endonuclease genes. Proc Natl Acad Sci. 2011;108(43):E874–80.

    CAS  PubMed  Article  Google Scholar 

  21. 21.

    Godfray HCJ, North A, Burt A. How driving endonuclease genes can be used to combat pests and disease vectors. BMC Biol. 2017;15(1):81.

    PubMed  PubMed Central  Article  Google Scholar 

  22. 22.

    North A, Burt A, Godfray HCJ. Modelling the spatial spread of a homing endonuclease gene in a mosquito population. J Appl Ecol. 2013;50(5):1216–25.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  23. 23.

    Unckless RL, Clark AG, Messer PW. Evolution of resistance against CRISPR/Cas9 gene drive. Genetics. 2017;205(2):827–41.

    PubMed  Article  Google Scholar 

  24. 24.

    Eckhoff PA, Wenger EA, Godfray HCJ, Burt A. Impact of mosquito gene drive on malaria elimination in a computational model with explicit spatial and temporal dynamics. Proc Natl Acad Sci. 2017;114(2):E255–64.

    CAS  PubMed  Article  Google Scholar 

  25. 25.

    Lambert B, North A, Burt A, Godfray HCJ. The use of driving endonuclease genes to suppress mosquito vectors of malaria in temporally variable environments. Malar J. 2018;17(1):154.

    PubMed  PubMed Central  Article  Google Scholar 

  26. 26.

    Hancock PA, White VL, Callahan AG, Godfray CHJ, Hoffmann AA, Ritchie SA. Density-dependent population dynamics in Aedes aegypti slow the spread of wMel Wolbachia. J Appl Ecol. 2016;53(3):785–93.

    Article  Google Scholar 

  27. 27.

    Muriu SM, Coulson T, Mbogo CM, Godfray HCJ. Larval density dependence in Anopheles gambiae ss, the major African vector of malaria. J Anim Ecol. 2013;82(1):166.

    PubMed  Article  Google Scholar 

  28. 28.

    White MT, Griffin JT, Churcher TS, Ferguson NM, Basáñez M-G, Ghani AC. Modelling the impact of vector control interventions on Anopheles gambiae population dynamics. Parasites Vectors. 2011;4(1):153.

    PubMed  PubMed Central  Article  Google Scholar 

  29. 29.

    Khamis D, El Mouden C, Kura K, Bonsall MB. Optimal control of malaria: combining vector interventions and drug therapies. Malar J. 2018;17(1):174.

    PubMed  PubMed Central  Article  Google Scholar 

  30. 30.

    Alphey N, Bonsall MB. Interplay of population genetics and dynamics in the genetic control of mosquitoes. J R Soc Interface. 2014;11(93):20131071.

    PubMed  PubMed Central  Article  Google Scholar 

  31. 31.

    Kitron UD, Webb DW, Novak RJ. Oviposition behavior of Aedes triseriatus (Diptera: Culicidae): prevalence, intensity, and aggregation of eggs in oviposition traps. J Med Entomol. 1989;26(5):462–7.

    CAS  PubMed  Article  Google Scholar 

  32. 32.

    Fay RW, Perry AS. Laboratory studies of oviposilional preferences of Aedes aegypti. Mosq News. 1965;25(3):276–81.

    Google Scholar 

  33. 33.

    Chadee DD, Corbet PS, Greenwood JJD. Egg-laying yellow fever mosquitoes avoid sites containing eggs laid by themselves or by conspecifics. Entomol Exp Appl. 1990;57(3):295–8.

    Article  Google Scholar 

  34. 34.

    Corbet PS, Chadee DD. An improved method for detecting substrate preferences shown by mosquitoes that exhibit ‘skip oviposition.’ Physiol Entomol. 1993;18(2):114–8.

    Article  Google Scholar 

  35. 35.

    Chadee DD, Corbet PS. Seasonal incidence and diel patterns of oviposition in the field of the mosquito, Aedes aegypti (L.) (Diptera: Culicidae) in Trinidad, West Indies: a preliminary study. Ann Trop Med Parasitol. 1987;81(2):151–61.

    CAS  PubMed  Article  Google Scholar 

  36. 36.

    Colton YM, Chadee DD, Severson DW. Natural skip oviposition of the mosquito Aedes aegypti indicated by codominant genetic markers. Med Vet Entomol. 2003;17(2):195–204.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  37. 37.

    Chen H, Fillinger U, Yan G. Oviposition behavior of female Anopheles gambiae in western Kenya inferred from microsatellite markers. Am J Trop Med Hyg. 2006;75(2):246–50.

    CAS  PubMed  Article  Google Scholar 

  38. 38.

    Bayoh MN, Lindsay SW. Effect of temperature on the development of the aquatic stages of Anopheles gambiae sensu stricto (Diptera: Culicidae). Bull Entomol Res. 2003;93(5):375–81.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  39. 39.

    Molineaux L, Gramiccia G, World Health Organization. The Garki project: research on the epidemiology and control of malaria in the Sudan savanna of West Africa. Geneva: World Health Organization; 1980.

    Google Scholar 

  40. 40.

    Baber I, Keita M, Sogoba N, Konate M, Doumbia S, Traore SF, Ribeiro JMC, Manoukis NC. Population size and migration of Anopheles gambiae in the Bancoumana Region of Mali and their significance for efficient vector control. PLoS ONE. 2010;5(4):e10270.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  41. 41.

    Charlwood JD, Smith T, Kihonda J, Heiz B, Billingsley PF, Takken W. Density independent feeding success of malaria vectors (Diptera: Culicidae) in Tanzania. Bull Entomol Res. 1995;85(1):29–35.

    Article  Google Scholar 

  42. 42.

    Andrieu C, Doucet A. Particle filtering for partially observed Gaussian state space models. J R Stat Soc Series B (Stat Methodol). 2002;64(4):827–36.

    Article  Google Scholar 

  43. 43.

    Doucet A, Johansen AM. A tutorial on particle filtering and smoothing: fifteen years later. Handb Nonlinear Filter. 2009;12(656–704):3.

    Google Scholar 

  44. 44.

    Champer J, Liu J, Oh SY, Reeves R, Luthra A, Oakes N, Clark AG, Messer PW. Reducing resistance allele formation in CRISPR gene drive. Proc Natl Acad Sci. 2018;115(21):5522–7.

    CAS  PubMed  Article  Google Scholar 

  45. 45.

    Beaghton A, Hammond A, Nolan T, Crisanti A, Godfray HCJ, Burt A. Requirements for driving antipathogen effector genes into populations of disease vectors by homing. Genetics. 2017;205(4):1587–96.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  46. 46.

    Marshall JM, Buchman A, Akbari OS. Overcoming evolved resistance to population-suppressing homing-based gene drives. Sci Rep. 2017;7(1):3776.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  47. 47.

    Hammond AM, Kyrou K, Bruttini M, North A, Galizi R, Karlsson X, Kranjc N, Carpi FM, D’Aurizio R, Crisanti A. The creation and selection of mutations resistant to a gene drive over multiple generations in the malaria mosquito. PLoS Genet. 2017;13(10):e1007039.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  48. 48.

    Heller R, Milinski M. Optimal foraging of sticklebacks on swarming prey. Anim Behav. 1979;27:1127–41.

    Article  Google Scholar 

  49. 49.

    Tollrian R, Dodson SI. Inducible defenses in Cladocera: constraints, costs, and multipredator environments. In: The ecology and evolution of inducible defenses. Princeton: Princeton University Press; 1999. p. 177–202.

    Google Scholar 

  50. 50.

    Fillinger U, Sombroek H, Majambere S, van Loon E, Takken W, Lindsay SW. Identifying the most productive breeding sites for malaria mosquitoes in The Gambia. Malar J. 2009;8(1):62.

    PubMed  PubMed Central  Article  Google Scholar 

  51. 51.

    Nayar J, Sauerman D Jr. The effects of nutrition on survival and fecundity in Florida mosquitoes: Part 2. Utilization of a blood meal for survival. J Med Entomol. 1975;12(1):99–103.

    CAS  PubMed  Article  Google Scholar 

  52. 52.

    Hurd H, Hogg J, Renshaw M. Interactions between bloodfeeding, fecundity and infection in mosquitoes. Parasitol Today. 1995;11(11):411–6.

    Article  Google Scholar 

  53. 53.

    Lyimo EO, Takken W, Koella JC. Effect of rearing temperature and larval density on larval survival, age at pupation and adult size of Anopheles gambiae. Entomol Exp Appl. 1992;63(3):265–71.

    Article  Google Scholar 

  54. 54.

    Reiskind MH, Lounibos LP. Effects of intraspecific larval competition on adult longevity in the mosquitoes Aedes aegypti and Aedes albopictus. Med Vet Entomol. 2009;23(1):62–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  55. 55.

    Scott DE. The effect of larval density on adult demographic traits in Ambystoma opacum. Ecology. 1994;75(5):1383–96.

    Article  Google Scholar 

Download references

Acknowledgements

We would like to thank the anonymous reviewers for their helpful and constructive feedback.

Funding

Funding for this work was provided by the Bill and Melinda Gates Foundation. The funders had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

Affiliations

Authors

Contributions

AM and NF devised and carried out the analysis; AM, NF and AG wrote the paper. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Aaron L. Morris.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

 Detailed model fitting methods and additional results.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Morris, A.L., Ghani, A. & Ferguson, N. Fine-scale estimation of key life-history parameters of malaria vectors: implications for next-generation vector control technologies. Parasites Vectors 14, 311 (2021). https://doi.org/10.1186/s13071-021-04789-0

Download citation

Keywords

  • Malaria
  • Mosquitos
  • Public health
  • Gene drive
  • Vector control
  • Population biology
  • Parameter estimation
  • Modelling