Skip to main content

Modelling the Wolbachia incompatible insect technique: strategies for effective mosquito population elimination

Abstract

Background

The Wolbachia incompatible insect technique (IIT) shows promise as a method for eliminating populations of invasive mosquitoes such as Aedes aegypti (Linnaeus) (Diptera: Culicidae) and reducing the incidence of vector-borne diseases such as dengue, chikungunya and Zika. Successful implementation of this biological control strategy relies on high-fidelity separation of male from female insects in mass production systems for inundative release into landscapes. Processes for sex-separating mosquitoes are typically error-prone and laborious, and IIT programmes run the risk of releasing Wolbachia-infected females and replacing wild mosquito populations.

Results

We introduce a simple Markov population process model for studying mosquito populations subjected to a Wolbachia-IIT programme which exhibit an unstable equilibrium threshold. The model is used to study, in silico, scenarios that are likely to yield a successful elimination result. Our results suggest that elimination is best achieved by releasing males at rates that adapt to the ever-decreasing wild population, thus reducing the risk of releasing Wolbachia-infected females while reducing costs.

Conclusions

While very high-fidelity sex separation is required to avoid establishment, release programmes tend to be robust to the release of a small number of Wolbachia-infected females. These findings will inform and enhance the next generation of Wolbachia-IIT population control strategies that are already showing great promise in field trials.

Background

Since the 1950s, the sterile insect technique (SIT) has been a popular biological control method for the management of pest insect populations [1]. The SIT takes advantage of the reproductive biology of a species to suppress or eliminate populations, applied through the inundative release of sterile individuals. Traditional approaches using radiation and chemosterilants have succeeded in the elimination of a variety of insect pests, including the New World screwworm in the Americas and Libya [2], the tsetse fly in Zanzibar [3] and the Mediterranean fruit fly in Mexico and Guatemala [4]. In a modern context, most SIT programmes now form a component of integrated pest management programmes, generating billions of dollars in savings for agricultural commerce annually, through the reduced environmental impacts of pesticide use [1]. The knowledge gained from these historical campaigns now provides baseline information for a renewed interest in area-wide SIT control programmes.

Rapid advances in fields such as molecular biology, genetics and computer science have seen a resurgence in SIT, particularly in attempting to prevent mosquito-borne diseases such as dengue, chikungunya, Zika and yellow fever [5, 6]. The incompatible insect technique (IIT) is closely related to SIT, but rather than releasing sterilised insects, this approach relies upon Wolbachia-infected male mosquitoes that are incapable of producing viable offspring after mating with a wild-type female [7]. Endosymbiotic bacteria, such as Wolbachia, may exhibit a biological mechanism known as cytoplasmic incompatibility (CI) [7]. When CI occurs, male mosquitoes infected with Wolbachia are not sterile per se, since they still produce offspring when mated with Wolbachia-infected females. However, in crosses with wild-type females not containing Wolbachia (or a compatible strain of Wolbachia), mating fails to produce viable offspring as eggs fail to hatch. While other methods of sterilisation, such as radiation and genetic modification, impose a large fitness burden [8] or suffer from complicated regulatory pathways [9], the use of a Wolbachia IIT approach offers a promising solution for controlling insect populations without these limitations [10].

Sterile mosquito strategies typically release only male mosquitoes, in part because it is the females that bite humans and spread disease. Introducing large numbers of females could conceivably aid disease spread [11] or contribute to nuisance biting. A challenge faced by the mass rearing of mosquitoes in IIT programmes is accurately separating males from females (sex separation). However, advances in the mechanical sorting of pupal and adult mosquitoes, as well as the availability of added safeguards such as irradiation, have seen renewed interest in controlling insect populations via the Wolbachia-IIT process [12, 13]. Whilst new methods can have a high degree of accuracy, occasional errors in sorting systems can lead to the release of Wolbachia-infected females. For IIT programmes that rely on Wolbachia-related CI effects, the release of female mosquitoes may have the unintended effect of population replacement rather than elimination. For a Wolbachia infection to become established within a population, it is generally thought that approximately 20% (depending upon the Wolbachia strain) of the population must become infected [14, 15]. Below this threshold, also known as the unstable equilibrium threshold (UET), the Wolbachia-infected genotype is likely to die out.

Historically, SIT and IIT control strategies have been underpinned by mathematical models, used to predict population dynamics of the target species. These models provide support tools for planning when, where and how many incompatible insects will be released into the wild population. The models of Knipling [16, 17] provided the original mathematical frameworks on which successful insect control campaigns have been based. Knipling’s models used discrete-time dynamics and a simple model of geometric population growth, where the resulting suppression depends on an overflooding ratio of sterile males to wild males [16, 17]. More recently, sophisticated SIT/IIT models have involved the modelling of genotypes and abundance via delay differential equations [18], spatio-temporal advection-diffusion-reaction models with multiple life stages [19] and complex agent-based simulations [20]. The majority of models used to study insect elimination processes via SIT/IIT are deterministic, ignoring the demographic randomness that is inherent to real populations. Demographic stochasticity can be a significant source of uncertainty in populations, particularly when populations are small and single births or deaths can have large effects on population dynamics [21]. Consequently, deterministic models do not quantify uncertain outcomes, such as the probability that a gene or Wolbachia symbiont will succeed or fail to enter a population. Models, such as that of Jansen et al. [22] or Magori et al. [20], are at a minimum stochastic but as we discuss below have certain attributes that affect their suitability for investigating important, general questions regarding the successful implementation of Wolbachia-IIT programmes.

Simulations from stochastic models yield ensembles of population trajectories which can be used to estimate probabilities of particular events (e.g. whether population elimination will occur). Where, in addition to stochasticity in population dynamics, there is uncertainty in the initial conditions or parameters that govern the biological system, these quantities can be sampled from probability distributions which summarise the scientific literature or from prior expert knowledge (called prior distributions). Combining a stochastic population model with prior distributions allows us to estimate the probabilities of events, whilst acknowledging uncertainties in parameters and demographic stochasticity.

