Mathematical analysis of a Wolbachia invasive model with imperfect maternal transmission and loss of Wolbachia infection.

Arboviral infections, especially dengue, continue to cause significant health burden in their endemic regions. One of the strategies to tackle these infections is to replace the main vector agent, Ae. aegypti, with the ones incapable of transmitting the virus. Wolbachia, an intracellular bacterium, has shown promise in achieving this goal. However, key factors such as imperfect maternal transmission, loss of Wolbachia infection, reduced reproductive capacity and shortened life-span affect the dynamics of Wolbachia in different forms in the Ae. aegypti population. In this study, we developed a Wolbachia transmission dynamic model adjusting for imperfect maternal transmission and loss of Wolbachia infection. The invasive reproductive number that determines the likelihood of replacement of the Wolbachia-uninfected (WU) population is derived and with it, we established the local and global stability of the equilibrium points. This analysis clearly shows that cytoplasmic incompatibility (CI) does not guarantee establishment of the Wolbachia-infected (WI) mosquitoes as imperfect maternal transmission and loss of Wolbachia infection could outweigh the gains from CI. Optimal release programs depending on the level of imperfect maternal transmission and loss of Wolbachia infection are shown. Hence, it is left to decision makers to either aim for replacement or co-existence of both populations.

Introduction

Vector borne diseases such as dengue, Zika, chikungunya and yellow fever are of global health concern. For instance, dengue has a widespread geographical distribution with around 3.9 billion people at risk and an annual estimate of 390 million new dengue infections (Bhatt et al., 2013; Kyle & Harris, 2008). The major vector responsible for the transmission of dengue and other arboviral infections is the female Ae. aegypti mosquito. Although, the female Ae. albopictus mosquitoes also contributes. The risk of mortality due to dengue infection is low but is modified by the serotype of the infecting dengue and an individual's infection history (particularly their immune response to different serotypes) (Kyle & Harris, 2008). Despite the low risk of mortality, the large number of confirmed dengue cases and associated morbidity make dengue a substantial contributor to the global health burden. The World Health Organization (WHO) global target for dengue by 2020 is to reduce morbidity and mortality by at least 25% and 50% respectively (Organizationet al., 2012). Integrated vector management is one among many potential control strategies being considered. Controlling the mosquito vectors appears to be promising but it comes with challenges and great cost (Manrique-Saide et al., 2015; Ooi, Goh, & Gubler, 2006).

Rather than preventing human-vector contacts, replacing the population of Ae. aegypti mosquitoes with another variant that is incapable of viral transmission has been successfully applied to reduce dengue infections (Hoffmann et al., 2011) and the approach appears promising for other mosquito-borne infections such as chikungunya, malaria, West-Nile virus, and zika virus (Dutra et al., 2016; Glaser & Meola, 2010; Gomes & Barillas-Mury, 2018; Moreira et al., 2009). Wolbachia, an intracellular insect bacterium, has the capacity to inhibit dengue virus proliferation inside the Ae. aegypti mosquitoes and can spread via maternal (vertical) transmission (Moreira et al., 2009; Turley, Moreira, O’Neill, & McGraw, 2009). Depending on the strain of Wolbachia, mosquitoes infected with Wolbachia have a reproductive advantage over those uninfected via cytoplasmic incompatibility (CI) the mechanism that prevents the embryo maturing following mating between Wolbachia infected (WI) males and Wolbachia uninfected (WU) females (Ant, Herd, Geoghegan, Hoffmann, & Sinkins, 2018; Turelli & Hoffmann, 1995). This advantage alone may not guarantee that WI mosquitoes will replace the Ae. aegypti population, as Wolbachia infection leads to a fitness cost to its host (Fine, 1978; Walker et al., 2011) and also, there are reports regarding Wolbachia infection and the loss of cytoplasmic incompatibility (Fine, 1978; Ross, Ritchie, Axford, & Hoffmann, 2019; Turelli & Hoffmann, 1995) as a result of unfavorable conditions that lead to loss of Wolbachia infection in infected adult mosquitoes. Another factor that could prevent WI mosquitoes from dominating the Ae. aegypti population is imperfect maternal transmission (Fine, 1978; Turelli, 2010; Turelli & Hoffmann, 1995; Yeap et al., 2011). Hence, having a full understanding of the interplay between key parameters in Wolbachia introduction is necessary to ensure the success of the strategy if it is to be used on a large scale.

