 Research article
 Open Access
 Published:
Modelling the Wolbachia incompatible insect technique: strategies for effective mosquito population elimination
BMC Biology volume 18, Article number: 161 (2020)
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 vectorborne diseases such as dengue, chikungunya and Zika. Successful implementation of this biological control strategy relies on highfidelity separation of male from female insects in mass production systems for inundative release into landscapes. Processes for sexseparating mosquitoes are typically errorprone and laborious, and IIT programmes run the risk of releasing Wolbachiainfected females and replacing wild mosquito populations.
Results
We introduce a simple Markov population process model for studying mosquito populations subjected to a WolbachiaIIT 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 everdecreasing wild population, thus reducing the risk of releasing Wolbachiainfected females while reducing costs.
Conclusions
While very highfidelity sex separation is required to avoid establishment, release programmes tend to be robust to the release of a small number of Wolbachiainfected females. These findings will inform and enhance the next generation of WolbachiaIIT 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 areawide 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 mosquitoborne 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 Wolbachiainfected male mosquitoes that are incapable of producing viable offspring after mating with a wildtype 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 Wolbachiainfected females. However, in crosses with wildtype 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 WolbachiaIIT process [12, 13]. Whilst new methods can have a high degree of accuracy, occasional errors in sorting systems can lead to the release of Wolbachiainfected females. For IIT programmes that rely on Wolbachiarelated 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 Wolbachiainfected 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 discretetime 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], spatiotemporal advectiondiffusionreaction models with multiple life stages [19] and complex agentbased 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 WolbachiaIIT 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 wildtype mosquitoes, but this process is not fitforpurpose in SIT/IIT programmes, since it does not allow for a declining population.
Markov population processes (MPPs) [23] are a class of continuoustime Markov chain models suitable for exploring the effects of WolbachiaIIT 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 largescale Wolbachia IIT programme in North Queensland, Australia [26], it was necessary to explore hypothetical population replacement scenarios where Wolbachiainfected females were released unintentionally. Here, we use an MPP model to address key questions related to the effective implementation of a WolbachiaIIT population suppression programme, namely:

(i)
What degree of fidelity in sex separation is required to ensure that a WolbachiaIIT population suppression programme does not result in Wolbachia establishment and the replacement of the wildtype (noninfected) population?

(ii)
Does the release of a small number of females during a WolbachiaIIT population suppression programme imply a high probability of Wolbachia establishment?

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