Markov chains are ubiquitous for studying population dynamics, and such models have already been used to study the dynamics of Wolbachia infection within a population. For example, Jansen et al. [22] used a Moran process to estimate the probability of Wolbachia establishment after the release of a single infected female into a population of wild-type mosquitoes, but this process is not fit-for-purpose in SIT/IIT programmes, since it does not allow for a declining population.

Markov population processes (MPPs) [23] are a class of continuous-time Markov chain models suitable for exploring the effects of Wolbachia-IIT population suppression, and therefore the risk of Wolbachia establishment. The Markov property dictates that in an MPP, the times between events in the population follow an exponential distribution. Consequently, Markov population models typically assume that individual lifetimes are exponentially distributed. This is a reasonable assumption for mosquitoes in their natural habitat, where survival is frequently reported in terms of constant daily mortality [24, 25]. For phenomena where lifetimes are not exponentially distributed, a sequence of states can be chained together so that the total time spent in those collective states is the sum of exponentially distributed random variables.

In preparation for a large-scale Wolbachia IIT programme in North Queensland, Australia [26], it was necessary to explore hypothetical population replacement scenarios where Wolbachia-infected females were released unintentionally. Here, we use an MPP model to address key questions related to the effective implementation of a Wolbachia-IIT population suppression programme, namely:

  1. (i)

    What degree of fidelity in sex separation is required to ensure that a Wolbachia-IIT population suppression programme does not result in Wolbachia establishment and the replacement of the wild-type (non-infected) population?

  2. (ii)

    Does the release of a small number of females during a Wolbachia-IIT population suppression programme imply a high probability of Wolbachia establishment?

  3. (iii)

    What types of release strategies are most effective for achieving population elimination without Wolbachia establishment?

  4. (iv)

    How does the establishment probability vary with the initial density of the Wolbachia infection in the population?

We recognise the important trade-off between increasing model complexity and a reduced ability to defensibly parameterize the model using what is known from the scientific literature so that it can be applied more generally. To overcome this, we present a novel MPP model that relies upon a relatively small set of parameters and is specifically designed to study Wolbachia-IIT population dynamics at the suburban block level. Importantly, our model is structured so that the parameters can be defensibly estimated from previous field studies (see Additional File A: Table A1.1) [14, 15, 24, 27,28,29,30,31,32].

Methods

Simulation of Wolbachia-IIT programmes

Our studies focus on simulating the repeated release of Aedes aegypti, infected with wAlbB Wolbachia, into a hypothetical population that is intended to represent a typical suburban block for Innisfail, northern Australia (− 17.5226° S, 146.0285° E). Our intention is to obtain results that resemble a well-mixed biological population at the suburban block scale. Results obtained can be extrapolated to larger landscapes assuming independence between blocks. To allow for some heterogeneity between suburban block scale populations, our simulations of Wolbachia-IIT programmes use adult equilibrium populations (prior to wAlbB releases) in the range 100 to 300 wild-type individuals/block. If we consider that a typical suburban block in Innisfail has about 20 houses, this represents between 5 and 15 adult Ae. aegypti mosquitoes per house, which is typical for North Queensland [25]. For each simulated IIT programme, we modelled 1 year of data with wAlbB releases occurring on Monday, Wednesday and Friday. After a year of wAlbB releases, we continued to model the population for a further year without releases, to monitor for any additional changes in population structure. After this 2-year period, we classify each simulated trajectory as having an IIT endpoint of (i) eliminated, (ii) established, (iii) indeterminate Wolbachia negative or (iv) indeterminate Wolbachia positive (see Table 1 for definitions). We assess the effectiveness of a Wolbachia-IIT programme by calculating the proportions of simulations that reach these four IIT endpoints, with the elimination endpoint being the obvious goal.

Table 1 Glossary of modelling terminology

We simulated Wolbachia-IIT programmes using three release strategies (defined in Table 1): a constant release strategy, an adaptive release strategy and a crude adaptive release strategy. For each of these release strategies, we also considered two overflooding ratios of 5:1 and 15:1 (Table 1). These overflooding ratios scale the numbers of wAlbB individuals released under each of these three strategies. Release strategies and overflooding ratios in combination with four levels of se separation fidelity are determined by the female contamination probability (FCP; Table 1). We consider FCPs of 10−4, 10−5, 10−6 and 10−7, which span ranges that are theoretically achievable using next-generation mechanical sex separation approaches. We used 1000 simulations per scenario to assess the probability of wAlbB establishment and successful elimination. The additional simulations at the lowest FCP are to ensure that at least some population trajectories result in the release of wAlbB females. Each set of simulations was used to study the probabilities of (i) unintended wAlbB establishment and (ii) wild-type elimination. Large numbers of simulations (15,000) per scenario were also used in conjunction with importance sampling (see Robert and Casella [28]) to efficiently estimate establishment and elimination probabilities and their standard errors (see Additional File A: Tables A2.1 and A2.2).

In our simulations, we do not assume that parameters are known precisely. Instead, we use a uniform probability distribution over plausible ranges for each of the parameters in our model (discussed in the subsequent section). For each simulation, a set of parameters is sampled from within these plausible ranges and checked to ensure that certain mandatory constraints are met (see Additional File B: B2.1) [23, 33,34,35,36]. In doing so, we ensure that our modelling results acknowledge parametric uncertainty and consider a plausible degree of variation that might be exhibited in a natural population. As such, simulations that result in a high probability of elimination can be considered robust to a range of parameterisations and population outcomes.

The Wolbachia-IIT MPP model

The states through which individuals in our model can progress are depicted in Fig. 1. In total, there are 2k + 8 different classes of individuals in the model, and our MPP tracks the number of individuals in each class. For the wild-type and the wAlbB mosquito populations, there are k states which represent “future adults” (see Additional File B: B2.1-B2.3) denoted IX, 1, IX, 2, …, IX, k, where X {wild, wAlb}; adult males denoted MX; unmated females denoted FX; and mated females denoted FX × Y, where Y {wild, wAlb}. The precise mathematical details of the MPP are provided in Additional File B: B1-B4, but we outline the general mechanics below. The model (as an installable R package) is publicly available at https://github.com/dpagendam/debugIIT.