Mathematical modelling plays a significant role in understanding the impact of variables involved in the dynamics of a particular infectious disease and has been used in the decision-making process that guides the application of some typical control strategies (Martcheva, 2015). Different mathematical models have been developed to simulate the introduction of Wolbachia into Ae. aegypti populations (Campo-Duarte, Vasilieva, Cardona-Salgado, & Svinin, 2018; Crain et al., 2011; Ndii, Hickson, & Mercer, 2012; Rafikov, Meza, Correa, & Wyse, 2019; Schraiber et al., 2012; Xue, Manore, Thongsripong, & Hyman, 2017; Zheng, Tang, Yu, & Qiu, 2018), with each specifying conditions that enable WI mosquitoes to dominate. Caspari and Watson (Caspari & Watson, 1959) demonstrated the importance of cytoplasmic incompatibility on the population replacement between WU and WI mosquitoes. Ndii et al. developed a deterministic compartmental model for the competition between the two mosquitoes populations and derived the steady-state solutions showing key parameters that could influence the competition between the two populations (Ndii et al., 2012). Xue et al. adopted similar approach as (Ndii et al., 2012) by incorporating sex structure into the compartmental models and showed that the endemic Wolbachia steady-state solution can be established by releasing a sufficiently large number of Wolbachia infected mosquitoes (Xue et al., 2017). Using ordinary differential equations to model the competitions between WU and WI mosquitoes, Zheng et al. showed that the successful replacement of WU mosquitoes with WI ones would depend on the strains of Wolbachia used and require a careful release design (Zheng et al., 2018). The idea of designed release methods was further emphasised by Qu et al. where they extended the model in (Xue et al., 2017) to include female mosquitoes mating once (Qu, Xue, & Hyman, 2018). The model of Li and Liu places emphasis on the combinations of birth and death rate functions, Wolbachia strain and the number of WI mosquitoes released (Li & Liu, 2017). All these modelling works pointed to the possibility of WI mosquitoes replacing the uninfected ones.

In this paper, we consider the impacts of imperfect maternal transmission and loss of Wolbachia infection by investigating the asymptotic dynamics of the Wolbachia invasive model and determining the necessary and sufficient conditions for Wolbachia invasion. These two factors (imperfect maternal transmission and loss of Wolbachia infection) have not been considered by previous models and the derivation of both asymptotic and global stability of the possible equilibrium points are the novel results of this work. With the derivation of the global stability, appropiate control strategies can be adopted to ensure that WI mosquitoes can replace uninfected ones or at least become more abundant than the uninfected ones. We consider these strategies via optimal control.

Model formulation

We consider the Ae. aegypti mosquito population that is responsible for the transmission of most arboviral infections, in particular dengue virus. Similar to (Xue et al., 2017), the Ae. aegypti mosquito population is divided into two major subpopulations: those with Wolbachia infection (w); and those without Wolbachia infection (). We further denote the number of mosquitoes in the aquatic stage (egg, larvae, and pupae) as for WU and for WI, and also by considering the experimental work in (Arrivillaga, 2004), we set the ratio of adult male to female mosquito to one and denote adult WU female mosquitoes as and adult WI female mosquitoes as . Hence, our Wolbachia invasive model with imperfect maternal transmission and loss of Wolbachia infection is:

In equations (2.3), (2.4), we assumed logistic growth in the aquatic stage and a mosquito in the aquatic stage can die at the rate or maturates to become adult mosquito at the rate ψ. The adult mosquito can die at the rate if WU and if WI. Also, an adult WI mosquito losses their Wolbachia infection at the rate σ. There is a carry capacity K on the aquatic stage implying that is always non-negative.

