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

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:

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?

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

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

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–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

Term	Definition
Constant release strategy	A release strategy where the same number of Wolbachia-infected mosquitoes are released at every release event for the duration of the IIT programme. This strategy does not use population monitoring to alter the number of individuals released. The number of Wolbachia-infected mosquitoes released per event is equal to the overflooding ratio multiplied by the male population at the commencement of the IIT programme.
Adaptive release strategy	A release strategy where the population size is assumed to be known at each release event. The number of Wolbachia-infected individuals released per event is equal to the overflooding ratio multiplied by the current male population.
Crude adaptive release strategy	A release strategy where releases are identical to the constant release strategy, but in three phases. In the first phase, the number of Wolbachia-infected mosquitoes released per release event is equal to the overflooding ratio multiplied by the male population at the commencement of the IIT programme. The second phase begins once the population reaches 50% of its initial size, at which point, the numbers of mosquitoes released per event is halved. Phase 3 begins when the population reaches 10% of its initial size, at the numbers released per event are reduced to 10% of that used in the first phase.
Eliminated	A wild-type population is considered to have been eliminated when there are no wild-type individuals alive in any of the life stages present in the model (i.e. no adults or future adults) at the end of the simulation.
Established	A Wolbachia population is (conservatively) defined to be established if the proportion of adults in the population that are infected with Wolbachia at the end of the simulation is greater than 20%.
Female contamination probability (FCP)	The probability that each mosquito released into a population at a release event is a female. The release of females may occur due to errors in the sex separation process and are modelled as events that follow independent, Bernoulli-distributed random variables for each individual.
Future adults	An abstract pool of individuals used in the model to represent non-adult life forms that will, with certainty, survive to become adults. See Additional File B: B2.1-B2.3 for further details.
IIT endpoint	The final state of a population that has been subjected to a Wolbachia-IIT programme. There are four possible endpoints: (i) elimination of the wild-type mosquito population, (ii) establishment of Wolbachia into the population, (iii) indeterminate Wolbachia negative and (iv) indeterminate Wolbachia positive. These four outcomes are defined in greater detail below.
Indeterminate Wolbachia positive	A population is considered to have reached an IIT endpoint of indeterminate Wolbachia positive if it is neither eliminated nor established, but where there are one or more individuals in one of the life stages (i.e. adults or future adults) that are infected with Wolbachia at the end of the simulation (i.e. 365 days post-final release).
Indeterminate Wolbachia negative	A population is considered to have reached an IIT endpoint of indeterminate Wolbachia negative if it is neither eliminated nor established and there are no individuals of any life stages (i.e. no adults or future adults) that are infected with Wolbachia at the end of the simulation (i.e. 365 days post-final release).
Overflooding ratio	The number of Wolbachia males released into the population at each release event divided by some baseline population (e.g. the initial wild-type male population size).
Unstable equilibrium threshold (UET)	The frequency of Wolbachia infection in a mosquito population above which a particular strain is likely to spread [15].
Wild-type population	A population consisting entirely of individuals that are not infected with Wolbachia.

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–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 I, I, …, I, where X ∈ {wild, wAlb}; adult males denoted M; unmated females denoted F; and mated females denoted F, 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.

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” (I), 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, 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 or wAlbB with probability . In other words, the rate at which females are mated by males is not density-dependent, but is mate-type. The parameter c 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 F and F can give birth to wAlbB future adults. Each mated female adds new “future adults” to state I at rate , where 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 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-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:

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.

Fig. 2.

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

Does the release of a small number of females during a Wolbachia-IIT 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 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.

Fig. 3.

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

What types of release strategies are most effective?

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

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.

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), 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.

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.