Fig. 1.
figure1

Model state transitions of wild-type and wAlbB-infected individuals from birth to death

The modelling approach is designed to be parsimonious and is premised on the idea that, for a population to persist, we expect individuals to produce new adults at a rate that balances deaths in the population. To accommodate this, our model includes an abstract pool of “future adults” (see Table 1), that is, a cohort of individuals that are immature, but will survive to become new adults in the population. Mated females in the population give birth to new future adults, and after a gamma-distributed, random period, each future adult emerges as a new adult that can mate and potentially give rise to new future adults. The rate at which mated females produce future adults is modelled to decrease as the pool of future adults increases in size and emulates density-dependent larval mortality without resorting to modelling the entire larval pool.

In Fig. 1, we see that the birth of a new individual introduces it into the first class of the “future adults” (IX, 1), and individuals in each of the k future adult classes transition out of each state at rate γ. For details regarding the choice of k states and its use in creating biologically realistic models via what is known as the “Linear Chain Trick”, we refer the reader to Additional File B:B1. When an individual transitions out of state IX, k, it is allocated to the pool of males or unmated females with equal probability. The rate at which females are mated is assumed to be independent of the density of males in the population, and provided that there are males in the population, each female transitions to a mated state at rate η. We assume that females only mate once and that her mate is either wild-type with probability \( \frac{M_{\mathrm{wild}}}{M_{\mathrm{wild}}+{c}_{w\mathrm{Alb}}{M}_{w\mathrm{Alb}}} \) or wAlbB with probability \( \frac{c_{w\mathrm{Alb}}{M}_{w\mathrm{Alb}}}{M_{\mathrm{wild}}+{c}_{w\mathrm{Alb}}{M}_{w\mathrm{Alb}}} \). In other words, the rate at which females are mated by males is not density-dependent, but is mate-type. The parameter cwAlb is the mating competitiveness coefficient or Fried’s Index of wAlbB-infected individuals relative to wild-type individuals [32]. We assume that CI between wAlbB males and wild-type females is 100% effective, so that only Fwild × wild females can give birth to wild-type future adults. Mated females of type FwAlb × wild and FwAlb × wAlb can give birth to wAlbB future adults. Each mated female adds new “future adults” to state IX, 1 at rate \( \frac{\lambda \left({I}_{\mathrm{max}}-{I}_{\mathrm{total}}\right)}{I_{\mathrm{max}}} \), where \( {I}_{\mathrm{total}}={\sum}_{i=1}^k\left({I}_{i,\mathrm{wild}}+{I}_{i,w\mathrm{Alb}}\right) \) and Imax is a parameter that acts as a population ceiling on the number of future adults allowed within the population (see Additional File B: B4.1-B4.2). The parameter λ is the per capita rate at which mated females produce future adults in an empty niche (i.e. the intrinsic rate of population growth). This growth rate is modulated by a density-dependent term \( \frac{\left({I}_{\mathrm{max}}-{I}_{\mathrm{total}}\right)}{I_{\mathrm{max}}} \) which can be thought of conceptually as a hypothetical proportion of niches for future adults which are currently unoccupied. Ultimately, this results in future adults being produced at a higher rate (respectively lower rate) when the density of future adults is low (respectively high). In essence, this provides a simple representation of higher larval mortality at higher larval population densities [29, 37]. Additional File B: B2.2 outlines how the abstract parameter, Imax, can be identified from a small number of easily estimated quantities presented in Additional File A: Table A1.1.

Our decision to model the female mating rate as a density-independent process stems from the assumption that adult males and females have little difficulty finding each other at the level of a suburban block and that the mating rate is not entirely governed by random interactions and may involve behavioural bottlenecks. Therefore, we argue that the rate of female mate selection is not simply increased monotonically by the presence of more males. There is also evidence to suggest that adult males and females can find each other by means of sensory cues [38], so that the female mating rate does not necessarily greatly diminish at low adult male densities. In modelling the probability of wAlbB establishment, we do introduce a conservative assumption that any wAlbB females accidentally released will have been pre-mated with a wAlbB male. This reflects an assumption of prolonged exposure of each wAlbB female to many wAlbB males in isolation.

We note that the mathematical model presented in Additional File B: B1-B4 is more general than that presented here and accommodates populations with multiple strains of Wolbachia infection. Additional File B: sections B1 and B2 also describe how the model relies on a set of primary and secondary parameters. Primary parameters are quantities that the user must specify to run the model, whilst secondary parameters are not defined by the user but are uniquely defined by the values of the primary parameters and the equilibrium dynamics of the model. Parameter ranges used as uniform prior probability distributions for simulations are provided in Additional File A: Table A1.1.

Most of the parameters listed in Additional File A: Table A1.1 can be estimated from field or laboratory data, and we provide references that were used for these purposes. Some of these parameters require some additional detail, namely k and γ, which together dictate that the mean time between a female laying an egg containing a potential future adult and the emergence of that new adult can range between 10 and 40 days (consistent with Hancock et al. [29]) with the distribution for the time spent in immature form following a gamma distribution. In addition, pmated was estimated in a series of (unpublished) experiments conducted over five different weeks using the rhodamine B marking approach of Johnson et al. [39], where wild-type females were captured and the spermatheca dissected to ascertain whether mating had taken place.

Unstable equilibrium threshold and the probability of wAlbB establishment

To ensure that our model was biologically defensible, we checked for the existence of a UET of the type discussed by Hoffmann et al. [40] and demonstrated in laboratory populations by Axford et al. [15]. These studies suggest that there should be a high probability of wAlbB becoming established when the number of wAlbB-infected adults introduced exceeds a wAlbB frequency of approximately 20%. There have been a number of theoretical investigations employing deterministic models to study the UET [41,42,43]. To examine the existence of a UET in our stochastic system, we ran simulations of our model using a wild-type population with an equilibrium of 400 adult individuals (to maintain similarity to the cage experiments of Axford et al. [15]). We generated simulations of our population where wAlbB adults in a 1:1 sex ratio were introduced to this population, representing between 5 and 40% (in intervals of 5%) of the total adult population. For each proportion of introduced wAlbB adults, we simulated 100 population trajectories over a 5-year period. Simulations were conducted under two further scenarios: the first where none of the released wAlbB females had been mated prior to release, and the second where all wAlbB females were assumed to have been mated by wAlbB males prior to release. The establishment probabilities for each set of 100 simulations were estimated as the proportions of trajectories where the wild-type populations were completely replaced by the wAlbB-infected strain. Confidence intervals around the estimated establishment probabilities were derived using profile likelihood (i.e. by inverting the likelihood ratio test statistic).