Unlike the Wolbachia invasive model in (Xue et al., 2017), we assumed that the proportions of offspring due to imperfect maternal transmission are different between adult Ae . aegypti males and females of different Wolbachia infection status (Turelli & Hoffmann, 1995; Bian et al., 2013). We further include the possibility of a decline in the level of Wolbachia infection by allowing some WI to become WU at a constant per capital rate σ. The system of differential equations (2.2), (2.3), (2.1) explicitly includes the possibility of uninfected and infected offspring being produced by WI female mosquitoes (Hoffmann et al., 2011; Turelli, 2010; Turelli & Hoffmann, 1995; Yeap et al., 2011), through the and terms. This possibility is excluded in the model of Wolbachia introduction studied in (Ndii et al., 2012). However, they adjusted for leakage by assuming there is a waiting time before WI offspring mature into either WI or WU adults. It is biologically plausible that some offspring are born uninfected by WI females (Turelli & Hoffmann, 1995) and the approach by (Ndii et al., 2012) is another way of modelling the effect of losing incompatibility between WU females and WI males, which we have incorporated with the σ term (see Fig. 1). The descriptions of the parameters in the system of differential equations (2.2), (2.3), (2.1) are shown in Table 1.

Fig. 1

Schematic flow diagram of the The WU population (blue color) is produced as a result of mating between adult WU females and males, WI females and WU males, WI females and males and loss of Wolbachia infection by adult WI females. The WI population (red color) is produced by mating between adult WI females and males and cross-breeding between WI females and WU males. Due to CI, the offsprings as a result of mating between adult WU females and WI males are not viable.

Table 1

Parameter description and values for model (2.1–2.4).

Parameter	Description	Estimate[Range]	Unit	References
K	Carrying capacity of the aquatic stage	10^6[10^4,10^8]	Aquatic mosquito	Assumed
ϕ_w¯	Per capita egg laying rate for Wolbachia uninfected mosquitoes	13[12−18]	Eggs per day	(Hoffmann et al., 2014; McMeniman et al., 2009; McMeniman and O’Neill, 2010)
ϕ_w	Per capita egg laying rate for Wolbachia-infected mosquitoes	11[8−12]	Eggs per day	(Hoffmann et al., 2014; Walker et al., 2011)
ρ_1	The fraction of eggs that are WU as a result of mating between adult WI female and male mosquitoes	0.05[0−0.11]	Dimensionless	Walker et al. (2011)
ρ_2	The fraction of eggs that are WU as a result of mating between adult WU male and WI female mosquitoes	0.05[0−0.1]	Dimensionless	Walker et al. (2011)
σ	Per capita loss of Wolbachia infection	0.04[0−0.1]	Per day	Assumed
b	Fraction of eggs that are female	0.5[0.34−0.6]	Dimensionless	(Arrivillaga, 2004; Lounibos & Escher, 2008)
ψ	Per capita maturation rate	0.11[0.1−0.12]	Per day	(Hoffmann et al., 2014; Walker et al., 2011)
μ_a	Per capita aquatic death rate	0.02	Per day	Xue et al. (2017)
μ_w¯	Per capita death rate of WU mosquitoes	0.061[0.02−0.09]	Per day	(McMeniman et al., 2009; Styer, Minnick, Sun, & Scott, 2007)
μ_w	Per capita death rate of WI mosquitoes	0.068[0.03−0.14]	Per day	(Styer et al., 2007; Walker et al., 2011)

We analyse the Wolbachia invasive model (2.1–2.4) for the conditions that will enable the WI mosquitoes to propagate following their introduction into an Ae. aegypti population that is naïve to Wolbachia infection. This is done as follows. Given a system of autonomous ordinary differential equations,

The asymptotic behavior of the solutions of (2.5) starting near an equilibrium solution are determined by the eigenvalues of the associated Jacobian matrix defined as evaluated at . We adopt this approach in this paper to understand Wolbachia propagation in the Ae. aegypti population. Also, the simulations of this model using published parameter values were done in MATLAB R2017a (Release M(2017) The MathWorks Inc, Natick, MA, USA).

Model analysis

Wolbachia invasive model with and

When and in equations (2.2), (2.3), (2.1), we assumed imperfect maternal transmission is only between WI females and WU males and that there is no loss of Wolbachia infection in adult infected mosquitoes. The Wolbachia invasive model (2.1–2.4) with and is biologically meaningful (see Appendix A). That is, all solutions with non-negative initial conditions will remain non-negative for future times. Showing this for (0, 1] and is also straight-forward.