(iv)
How does the establishment probability vary with the initial density of the Wolbachia infection in the population?
We recognise the important tradeoff 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 WolbachiaIIT 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 WolbachiaIIT 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 wellmixed 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 WolbachiaIIT programmes use adult equilibrium populations (prior to wAlbB releases) in the range 100 to 300 wildtype 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 2year 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 WolbachiaIIT programme by calculating the proportions of simulations that reach these four IIT endpoints, with the elimination endpoint being the obvious goal.
We simulated WolbachiaIIT 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 nextgeneration 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) wildtype 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 WolbachiaIIT 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 wildtype and the wAlbB mosquito populations, there are k states which represent “future adults” (see Additional File B: B2.1B2.3) denoted I_{X, 1}, I_{X, 2}, …, I_{X, k}, where X ∈ {wild, wAlb}; adult males denoted M_{X}; unmated females denoted F_{X}; and mated females denoted F_{X × Y}, where Y ∈ {wild, wAlb}. The precise mathematical details of the MPP are provided in Additional File B: B1B4, but we outline the general mechanics below. The model (as an installable R package) is publicly available at https://github.com/dpagendam/debugIIT.
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 gammadistributed, 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 densitydependent 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” (I_{X, 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 I_{X, 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 wildtype 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 densitydependent, but is matetype. The parameter c_{wAlb} is the mating competitiveness coefficient or Fried’s Index of wAlbBinfected individuals relative to wildtype individuals [32]. We assume that CI between wAlbB males and wildtype females is 100% effective, so that only F_{wild × wild} females can give birth to wildtype future adults. Mated females of type F_{wAlb × wild} and F_{wAlb × wAlb} can give birth to wAlbB future adults. Each mated female adds new “future adults” to state I_{X, 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 I_{max} is a parameter that acts as a population ceiling on the number of future adults allowed within the population (see Additional File B: B4.1B4.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 densitydependent 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, I_{max}, 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 densityindependent 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 premated 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: B1B4 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, p_{mated} 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 wildtype 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 wAlbBinfected 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 wildtype 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 5year 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 wildtype populations were completely replaced by the wAlbBinfected 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 E52670 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:

(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.

(ii)
Does the release of a small number of females during a WolbachiaIIT population suppression programme imply a high probability of Wolbachia establishment?
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 WolbachiaIIT 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}.

(iii)
What types of release strategies are most effective?
To ensure a very low probability (< 0.01) of any blocklevel 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).

(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 wildtype equilibrium of 400 individuals (i.e. K_{wild} = 400), we simulated trajectories with different initial proportions of wAlbB and wildtype individuals. We display the first ten simulated population trajectories summarised as the proportion of adults infected with wAlbB for both frequency scenarios where no prerelease 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 prerelease mating and between 10 and 15% where prerelease 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].
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 welldeveloped implementation pathway, with lower barriers to regulatory approval and high community acceptance [45]. However, the release of Wolbachiainfected 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 WolbachiaIIT programmes and highlights the benefits of adopting stochastic rather than deterministic modelling frameworks. Importantly, our model highlights the attributes of WolbachiaIIT programmes that yield high probabilities of eliminating a population without Wolbachia establishment: (i) very highfidelity sex separation and (ii) overflooding ratios and release strategies that employ only moderate numbers of Wolbachiainfected 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 longterm outlook when conducting such IIT programmes. It follows that releasing an excess of Wolbachiainfected 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 areawide 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 densityindependent and only the mate type is densitydependent (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, realtime 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 densitydependent 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, noninteracting 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 blocklevel 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 blocklevel 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 Wolbachiainfected 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 postrelease 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 premated wAlbB individuals, the majority of simulations resulted in establishment where the initial wAlbB frequency was above 20%, and (ii) for wAlbB individuals without premating, 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 nonoverlapping 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 densitydependent 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 highfidelity 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 nonacademics: 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.
Dyck VA, Hendrichs J, Robinson AS. Sterile insect technique: principles and practice in areawide integrated pest management. Dordrecht: Springer; 2006.
 2.
VargasTerá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 areawide integrated pest management. Dordrecht: Springer Netherlands; 2005. p. 629–50.
 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.
 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.
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.
 6.
Kittayapong P, Kaeothaisong No, 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.
 7.
Laven H. Eradication of Culex pipiens fatigans through cytoplasmic incompatibility. Nature. 1967;216(5113):383–4.
 8.
Helinski MEH, Parker AG, Knols BGJ. Radiation biology of mosquitoes. Malar J. 2009;8(2):S6.
 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 SubSaharan Africa: recommendations of a Scientific Working Group^{†}. Am J Trop Med Hyg. 2018;98(6_Suppl):1–49.
 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.
 11.
Alphey L, Benedict M, Bellini R, Clark GG, Dame DA, Service MW, et al. Sterileinsect methods for control of mosquitoborne diseases: an analysis. Vector Borne Zoonotic Dis. 2010;10(3):295–311.
 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 southwestern Indian Ocean. PLoS Negl Trop Dis. 2011;5(12):e1440.
 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.
 14.
Xi Z, Khoo CCH, Dobson SL. Wolbachia establishment and invasion in an Aedes aegypti laboratory population. Science. 2005;310(5746):326–8.
 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.
 16.
Knipling EF. Possibilities of insect control or eradication through the use of sexually sterile males. J Econ Entomol. 1955;48(4):459–62.
 17.
Knipling EF. Sterilemale 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.
Hancock PA, Sinkins SP, Godfray HCJ. Population dynamic models of the spread of Wolbachia. Am Nat. 2011;177(3):323–33.
 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.
 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.
 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.
 22.
Jansen VAA, Turelli M, Godfray HCJ. Stochastic spread of Wolbachia. Proc Biol Sci. 2008;275(1652):2769–76.
 23.
Kingman JFC. Markov population processes. J Appl Probab. 1969;6(1):1–18.
 24.
Muir LE, Kay BH. Aedes aegypti survival and dispersal estimated by markreleaserecapture in northern Australia. Am J Trop Med Hygiene. 1998;58(3):277–82.
 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.
 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.
 27.
Hesterberg TC. Estimates and confidence intervals for importance sampling sensitivity analysis. Math Comput Model. 1996;23(8):79–85.
 28.
Robert C, Casella G. Monte Carlo statistical methods. New York: Springer Science & Business Media; 2013.
 29.
Hancock PA, White VL, Callahan AG, Godfray CHJ, Hoffmann AA, Ritchie SA. Densitydependent population dynamics in Aedes aegypti slow the spread of wMel Wolbachia. J Appl Ecol. 2016;53(3):785–93.
 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.
 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.
 32.
Pagendam D, Snoad N, Yang WH, 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.
 33.
Doob JL. Markoff chainsdenumerable case. Trans Am Math Soc. 1945;58(3):455–73.
 34.
Fried M. Determination of sterileinsect competitiveness. J Econ Entomol. 1971;64(4):869–72.
 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.
 36.
Ross JV, Pagendam DE, Pollett PK. On parameter estimation in population models II: multidimensional processes and transient dynamics. Theor Popul Biol. 2009;75(2):123–32.
 37.
Dye C. Competition amongst larval Aedes aegypti: the role of interference. Ecol Entomol. 1984;9(3):355–7.
 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.
 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 markreleaserecapture experiments and estimating efficacy of sterile male releases. PLoS Negl Trop Dis. 2017;11(9):e0005902.
 40.
Hoffmann AA, Montgomery BL, Popovici J, IturbeOrmaetxe I, Johnson PH, Muzzi F, et al. Successful establishment of Wolbachia in Aedes populations to suppress dengue transmission. Nature. 2011;476(7361):454–7.
 41.
Turelli M, Hoffmann AA. Rapid spread of an inherited incompatibility factor in Californian Drosophila. Nature. 1991;353(6343):440–2.
 42.
Turelli M. Evolution of incompatibilityinducing microbes and their hosts. Evolution. 1994;48(5):1500–13.
 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.
 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.
Murray JV, Jansen CC, De Barro P. Risk associated with the release of Wolbachiainfected Aedes aegypti mosquitoes into the environment in an effort to control dengue. Front Public Health. 2016;4:43.
 46.
Kittayapong P, Ninphanomchai S, Limohpasmanee W, Chansang C, Chansang U, Mongkalangoon P. Combined sterile insect technique and incompatible insect technique: the first proofofconcept to suppress Aedes aegypti vector populations in semirural settings in Thailand. PLoS Negl Trop Dis. 2019;13(10):e0007771e.
 47.
Crawford JE, Clarke DW, Criswell V, Desnoyer M, Cornel D, Deegan B, et al. Efficient production of male Wolbachiainfected Aedes aegypti mosquitoes enables largescale suppression of wild populations. Nat Biotechnol. 2020;38(4):482–92.
 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.
 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.
Acknowledgements
Not applicable
Funding
This work was supported by the Australian National Health and Medical Research Council (NHMRC 1082127).
Author information
Affiliations
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
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.1A2.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 blocklevel 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 (I_{max}) 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.
About this article
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/s12915020008870
Received:
Accepted:
Published:
Keywords
 Incompatible insect technique
 Wolbachia
 Establishment risk
 Elimination
 Stochastic model
 Simulation