Results

We ran a total of 114,000 individual simulations each taking an average of 17 min to run on a single CPU core, with 90,000 of those simulations for importance sampling estimates of establishment and elimination probabilities. This amounted to a total of 1345 days of HPC compute time with runs parallelised to perform up to 5000 concurrent simulations at maximum capacity. Each simulation was allocated 8 GB of RAM and was run on an Intel Xeon E5-2670 V3 processor running at 2.6 GHz.

In addressing the four questions proposed by this study, we focused on a small set of plots and tables that demonstrate the following key results:

  1. (i)

    What degree of fidelity in sex separation is required?

To ensure a very low probability (< 0.01) of wAlbB establishment on a single suburban block, an FCP of 10−6 or less is advisable, but this probability is also a function of the number of mosquitoes released and is affected by the release strategy and overflooding ratio. The proportions of simulations whose IIT endpoints were classed as either successful elimination, wAlbB establishment, indeterminate wAlbB negative or indeterminate wAlbB positive are shown in Fig. 2. These proportions are plotted for different sex separation fidelities, release strategies and target populations. We observed a clear decline in the incidence of establishment events with lower FCPs.

  1. (ii)

    Does the release of a small number of females during a Wolbachia-IIT population suppression programme imply a high probability of Wolbachia establishment?

Fig. 2.
figure2

Proportions of simulations that either resulted in wAlbB establishment (red), successful wild-type elimination (green), indeterminate wAlbB negative (blue) and indeterminate wAlbB positive (yellow) for scenarios based on different FCPs and release strategies at 5:1 (top row) and 15:1 (bottom row) overflooding. Results suggest that a lower FCP leads to higher elimination rates. At the lower FCP, the best results appear to be for the constant and crude adaptive release strategies

To address the question of whether the accidental release of a small number of wAlbB females is likely to result in a high probability of establishment, we classified simulations by their IIT endpoints. Figure 3 dissects some of these simulated trajectories (15:1 overflooding, with FCP of 10−6 and 10−7) in this way. The plots serve to provide useful statistics (minimum, median and maximum) for the numbers of males and females released for each of the four Wolbachia-IIT endpoints. Additionally, Fig. 2 and Additional File A: Tables A3.1 and A3.2 provide importance sampling estimates of the wAlbB establishment probabilities. Notably, this statistical approach provides statistically efficient estimates of the establishment probability, even when simulated establishment events are rarely encountered at the FCP of 10−7.

  1. (iii)

    What types of release strategies are most effective?

Fig. 3.
figure3

Statistics for the number of males and females released using a 15:1 overflooding and FCP of 10−7 (top row) and 10−4 (bottom row). Rectangles span the minima and maxima for each set of simulations; horizontal red and blue lines show the medians for each set of simulations, and numbers in brackets show the number of simulations (out of 1000) that resulted in each IIT endpoint

To ensure a very low probability (< 0.01) of any block-level wAlbB establishment over an area spanning as many as 100 suburban blocks, it is advisable to adopt the crude adaptive release strategy with an FCP of 10−7 or less using a 5:1 overflooding ratio (Additional File A: Fig. A3.1 and Table A2.1). Furthermore, the adaptive release strategy demonstrated a high likelihood of achieving an endpoint of indeterminate wAlbB negative across all FCPs and overflooding ratios studied. Of the three release strategies examined, the best IIT outcomes (i.e. high probability of elimination and low probability of establishment) were achieved using a crude adaptive release strategy with a 5:1 overflooding ratio at lower FCPs (Fig. 2 and Additional File A: Tables A2.1 and A2.2).

  1. (iv)

    How does the establishment probability vary with the initial density of the Wolbachia infection in the population?

The stochastic population model developed and used in this study demonstrates a similar UET to that observed previously in laboratory cage experiments (Fig. 4 and Additional File A: Fig. A3.2). We demonstrated that our model (and its parameterisation) yields simulated population dynamics consistent with observations of unstable equilibria in biological systems. In a population having a stable wild-type equilibrium of 400 individuals (i.e. Kwild = 400), we simulated trajectories with different initial proportions of wAlbB and wild-type individuals. We display the first ten simulated population trajectories summarised as the proportion of adults infected with wAlbB for both frequency scenarios where no pre-release mating has occurred (Fig. 4). In those simulations where there was an initial 25% wAlbB frequency, 62 out of 100 simulations ended in wAlbB establishment, whereas at the smaller initial frequency of 10% wAlbB, only 25 out of 100 trajectories established. The point at which the establishment probability exceeded a 50% probability occurred at between 15 and 20% wAlbB for no pre-release mating and between 10 and 15% where pre-release mating was present. To quantify the differences in wAlbB establishment probabilities between mated and unmated introduced females, Additional File A: Fig. A3.2 displays the estimated establishment probabilities for each wAlbB frequency level using 100 simulated populations at each of these levels. The datasets supporting the conclusions of this article are available in the CSIRO Data Access Portal repository [44].

Fig. 4.
figure4

A random sample of ten simulated trajectories of the cage experiment scenario (no wAlbB pre-release mating) at wAlbB initial proportions of 10% (left) and 25% (right)

Discussion

Wolbachia IIT shows promise as an alternative for radiation and genetic modification (e.g. RIDL, gene drive) methods for insect population suppression. This is primarily due to the lower fitness costs of certain Wolbachia strains [15] and a well-developed implementation pathway, with lower barriers to regulatory approval and high community acceptance [45]. However, the release of Wolbachia-infected females has the potential to jeopardise suppression activities by replacing the local mosquito population and rendering the intervention unsuccessful. Taking a “do no harm” philosophy essential before implementing a novel field intervention, it was necessary to explore these concerns within an appropriate modelling framework.