The Wolbachia invasive model (2.1)–(2.4) with and has four steady states: - where there are no mosquitoes; - where the WU mosquitoes dominate and leads to the extinction of infected ones; - where only the WI mosquitoes exist; and - where both WU and WI mosquitoes coexist. It is important for the control of arboviral infections that are transmitted by Ae. aegypti mosquitoes to determine the nature of these stability points.

No mosquitoes

The point is trivial but not interesting as it is not realistic. However, we can gain insights about the nature of this steady state solution by examining a special case when there is no interaction between WU and WI mosquitoes. We derivedand

Which are the thresholds that determine whether each population will persist or extinguish in the absence of interactions. The thresholds in equations (2.3), (2.4) are derived from the stability conditions of the associated Jacobian matrix when no interaction exists between the uninfected and infected mosquitoes. That is, individual populations do not depend on each other. Equivalent expressions were given in (Xue et al., 2017) for the dynamics that explicitly include the male mosquito compartments. Hence, for these models (2.1)–(2.4), the two populations are extinguished whenever and (see Fig. 2) as the reproductive terms cannot sustain the populations. Also, since the solutions are always non-negative for non-negative initial data, the solutions tend to the no-mosquito equilibrium point. However and except for the biological implications of using insecticides, applying insecticides and destroying breeding sites have been effective method in reducing mosquito populations (Amer & Mehlhorn, 2006).

Fig. 2

No mosquito equilibrium point. In this simulation, we set , , , , , , , , and .

WU mosquitoes-only

The WU equilibrium point is given as where

Hence is necessary for the existence of this equilibrium point, otherwise, there will be no WU mosquitoes. Let us define the invasive reproductive number associated with the number of secondary offspring that would be WI due to the introduction of a typical WI adult mosquito into a population of WU adult mosquitoes. In a similar way to (Xue et al., 2017), we derived as follows: the WI compartments can be divided into the rate of appearance of new mosquitoes with Wolbachia infection () and other transition rates such as progression into adult mosquitoes with Wolbachia infection and death rates :

Next, we introduce the matrices F and V with components and , where the represent the infected compartments and . Hence,and the next-generation matrix is

Hence the invasive reproductive number iswhere is the spectral radius of M. The factor shows the effect of the proportion of aquatic stage mosquitoes that are WI as a result of mating between WU male mosquitoes and WI female mosquitoes on the likelihood of the WI mosquitoes to replace the WU ones.

The Jacobian at is given as

The characteristic equation for this Jacobian is given aswhere

The equilibrium point is locally asymptotically stable whenever and . This implies for this case that WI mosquitoes will not spread following their introduction if those conditions are satisfied. The condition is the same as the only condition given by Ndii et al. (Ndii et al., 2012) for the stability of this point. We have additional condition () which states that the invasive reproductive rate of WI mosquitoes when introduced into a background of WU mosquitoes be less than one.

WI mosquitoes-only

The equilibrium point associated with WI mosquitoes only is

As pointed out earlier that WI population dies out when , the equilibrium point () is expected to be unstable when. Thus, the corresponding Jacobian is defined as

Hence, the characteristic equation iswhere , , , and . The eigenvalues of the quartic characteristic equation are negative or have negative real parts if and . The condition gives the fitness level of the adult WI that is sufficient for spread. As Wolbachia infection decreases the fitness of infected mosquitoes, external support will be needed for WI mosquitoes to propagate (Walker et al., 2011).

Following from the expression of the basic reproduction number (equation (3.10)), the equilibrium point for adult female mosquitoes can be written as

The expression above shows that the WI-mosquito-only equilibrium can exist when. The existence of endemic equilibria for is an indicator of a backward bifurcation in the conventional infectious diseases modelling papers (Cui, Sun, & Zhu, 2008; Dushoff, Huang, & Castillo-Chavez, 1998). However, it is interesting to know that this equilibrium point is unstable whenever and it is locally asymptotically stable even if , in as much as and . For , both the and equilibrium points are locally asymptotically stable in as much as for , and and for .

Provided , the WI mosquitoes only equilibrium point is globally asymptotically stable whenever and .

Proof. When then and this implies. Define a Lyapunov function V as

Differentiating equation (3.20) with respect to time, we have

Substituting the expression for the differential equations (2.3), (2.4) we have,where , and

