Analysis of a household-scale model for the invasion of Wolbachia into a resident mosquito population

TL;DR

The paper addresses Wolbachia invasion in Ae. aegypti at the household scale by formulating a stochastic continuous-time Markov chain model with a natural carrying-capacity bound. It links deterministic bistability and the invasion threshold to stochastic outcomes via quasi-stationary distributions, enabling computation of invasion probabilities, times to invasion, and reversion risks under two vertical-transmission regimes. Key findings show that deterministic thresholds translate into probabilistic invasion outcomes (e.g., only ~23% of releases displaced wildtype under certain parameters; larger releases are needed for 90% likelihood of Wolbachia fixation), and that reversion due to imperfect vertical transmission can substantially shorten effective household protection. The approach provides practical metrics for release design and highlights how stochasticity in small populations can alter control strategies, with future work extending to connected households and disease dynamics.

Abstract

In areas infested with Aedes aegypti mosquitoes it may be possible to control dengue, and some other vector-borne diseases, by introducing Wolbachia-infected mosquitoes into the wildtype population. Thus far, empirical and theoretical studies of Wolbachia release have tended to focus on the dynamics at the community scale. However, Ae. aegypti mosquitoes typically dwell in and around the same houses as the people they bite and it can be insightful to explore what happens at the household scale where small population sizes lead to inherently stochastic dynamics. Here we use a continuous-time Markov framework to develop a stochastic household model for small populations of wildtype and Wolbachia-infected mosquitoes. We investigate the transient and long term dynamics of the system, in particular examining the impact of stochasticity on the Wolbachia invasion threshold and bistability between the wildtype-only and Wolbachia-only steady states previously observed in deterministic models. We focus on the influence of key parameters which determine the fitness cost of Wolbachia infection and the probability of Wolbachia vertical transmission. Using Markov and matrix population theory, we derive salient characteristics of the system including the probability of successful Wolbachia invasion, the expected time until invasion and the probability that a Wolbachia-infected population reverts to a wildtype population. These attributes can inform strategies for the release of Wolbachia-infected mosquitoes. In addition, we find that releasing the minimum number of Wolbachia-infected mosquitoes required to displace a resident wildtype population according to the deterministic model, only results in that outcome about 20% of the time in the stochastic model; a significantly larger release is required to reach a steady state composed entirely of Wolbachia-infected mosquitoes 90% of the time.

Figures (34)

• Figure 1: Larval intraspecific competition functions $F(N)$ and $G(N)$; see equations \ref{['eq:hughes_F']} and \ref{['eq:qu']}. Parameter values are $C=30$, $k=0.3$ and $h=0.76$.
• Figure 2: Stability and attraction of the steady states of system \ref{['eq:hughes']}--\ref{['eq:hughes2']}, with larval density function \ref{['eq:hughes_F']}, depending on the fitness cost of Wolbachia infection $1-\phi$ and the initial proportion of the mosquito population that is Wolbachia-infected $N_w(0)/N_0$. Here $u,v=1$ and $N_0 = 10$. Indicated on the diagram are the regions where the system is bistable and either the wildtype-only ($E1$) or Wolbachia-only ($E2$) steady state is attracting, and where the unique non-trivial steady state is wildtype-only, and stable. The invasion threshold under the bistable region is labelled in blue text. All of our stability analysis plots are produced using the open source software bSTAB stender2022bstab.
• Figure 3: Stability and attraction of steady states of system \ref{['eq:hughes']}--\ref{['eq:hughes2']} depending on the initial numbers of female wildtype $N_m(0)$ and Wolbachia-infected $N_w(0)$ mosquitoes. Here $u,v=1$ and $1-\phi=0.15$, with larval density function \ref{['eq:hughes_F']}. Indicated on the diagram are the regions where the wildtype-only ($E1$) or Wolbachia-only ($E2$) steady states are attracting. The grey triangle covers populations outside our upper bound on the total female mosquito population size. The dashed blue line indicates where $N_m(0)+N_w(0)=10$, corresponding to $1-\phi = 0.15$ in Figure \ref{['fig:basin_attr_Dye']}.
• Figure 4: Schematic of the communicating classes in the $3$ mosquito model. The states in the mixed class $\mathcal{S}_{3}$ are shaded green, the states in the wildtype-only class $\mathcal{S}_{1}$ are orange and the states in the Wolbachia-only class $\mathcal{S}_{2}$ are blue.
• Figure 5: Probability distribution of all transient states over time, conditioned on non-extinction and initial state being $(1,1)$ with probability $1$. Parameter values are as stated in Subsection \ref{['subsec:parameterisation']}. (a) Over $575$ days and (b) enlargement of the first $100$ days. The dashed lines indicate the QSD values for the wildtype-only states.