The development of the MPP presented here has allowed a deeper understanding of the uncertainties associated with female contamination frequencies in Wolbachia-IIT programmes and highlights the benefits of adopting stochastic rather than deterministic modelling frameworks. Importantly, our model highlights the attributes of Wolbachia-IIT programmes that yield high probabilities of eliminating a population without Wolbachia establishment: (i) very high-fidelity sex separation and (ii) overflooding ratios and release strategies that employ only moderate numbers of Wolbachia-infected males. The former point highlights a need for technologies that can effectively discriminate between sexes and sort mosquitoes without reducing their fitness. The latter highlights a need to understand the long-term outlook when conducting such IIT programmes. It follows that releasing an excess of Wolbachia-infected male mosquitoes will undermine a programme by increasing the potential number of females accidentally released, while increasing the demand on rearing facilities. This is particularly the case as the first area-wide SIT/IIT programmes have typically employed a constant or increasing release rate over time [46, 47]. If the rate at which females are mated is density-independent and only the mate type is density-dependent (as we have assumed in our model), then very high overflooding ratios hold little strategic value over moderate ones. Figure 3 shows that releasing very large numbers of mosquitoes (e.g. 15:1 overflooding under the constant release strategy compared to the crude adaptive release strategy) tends to result in a greater likelihood of wAlbB establishment.

Of the three release strategies and two overflooding ratios simulated, the best performing approach was the crude adaptive release strategy in conjunction with a 5:1 overflooding ratio. The adaptive release strategy appeared to be ineffective at both the 5:1 and 15:1 overflooding ratios, as the number of males released over 1 year resulted in many outcomes of “indeterminate Wolbachia negative” (population suppression without elimination or Wolbachia establishment). It is possible that this approach may achieve satisfactory elimination when releases are scheduled over a longer period of time; however, given the observed success of the crude adaptive release strategy at 5:1 overflooding, we did not pursue this further. Furthermore, implementing the adaptive release strategy seems difficult in practice, requiring accurate, real-time surveillance and population size estimates to tailor release numbers appropriately. The constant release strategy also showed encouraging results at the lowest FCP and 5:1 overflooding ratio. However, we estimate the establishment probability for the crude adaptive release strategy to be an order of magnitude lower than for constant overflooding and is therefore preferable for mitigating the risk of Wolbachia establishment (particularly when applied over larger urban areas). An aspect that was not investigated in our modelling is the probability of establishment when a female is released towards the end of an IIT programme, when compared to earlier stages. This reflects a reduced barrier to entering a niche when larval densities are lower because of reduced density-dependent mortality.

The simulations in this study were intended to represent suburban blocks in North Queensland, Australia, which typically consist of ~ 20 houses. Since a population suppression programme would typically involve treating larger urban areas, the probability of one or more blocks resulting in wAlbB establishment will increase over regions consisting of multiple blocks. Under a simplifying assumption that all blocks in a treated region are independent, non-interacting subpopulations, the probability of wAlbB establishing on one or more blocks is calculated as π = 1 − (1 − p)n, where p is the block level establishment probability and n is the number of blocks treated with wAlbB. For a hypothetical region consisting of 100 blocks, a block-level wAlbB establishment rate of less than 10−4 would be required in order to ensure that the overall probability of one or more blocks having an establishment event was less than 0.01. Based on this, the most suitable release strategy investigated would be a crude adaptive release strategy with an FCP of 10−7 and an overflooding ratio of 5:1 or 15:1, which have block-level wAlbB establishment rates of 1.546 × 10−5 and 4.159 × 10−5, respectively.

We believe that our estimates of the probability of establishment risk are conservative, since we assumed that all wAlbB females were mated and yielded viable offspring. Through these simulations, we hope to debunk the legitimate concern amongst biologists that IIT programmes may not be robust to the release of a small number of Wolbachia-infected females and the subsequent risk of establishment. The probabilistic evidence that we have presented here suggests that Wolbachia establishment depends on many factors, including the population size at the time of accidental female release, the release strategy and the overflooding ratio employed. In many of our simulated IIT programmes, where females were accidentally released, Wolbachia failed to establish, primarily due to the presence of the UET coupled with demographic stochasticity. To demonstrate this point, Fig. 3 shows that as many as 23 females were released in one of the trajectories at the 10−4 FCP, and yet, elimination was achieved, and the population was successfully eradicated. Major factors affecting establishment in the population are whether a female survives long enough post-release to produce new individuals and the population density of larvae that reduces the probability of an individual female replacing herself. As such, it is important to acknowledge that if a small number of wAlbB females were detected in the population, then the chance of establishment would be further decreased if releases ceased for a period.

To increase confidence in the population dynamics exhibited in our MPP model, we endeavoured to replicate a scenario similar to the cage experiments of Axford et al. [15]. The key findings from our simulations were as follows: (i) for pre-mated wAlbB individuals, the majority of simulations resulted in establishment where the initial wAlbB frequency was above 20%, and (ii) for wAlbB individuals without pre-mating, the majority of simulations established above an initial 15% wAlbB frequency. These results show a general agreement with the Axford et al. [15] laboratory cage experiments where a high probability (five out of five laboratory populations) of wAlbB establishment occurred when the initial frequency of wAlbB was 22%. Likewise, both this study and Axford et al. [15] show a much lower probability of establishment (one out of five populations in the latter study) when the initial frequency of wAlbB in the population was 10%.