From equation (3.22),and from equation (3.23),

Adding equations (2.3), (2.4) yields

Rearrangement and some manipulations give,

Thus,

Equation (3.28) can be re-written as

The second term of the right hand side of equation (3.29) is less than or equal to zero because . The first term is less than or equal to zero when . When , equation (3.29) is less than zero because . Hence, it follows from the Krasovkii-Lasalle Theorem (Krasovskii, 1959; LaSalle, 1960) that

It remains to show that . In this case, and . Hence, there exists a sufficiently small number and , such that and for all . It follows from equation (2.1) that for,where . Hence, by the comparison Theorem (Smith & Waltman, 1995) and letting

If then . Otherwise,and

Hence,

Thus, it is immediate that and . Hence, as . This concludes the Proof.

The condition was used in the derivation of equation (3.34) above. Also, the conditions of Theorem (3.1) make the only stable equilibrium point. We examined theorem (3.1) by simulating equations (2.2), (2.3), (2.4), (2.1) with parameters that satisfy these conditions (see Fig. 3).

Fig. 3

Global stability of 2.1)-(2.4). In this simulation, , , , , (we assumed external factor to compensate for decrease in fitness), , and the initial data are and The black dash line is the steady-state solution for .

Both mosquitoes

An interesting situation is to have both WI and uninfected mosquitoes in the Ae. aegypti population. In such case, we will want the majority of the mosquitoes to be infected with Wolbachia. For the systems of differential equations (2.2), (2.3), (2.4), (2.1), the co-existence equilibrium point is given aswhere, , , . From (3.36), it can be observed immediately that , , and either of these two conditions:

and

and

Must be true for the existence of this equilibrium point. The conditions show key parameter relationships for both WI and WU mosquitoes to sustain themselves. If , WU mosquitoes has a tolerable death rate that allows WI mosquitoes to survive and similar tolerance for . Hence all of the equilibrium points of the Wolbachia spread model (2.1)–(2.4) can co-exist when (Xue et al., 2017) and only the co-existence and WI-mosquito-only equilibrium points are locally stable when . To establish whether this co-existence equilibrium point is stable or not, we use the general Jacobian expressions (see appendix (Appendix B)) to derive its corresponding characteristic equation:where the coefficients are given by the following expressions:

One way to establish the nature of the equilibrium point is to apply the Lienard and Chipart criterion (Liénard & Chipart, 1914), or the popular Routh-Hurwitz Criteria (Martcheva, 2015). For this criterion, it is necessary and sufficient to show that the coefficients of the quartic equation (3.37) are greater than zero and that for the equilibrium point to be locally asymptotically stable. We investigate this equilibrium point by using a randomization method to show that there exists a parameter set for which the conditions above are satisfied and the quartic equation (3.37) has negative roots (i.e. the eigenvalues of the associated Jacobian have negative real parts). For the parameters in Table 1, we sample 10000 parameter combinations assuming uniform distributions for the ranges listed in Table 1 and check whether the equilibrium point is asymptotically stable or not when condition (1) or (2) is satisfied. We found for the set of parameters satisfying condition (1) that the equilibrium point is unstable. For condition (2), the equilibrium point is locally asymptotically stable but with unrealistic parameter set (Fig. 4a). Despite the unrealistic nature of the parameter set satisfying condition (2) above, that equilibrium point is globally asymptotically stable as demonstrated numerically (Fig. 4b). Also, Table 2 below lists the conditions for local asymptotic stability of the equilibrium points.

Fig. 4

The stability conditions of the co-existence equilibrium point. (a) The nature of the stability point changes with . The Routh- Hurtwitz conditions are satisfied for condition (2) of the co-existence equilibrium points. (b) We set , , , , , , to show its global stability. The initial data are and.

Table 2

Conditions for stability of the equilibrium points.

Equilibrium point	Stability conditions
E_1 (No Mosquitoes)	R_0w<1 and R_0w¯<1
E_2 (Only WU Mosquitoes)	R_0w|w¯<1 and R_0w¯>1
E_3 (Only WI Mosquitoes)	R_0w>1 and μ_w¯>ρ_2μ_w
E_4 (Both Mosquitoes)	ρ_2<1, R_0w|w¯>1,R_0w>1,R_0w¯>1
	and μ_w¯<ρ_2μ_w.