There are clear benefits to using stochastic models over deterministic models in population biology for addressing questions related to probability and risk. For estimating probabilities of events such as elimination or establishment, stochastic models provide a natural approach and capture important demographic stochasticity (especially when particular subpopulations of interest become small) that can drastically change the outcome for a population. However, in addressing our questions about wAlbB establishment risk, we required a model that was fit for purpose and easily parameterised. Existing stochastic models were deemed unsuitable. For example, the model of Jansen et al. [22] assumes that the population size remains constant over time (which is violated during IIT programmes) and that populations evolve in discrete time through non-overlapping generations, and Magori et al. [20] require the identification of many spatially explicit parameters (e.g. numbers and locations of breeding containers). The simple MPP model introduced in this work relied on a simple set of parameters that could be determined from the literature and past field data and was general enough to draw conclusions about a typical suburban block in Innisfail, Australia. We believe that this study demonstrates the usefulness of such MPP models for making important in silico observations about what factors might be important in designing an effective SIT/IIT strategy. Furthermore, our model accounts for lags between a female giving birth to a new individual and that individual’s eventual emergence as an adult, avoiding the complicated processes of tracking every individual in each life stage of a mosquito population or conditioning our results on detailed spatial information (e.g. location and number of containers) used in other mosquito modelling systems such as CimSim [48, 49] and Skeeter Buster [20]. Importantly, we avoid explicitly modelling all life stages because (i) the life history parameters of juveniles appear to be less well understood in the field and (ii) this would drastically increase the numbers of individuals that would need to be stochastically modelled in the population resulting in significant computational overheads. Furthermore, modelling every individual during each life cycle stage in a population may be of little real benefit since the vast majority of juveniles do not reach adulthood due to density-dependent larval mortality [29, 37]. We believe that the parsimonious approach adopted in our MPP model simplifies these processes and can be conceptually and computationally advantageous when running simulations over much larger spatial areas.

Conclusion

From the perspective of designing an effective IIT programme, our study suggests that to achieve elimination while avoiding establishment at large spatial scales (i.e. many suburban blocks), it is advisable to have a high-fidelity sex separation method with an FCP of 10−7 or less. Our results favour the crude adaptive release strategy at a 5:1 overflooding ratio to achieve both a low establishment probability and high elimination probability at the 10−7 FCP level. Our simulations indicate that releasing a relatively small number of wAlbB females during a programme does not automatically render it a failure.

Availability of data and materials

The datasets supporting the conclusions of this article are available in the CSIRO Data Access Portal repository [https://doi.org/10.25919/5f473485911aa]. The datasets and code supporting the conclusions of this article are included within the article [44] and its additional files (Additional File A and Additional File B).

Project name: Verily: Mosquito population Control

Project home page: https://research.csiro.au/mozzieproject/

Archived version: https://doi.org/10.25919/5f473485911aa

Operating system(s): Installable R package

Programming language: R

Licence: All Rights (including copyright) CSIRO 2020

Any restrictions to use by non-academics: The metadata and files are available to the public.

The population model code is publicly available as an installable R package at https://github.com/dpagendam/debugIIT.

References

  1. 1.

    Dyck VA, Hendrichs J, Robinson AS. Sterile insect technique: principles and practice in area-wide integrated pest management. Dordrecht: Springer; 2006.

  2. 2.

    Vargas-Terán M, Hofmann HC, Tweddle NE. Impact of screwworm eradication programmes using the sterile insect technique. In: Dyck VA, Hendrichs J, Robinson AS, editors. Sterile insect technique: principles and practice in area-wide integrated pest management. Dordrecht: Springer Netherlands; 2005. p. 629–50.

    Google Scholar 

  3. 3.

    Vreysen MJ, Saleh KM, Ali MY, Abdulla AM, Zhu ZR, Juma KG, et al. Glossina austeni (Diptera: Glossinidae) eradicated on the island of Unguja, Zanzibar, using the sterile insect technique. J Econ Entomol. 2000;93(1):123–35.

    CAS  Article  PubMed  Google Scholar 

  4. 4.

    Hendrichs J, Ortiz G, Liedo P, Schwarz A, editors. Six years of successful medfly program in Mexico and Guatemala. Fruit flies of economic importance. Proc. CEC/IOBC Int. Symp., November 1982; 1983.

  5. 5.

    Gilles JRL, Schetelig MF, Scolari F, Marec F, Capurro ML, Franz G, et al. Towards mosquito sterile insect technique programmes: exploring genetic, molecular, mechanical and behavioural methods of sex separation in mosquitoes. Acta Trop. 2014;132:S178–S87.

    Article  PubMed  Google Scholar 

  6. 6.

    Kittayapong P, Kaeothaisong N-o, Ninphanomchai S, Limohpasmanee W. Combined sterile insect technique and incompatible insect technique: sex separation and quality of sterile Aedes aegypti male mosquitoes released in a pilot population suppression trial in Thailand. Parasit Vectors. 2018;11(2):657.

    Article  CAS  PubMed  Google Scholar 

  7. 7.

    Laven H. Eradication of Culex pipiens fatigans through cytoplasmic incompatibility. Nature. 1967;216(5113):383–4.

    CAS  Article  PubMed  Google Scholar 

  8. 8.

    Helinski MEH, Parker AG, Knols BGJ. Radiation biology of mosquitoes. Malar J. 2009;8(2):S6.

    Article  PubMed  Google Scholar 

  9. 9.

    James S, Collins FH, Welkhoff PA, Emerson C, Godfray HCJ, Gottlieb M, et al. Pathway to deployment of gene drive mosquitoes as a potential biocontrol tool for elimination of malaria in Sub-Saharan Africa: recommendations of a Scientific Working Group. Am J Trop Med Hyg. 2018;98(6_Suppl):1–49.

    Article  PubMed  Google Scholar 

  10. 10.

    De Barro PJ, Murphy B, Jansen CC, Murray J. The proposed release of the yellow fever mosquito, Aedes aegypti containing a naturally occurring strain of Wolbachia pipientis, a question of regulatory responsibility. J Consum Prot Food S. 2011;6(1):33–40.

    Article  Google Scholar 

  11. 11.

    Alphey L, Benedict M, Bellini R, Clark GG, Dame DA, Service MW, et al. Sterile-insect methods for control of mosquito-borne diseases: an analysis. Vector Borne Zoonotic Dis. 2010;10(3):295–311.

    Article  PubMed  Google Scholar 

  12. 12.

    Atyame CM, Pasteur N, Dumas E, Tortosa P, Tantely ML, Pocquet N, et al. Cytoplasmic incompatibility as a means of controlling Culex pipiens quinquefasciatus mosquito in the islands of the south-western Indian Ocean. PLoS Negl Trop Dis. 2011;5(12):e1440.

    Article  PubMed  Google Scholar 

  13. 13.

    Zheng X, Zhang D, Li Y, Yang C, Wu Y, Liang X, et al. Incompatible and sterile insect techniques combined eliminate mosquitoes. Nature. 2019;572(7767):56–61.

    CAS  Article  PubMed  Google Scholar 

  14. 14.

    Xi Z, Khoo CCH, Dobson SL. Wolbachia establishment and invasion in an Aedes aegypti laboratory population. Science. 2005;310(5746):326–8.

    CAS  Article  PubMed  Google Scholar 

  15. 15.

    Axford JK, Ross PA, Yeap HL, Callahan AG, Hoffmann AA. Fitness of wAlbB Wolbachia infection in Aedes aegypti: parameter estimates in an outcrossed background and potential for population invasion. Am J Trop Med Hyg. 2016;94(3):507–16.

    Article  PubMed  Google Scholar 

  16. 16.

    Knipling EF. Possibilities of insect control or eradication through the use of sexually sterile males. J Econ Entomol. 1955;48(4):459–62.

    Article  Google Scholar 

  17. 17.

    Knipling EF. Sterile-male method of population control. Successful with some insects, the method may also be effective when applied to other noxious animals 1959;130(3380):902–904.

  18. 18.

    Hancock PA, Sinkins SP, Godfray HCJ. Population dynamic models of the spread of Wolbachia. Am Nat. 2011;177(3):323–33.

    Article  Google Scholar 

  19. 19.

    Dufourd C, Dumont Y. Impact of environmental factors on mosquito dispersal in the prospect of sterile insect technique control. Comput Math Appl. 2013;66(9):1695–715.

    Article  Google Scholar 

  20. 20.

    Magori K, Legros M, Puente ME, Focks DA, Scott TW, Lloyd AL, et al. Skeeter Buster: a stochastic, spatially explicit modeling tool for studying Aedes aegypti population replacement and population suppression strategies. PLoS Negl Trop Dis. 2009;3(9):e508.

    Article  PubMed  Google Scholar 

  21. 21.

    Hancock PA, Ritchie SA, Koenraadt CJM, Scott TW, Hoffmann AA, Godfray HCJ. Predicting the spatial dynamics of Wolbachia infections in Aedes aegypti arbovirus vector populations in heterogeneous landscapes. J Appl Ecol. 2019;56(7):1674–86.

    Article  Google Scholar 

  22. 22.

    Jansen VAA, Turelli M, Godfray HCJ. Stochastic spread of Wolbachia. Proc Biol Sci. 2008;275(1652):2769–76.

    PubMed  PubMed Central  Google Scholar 

  23. 23.

    Kingman JFC. Markov population processes. J Appl Probab. 1969;6(1):1–18.

    Article  Google Scholar 

  24. 24.

    Muir LE, Kay BH. Aedes aegypti survival and dispersal estimated by mark-release-recapture in northern Australia. Am J Trop Med Hygiene. 1998;58(3):277–82.

    CAS  Article  Google Scholar 

  25. 25.

    Ritchie SA, Montgomery BL, Hoffmann AA. Novel estimates of Aedes aegypti (Diptera: Culicidae) population size and adult survival based on Wolbachia releases. J Med Entomol. 2013;50(3):624–31.

    Article  PubMed  Google Scholar 

  26. 26.

    Staunton KM, Yeeles P, Townsend M, Nowrouzi S, Paton CJ, Trewin B, et al. Trap location and premises condition influences on Aedes aegypti (Diptera: Culicidae) catches using biogents sentinel traps during a ‘rear and release’ program: implications for designing surveillance programs. J Med Entomol. 2019;56(4):1102–11.

    Article  PubMed  Google Scholar 

  27. 27.

    Hesterberg TC. Estimates and confidence intervals for importance sampling sensitivity analysis. Math Comput Model. 1996;23(8):79–85.

    Article  Google Scholar 

  28. 28.

    Robert C, Casella G. Monte Carlo statistical methods. New York: Springer Science & Business Media; 2013.

  29. 29.

    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 

  30. 30.

    Costero A, Scott TW, Edman JD, Clark GG. Life table study of Aedes aegypti (Diptera: Culicidae) in Puerto Rico fed only human blood versus blood plus sugar. J Med Entomol. 1998;35(5):809–13.

    CAS  Article  PubMed  Google Scholar 

  31. 31.

    Sowilem MM, Kamal HA, Khater EI. Life table characteristics of Aedes aegypti (Diptera: Culicidae) from Saudi Arabia. Trop Biomed. 2013;30(2):301–14.

    PubMed  PubMed Central  Google Scholar 

  32. 32.

    Pagendam D, Snoad N, Yang W-H, Segoli M, Ritchie S, Trewin B, et al. Improving estimates of Fried’s Index from mating competitiveness experiments. J Agric Biol Environ Stat. 2018;23(4):446–62.

    Article  Google Scholar 

  33. 33.

    Doob JL. Markoff chains--denumerable case. Trans Am Math Soc. 1945;58(3):455–73.

    Article  Google Scholar 

  34. 34.

    Fried M. Determination of sterile-insect competitiveness. J Econ Entomol. 1971;64(4):869–72.

    Article  Google Scholar 

  35. 35.

    Gillespie DT. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comput Phys. 1976;22(4):403–34.

    CAS  Article  Google Scholar 

  36. 36.

    Ross JV, Pagendam DE, Pollett PK. On parameter estimation in population models II: multi-dimensional processes and transient dynamics. Theor Popul Biol. 2009;75(2):123–32.

    CAS  Article  PubMed  Google Scholar 

  37. 37.

    Dye C. Competition amongst larval Aedes aegypti: the role of interference. Ecol Entomol. 1984;9(3):355–7.

    Article  Google Scholar 

  38. 38.

    Johnson BJ, Ritchie SA. The siren’s song: exploitation of female flight tones to passively capture male Aedes aegypti (Diptera: Culicidae). J Med Entomol. 2015;53(1):245–8.

    Article  Google Scholar 

  39. 39.

    Johnson BJ, Mitchell SN, Paton CJ, Stevenson J, Staunton KM, Snoad N, et al. Use of rhodamine B to mark the body and seminal fluid of male Aedes aegypti for mark-release-recapture experiments and estimating efficacy of sterile male releases. PLoS Negl Trop Dis. 2017;11(9):e0005902.

    Article  CAS  PubMed  Google Scholar 

  40. 40.

    Hoffmann AA, Montgomery BL, Popovici J, Iturbe-Ormaetxe I, Johnson PH, Muzzi F, et al. Successful establishment of Wolbachia in Aedes populations to suppress dengue transmission. Nature. 2011;476(7361):454–7.

    CAS  Article  PubMed  Google Scholar 

  41. 41.

    Turelli M, Hoffmann AA. Rapid spread of an inherited incompatibility factor in Californian Drosophila. Nature. 1991;353(6343):440–2.

    CAS  Article  PubMed  Google Scholar 

  42. 42.

    Turelli M. Evolution of incompatibility-inducing microbes and their hosts. Evolution. 1994;48(5):1500–13.

    Article  PubMed  Google Scholar 

  43. 43.

    Barton NH, Turelli M. Spatial waves of advance with bistable dynamics: cytoplasmic and genetic analogues of Allee effects. Am Nat. 2011;178(3):E48–75.

    CAS  Article  PubMed  Google Scholar 

  44. 44.

    Pagendam D, Trewin B, Snoad N, Ritchie S, Hoffmann A, Staunton K, et al. Wolbachia establishment risk - model simulations. In: CSIRO, editor. 1.0 ed. CSIRO Data Access Portal: CSIRO; 2020.

  45. 45.

    Murray JV, Jansen CC, De Barro P. Risk associated with the release of Wolbachia-infected Aedes aegypti mosquitoes into the environment in an effort to control dengue. Front Public Health. 2016;4:43.

    Article  PubMed  Google Scholar 

  46. 46.

    Kittayapong P, Ninphanomchai S, Limohpasmanee W, Chansang C, Chansang U, Mongkalangoon P. Combined sterile insect technique and incompatible insect technique: the first proof-of-concept to suppress Aedes aegypti vector populations in semi-rural settings in Thailand. PLoS Negl Trop Dis. 2019;13(10):e0007771-e.

    Article  Google Scholar 

  47. 47.

    Crawford JE, Clarke DW, Criswell V, Desnoyer M, Cornel D, Deegan B, et al. Efficient production of male Wolbachia-infected Aedes aegypti mosquitoes enables large-scale suppression of wild populations. Nat Biotechnol. 2020;38(4):482–92.

    CAS  Article  PubMed  Google Scholar 

  48. 48.

    Focks DA, Haile DG, Daniels E, Mount GA. Dynamic life table model for Aedes aegypti (Diptera: Culicidae): analysis of the literature and model development. J Med Entomol. 1993;30(6):1003–17.

    CAS  Article  PubMed  Google Scholar 

  49. 49.

    Focks DA, Haile DG, Daniels E, Mount GA. Dynamic life table model for Aedes aegypti (Diptera: Culicidae): simulation results and validation. J Med Entomol. 1993;30(6):1018–28.

    CAS  Article  PubMed  Google Scholar 

Download references

Acknowledgements

Not applicable

Funding

This work was supported by the Australian National Health and Medical Research Council (NHMRC 1082127).

Author information

Affiliations

Authors

Contributions

Made substantial contributions to the conception: D.P., B.T., N.S., S.R., A.H., K.S., C.P. and N.B. Made substantial contributions to the analysis or interpretation of the data: D.P., B.T., S.R., A.H. and K.S. Made substantial contributions to the creation of software: D.P. Made substantial contributions to or substantively revised the drafts: D.P., B.T., N.S., S.R., A.H., K.S., C.P. and N.B. All authors read and approved the final manuscript.

Corresponding author

Correspondence to B. J. Trewin.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Competing interests

The authors declare that they have 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: Additional Analyses, Figures and Tables. A1. Modelling Population Parameter Ranges. Table A1.1

Primary model parameters and the lower and upper bounds of uniform prior distributions and references used to help define these. A2. Importance Sampling Estimates of Establishment and Elimination Probabilities. Table A2.1-A2.2 Importance sampling estimates for wAlbB establishment probabilities under the 5:1 and 15:1 overflooding ratio. A3. Additional File A Figures. Fig. A3.1 The probability of seeing one or more blocks with a wAlbB establishment event for increasingly large numbers of treated blocks at different block-level establishment probabilities: 0.1 (red); 0.01 (green); 0.001 (blue) and 0.0001 (black). Fig. A3.2 Estimates of wAlbB establishment probabilities across different wAlbB initial proportions.

Additional file 2: Model Description. B1. Conceptualising Mosquito Life Stages. Fig. B1.1.

Schematic outlining the states within the MPP model. B2. Derivation of Secondary Parameters from the Wildtype Equilibrium. B2.1 Enforcing Parametric Constraints. B2.2 Numbers of Future Adults (Juveniles) at Equilibrium. B2.3 Population Ceiling (Imax) for Future Adults. B2.4 Per Capita Mating Rate of Unmated Females. B3. Background Theory: Markov Population Process. B4. Population Dynamics and Model Transition Rates. B4.1 Birth of an Immature Future Adult into a Subpopulation. B4.2 Development of Immature Future Adults in Subpopulation X. B4.3 Emergence of Adults Within a Subpopulation. B4.4 Mating of Unmated Females Within a Subpopulation. B4.5 Death of Adults Within a Subpopulation.

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

Pagendam, D.E., Trewin, B.J., Snoad, N. et al. Modelling the Wolbachia incompatible insect technique: strategies for effective mosquito population elimination. BMC Biol 18, 161 (2020). https://doi.org/10.1186/s12915-020-00887-0

Download citation

Keywords

  • Incompatible insect technique
  • Wolbachia
  • Establishment risk
  • Elimination
  • Stochastic model
  • Simulation