When and , the Wolbachia invasive model (2.1)–(2.4) has three steady state solutions: , where and with the expression for terms in defined later. Here, we do not have the WI-only mosquito equilibrium point because of the per capita loss of Wolbachia infection rate (σ) that always replenishes the WU population.

The adjusted invasive reproductive number is:

As expected, the expression for the invasive reproductive number shows that the loss of Wolbachia infection reduces by the term and in turn, reduces the rate at which the WI mosquitoes invade the WU population. The imperfect maternal transmission between adult WI mosquitoes does not affect the adjusted reproductive number (no term in equation (3.42)).

First, we investigate the impact of the individual reproduction numbers on the dynamics of the general Wolbachia invasive model (2.1)–(2.4). As before, when both and are less than one the solutions tend to the no-mosquitoes equilibrium point (Fig. 5a). Similar to subsection (3.1.2), the Jacobian for the steady-state solution isand the point is locally asymptotically stable whenever and . We state the following Theorem to show the conditions for global stability of.

Fig. 5

Simulation of the general 2.1) -(2.4) for different steady state solutions. (a) Here, we set , , and . (b) The general reproduction number () is 0.51, , , and .(c), , , and .

Whenever , the WU mosquito-only equilibrium point () is the only non-trivial equilibrium point whenever .

Proof. For the equilibrium point , we have;

And using Equations (2.3), (2.4), we have

From (3.44) and (3.45), equation (3.46) becomeswhere

Thus if and , the quadratic equation (3.47) has a positive solution: , withand

Hence, when and we get a negative solution for. Since our general model is biological meaningful (see Appendix A), cannot exist for these conditions. That leaves has the only non-trivial equilibrium point.

Fig. 5b shows the numerical demonstration of Theorem 3.2. The condition can be interpreted as a condition that determines the mosquito fitness advantage of WU mosquitoes over WI mosquitoes. Hence, when this condition is satisfied, WI mosquitoes cannot replace WU except there is a deliberate attempt to wipe-out all WU and replace them with WI mosquitoes. This is not realistic.

If , either of the following two conditions guarantees a positive solution for;

,

And if and leave as the only possible stable point in the positive quadrant . If this point is locally asymptotically stable and no other solutions exist in the plane then it is globally asymptotically stable for any positive initial condition (see Poincaré-Bendixson Trichotomy Theorem and Fig. 5c). The condition is the condition for backward bifurcation (Fig. 6a) and it implies that has two equilibria points with the one with higher locally asymptotically stable (Fig. 6b and c) with . The derivation of the conditions for local stability of the point when is shown in appendix Appendix C.

Fig. 6

The backward bifurcation and local stability of the steady-state solution. We vary and set values the following key parameters as: , , and . (a) Shows the backward bifurcation of the general model association with equilibrium points and . (b) and (c) establish the local stability of the equlibrium point using Routh-Hurwitz conditions - (b), , , , and (c),.

Optimal Wolbachia release problem

Similar to (Campo-Duarte et al., 2018), we considered two release strategies: constant release rate and variable release rate. In (Campo-Duarte et al., 2018), the focus is to wipeout WU mosquitoes and Rafikov et al. (Rafikov et al., 2019) strategy is to have more WI mosquitoes than WU mosquitoes. Our modelling work has shown that the focus could be different depending on the dynamics of the Wolbachia strain in the Ae. aegypti population.

Release strategy when and

For this case, the obvious strategy is to replace the WU mosquitoes with the infected ones. The parameter sets to ensure the conditions for the existence of the co-existence equilibrium in this case are unrealistic and it will be very difficult to achieve as there are limitations on the parameters that can be controlled. Adjusting equation (2.4) for the constant release rate strategy, we havewhere δ is the per capita release rate. We need to force and . With this, the only locally asymptotically stable point is WI-only mosquito point. The two prior conditions indicate that:

Here, we set , the value of in Table 1. However, it can be large as desired depending on the resource constraint. Using the values in Table 1, we have . This bound on δ is novel as we know the allowable range to acheive our goal. Figure (7) shows the application of the constant release rate program for different rates for 480 days (same as the release program in Townsville, Australia (O’Neill et al., 2018)) and fixing for different time frames. In Figure (7a and d), the goal is not achieved as WI mosquitoes are immediately wipeout and for Figure (7b, c, e, and f), it takes time for WI mosquitoes to establish themselves.

Fig. 7

The δ and time frames. With final time of 480 days as shown by the grey area, (a) , (b) , and (c) . Setting , we varied the final time (grey area) for (d) days, (e) days, and (f) days. The other parameter values are as in Table 1. Note that , and assuming that WU mosquitoes are in the WU-only equilibrium point at .

For the variable release rate, we follow similar approach as (Campo-Duarte et al., 2018) by introducing control variable . Hence, the variable release problem is an optimal control problem with constraint on the both end points:

Subject to equations (2.2), (2.3), (2.1), and

is the cost per unit time associated wih this strategy and is a balancing cost. The quadratic term is to ensure we have a regular optimal control. With the terminal constraint (4.5), the global stability of the WI-mosquito only equilibrium point (Theorem 3.1) only requires . Following the release program of wMel strain of Wolbachia in Townsville, Australia (O’Neill et al., 2018), stage 1 requires 14 months ( 480 days) and it cost per km2. Additionally, there is a cost associated with staff involved in the program. Here, we set and . and solve the optimal control problem (4.3)–(4.5) using the function space conjugate gradient algorithm (Edge & Powers, 1976). See Appendix D for the optimal characterization and algorithm for this problem. Since is unknown, we solve the optimal control problem (4.3)–(4.6) for with a step of 50 days and select the minimum objective function value. This approach is a modification of the algorithm described in chapter 7, section 7.4 of (Bryson & Ho, 1975) for solving an optimal control problem with an unspecified final time. Fig. 8 shows the optimal solutions for the variable release strategy. For this strategy, we do not need to wipe out all the WU mosquitoes to achieve our aim.

Fig. 8

The . (a) Objective function value for different final times. The minimal cost is at the final time, , as indicated by the pink line, (b) The optimal control rate is set at the minimum value, and (c) The outcome of this program for . The novelty of the algorithm used in this computation is that need not to be zero at for successful replacement of WU mosquitoes. The parameter values used are in Table 1 with WU mosquitoes at the WU-only equilibrium point at .The grey line in (c) is the release program period.

For this case, we can only have a mix population and we will want more WI mosquitoes. Hence, equation (2.4) becomes

For the constant release rate strategy. Since we want more WI mosquitoes, , and . Thus, . Again, we set and using values in Table 1 with gives . We varied δ and implement this release program for two years. From Fig. 9, after stoping the release program, it does not take much time for the WU mosquitoes to regain their dominance. If this strategy is to be adopted it will require continous application until all arboviral infections are eliminated.

Fig. 9

The δ whenand. (a) , (b) and (c) . The other parameter values are as in Table 1 and the grey area showing the period of release.

For the optimal variable release rate strategy, we want at a particular time that the WI mosquitoes will be greater than the WU mosquitoes and maintains it. From Fig. 9, it shows that this will be an ongoing process and that we can make WI mosquitoes more abundant as we want (Fig. 9c) depending on the allowable budget. Hence, since the release program will be ongoing, the benefit of this program is well captured when we incorporate the human infection dynamics to see the level of reduction in arboviral infections as a result of this scheme. This is out of scope of this manuscript.

Discussion and conclusion

In this study, we developed and analysed a complex dynamical system of a two-type-mosquito population in the presence of imperfect maternal transmission and loss of Wolbachia infection in order to determine the necessary and sufficient conditions for the propagation of Wolbachia infection in an Ae. aegypti population. We derived the invasive reproductive numbers with or without the adult WI mosquitoes losing their Wolbachia infection and established the conditions for local and global stability of the equilibrium points. We further adopted our models to determine the optimal release program that will ensure that WI mosquitoes replace or are become more abundant than the WU ones. Our analyses showed that mosquitoes with Wolbachia infection can dominate, co-exist or die out depending on whether they are fitter than mosquitoes without Wolbachia infection. The results showed clearly which factors and conditions are necessary and sufficient for WI mosquitoes to persist.

From our modelling, and consistent with other modelling works (Ferguson et al., 2015; Ndii et al., 2012; Xue et al., 2017), if WI mosquitoes are introduced in small numbers, the WI mosquitoes will not disrupt and outbreed WU ones. However, continuous introduction of WI mosquitoes for a particular period of time (O’Neill et al., 2018) will rescale the invasive reproduction number and increase it to above one where the WI-only mosquito equilibrium point is globally stable. A model in (Xue et al., 2017) derived the conditions for WI mosquitoes to invade the Ae. aegypti population and showed that the Wolbachia can still spread despite the fact that the number of Wolbachia-infected offspring due to WI adult mosquitoes in the next generation was less than one. By assuming an equal sex ratio between male and female Ae. aegypti mosquitoes, imperfect maternal transmission and loss of Wolbachia infection, we were able to find this threshold explicitly. This threshold does not guarantee replacement of uninfected mosquitoes as there are two possible steady-state solutions with one stable and the other unstable. However, Wolbachia-infected mosquitoes can dominate if the defined invasive reproductive number is greater than one (this is not possible realistically except external measure such as deliberate introduction of WI mosquitoes is implemented) and other conditions stated in our results are satisfied. Then it is possible to completely replace the WU mosquitoes or have them in less proportion depending on the Wolbachia infection dynamics in the Ae. aegypti population as shown by the optimal control problem.

Our study has some limitations that may affect our conclusions. One, we assumed that the ratio of male to female Ae. aegypti mosquitoes is the same. This has been shown under a laboratory study and may not be necessarily true in a real-life situation (Arrivillaga, 2004). Whatever the ratio of males to females in the number of eggs laid by either WU or infected female Ae. aegypti mosquitoes, the main factors that determines Wolbachia take-over is the proportion of WI eggs in the next generation and the per capita death rate of WI mosquitoes. Two, most of the parameters in the associated invasive reproductive number are seasonally dependent (Ndii, Allingham, Hickson, & Glass, 2016; Yang, Macoris, Galvani, Andrighetti, & Wanderley, 2009). The dependency of key parameters on temperature is likely to affect WU and infected mosquitoes in a similar way. This means that we are likely to have Ae. aegypti populations with any of three possibilities: without Wolbachia; with Wolbachia; and co-existence. Adverse conditions for the female Ae. aegypti mosquitoes are likely to reduce the ability to reproduce and fertilize their eggs, and this is likely to push the population towards the no-mosquito equilibrium point rather than changing the proportion of the Ae. aegypti population with Wolbachia infection. Lastly, we have mimiced the transmission dynamics of the wMel strain of Wolbachia in this work. Other strains such as wMelPop and wAlbB have similar dynamics to the wMel strain but with some variabilty in their reproductive advantage that affects the rate of introduction to ensure replacement (Xue, Fang, & Hyman, 2018). However, the wAu strain does not have the advantage of CI but comes with a high virus transmision blocking potential (Ant et al., 2018). For the wAu strain, our modelling is not applicable as the lack of CI implies that this strain has to be combined with another Wolbachia strain that has the advantage of CI and superinfection (Ant et al., 2018). Further, the experimental modelling work by Ferguson et al. (Ferguson et al., 2015) showed different transmission settings that different strains of Wolbachia can be adopted to reduce dengue infections. However, with the evident of loss of cytoplasmic incompatibility under field conditions (Ross et al., 2019) this needs to be revisited and our modelling work can be a template.

In general, our modelling work in this study complements existing works (Campo-Duarte et al., 2018; Ndii et al., 2012; Qu et al., 2018; Rafikov et al., 2019; Xue et al., 2017) and bridges the gap between alternative ways that WU mosquitoes may have advantages over the WI infected ones. Hence, controlling dengue epidemics and other arboviral infections with Wolbachia is promising but implementing the strategy comes at a cost that requires careful evaluation. We have shown the potential outcomes of implementing such a strategy and the key parameters that could be targeted to achieve the desired objectives. Another question is what level of other vector control method is needed to aid Wolbachia propagation so progress made will not be lost. If Wolbachia has a strong and sustained effect in the Ae. aegypti population, it remains to be seen whether the Ae. albopictus will take over as a key vector agent for viral transmission or whether climate change can negate all gains from such an introduction. These are questions for future research.

Declaration of competing interest

The authors declare no conflict of interest.