COMPLEX WOLBACHIA INFECTION DYNAMICS IN MOSQUITOES WITH IMPERFECT MATERNAL TRANSMISSION

. Dengue, malaria, and Zika are dangerous diseases primarily transmitted by Aedes aegypti , Aedes albopictus , and Anopheles stephensi . In the last few years, a new disease control method, besides pesticide spraying to kill mosquitoes, has been developed by releasing mosquitoes carrying bacterium Wolbachia into the natural areas to infect the wild population of mosquitoes and block disease transmission. The bacterium is transmitted by infected mothers and the maternal transmission was assumed to be perfect in virtu-ally all previous models. However, recent experiments on Aedes aegypti and Anopheles stephensi showed that the transmission can be imperfect. In this work, we develop a model to describe how the imperfect maternal transmission aﬀects the dynamics of Wolbachia spread. We establish two useful identities and employ them to ﬁnd suﬃcient and necessary conditions under which the system exhibits monomorphic, bistable, and polymorphic dynamics. These analytical results may help ﬁnd a plausible explanation for the recent obser- vation that the Wolbachia strain w MelPop failed to establish in the natural populations in Australia and Vietnam.


1.
Introduction. Aedes aegypti, Aedes albopictus, and Anopheles stephensi, are the primary vectors of some life-threatening mosquito-borne diseases such as dengue, malaria, and Zika. Since there are no effective vaccines available for these diseases, current controls mostly rely on environmental management by destructing vector breeding sites physically or chemically. Unfortunately, the chemical method of spraying pesticides causes severe environmental pollution and induces mosquito's insecticide resistance.
The study on the interaction of Wolbachia and mosquitoes has flourished due to the groundbreaking work [1,2,37,38,39,40], which paved an avenue for a biologically safe method of mosquito-borne disease control [5,10,17,18]. As expected, data indicates that uninfected ova produced by the wRi infected females are susceptible to CI. However, the wRu in D. simulans and the wMel in D. melanogaster cause weak or non-CI. Hence, among-female variation of Wolbachia deposition in ova seems to be generally associated with among-male variation in the effect of Wolbachia on sperm in Drosophila. For systematic studies on Wolbachia spreading in D. simulans, we refer the readers to Hoffmann and Turelli's work since 1990s [15,16,30,31,32,33,34] in which both Wolbachia infidelity and partial CI were identified.
It was predominantly believed that the maternal transmission of Wolbachia in mosquitoes is perfect. However, a very recent study [28] demonstrates that high temperatures (26-37 • C) can induce both imperfect maternal transmission and incomplete CI of Wolbachia in Aedes aegypti, the primary vector of dengue. Imperfect maternal transmission was also observed in a recent study on the spread of Wolbachia strain LB1 in Anopheles stephensi, along with a confirmed complete CI in [2]. The findings in [2] show that self-crossing LB1 produces significantly lower hatching (∼52.4%) than the outcrossing of LB1 females with uninfected males (∼80%) and uninfected self-crossing (∼90%), see Fig.1(B) in [2]. With the potential of a Wolbachia-based intervention for vector control on mosquito-borne diseases, we develop a model with complete CI and imperfect maternal transmission, and present a detailed mathematical analysis on the dynamics of Wolbachia spread. The main purpose of this manuscript is to dissect the effect of imperfect maternal transmission on Wolbachia spreading, which has not been studied in the existing literature.
To proceed, we divide the mosquitoes into four classes: uninfected females U F , infected females I F , uninfected males U M , and infected males I M . Assume that U F : U M = I F : I M and hence the population can be formally considered hermaphroditic. Let I and U denote the total numbers of infected, and uninfected mosquitoes, respectively, i.e., I = I F + I M and U = U F + U M . Let µ ∈ (0, 1) be the proportion of uninfected ova produced by an infected female. Following [22,41], we let b U denote the average number of offspring produced by one uninfected female. Let δ U denote the constant death rate for uninfected mosquitoes. For generality, we introduce b I and δ I as the corresponding parameters for infected mosquitoes as Wolbachia infection usually modifies the fitness and the longevity [25,35] of the host. CI arises when an uninfected ova, born from infected females or uninfected females, is fertilized by the sperm from an infected male. Then, with complete CI between I M and U F [1,2,35,39] and random mating, no uninfected adult progeny will be observed from the crossing of I F × I M . Uninfected adult progeny still arise from the crossing I F × U M due to imperfect transmission. With random mating, the chance of a female crossing with U M is U M /(I M + U M ) = U/(I + U ), and the chance of a female crossing with I M is I M /(I M + U M ) = I/(I + U ). The latter one is also the probability of CI-occuring. Our model takes the form In (2), the first term accounts for the birth of uninfected progeny from infected mothers, and the second term accounts for the birth from uninfected mothers, who have a probability I/(I + U ) to mate with infected males, and a contribution to the birth reduced by a factor 1 − I/(I + U ) = U/(I + U ) due to complete CI.
If we apply the rescaling and rewrite d · /ds as d · /dt, then we can transform (1)-(2) into In [41], we studied the case of perfect maternal transmission (i.e., µ = 0), and included the average waiting time τ from parent mating to the emergence of reproductive progeny as a time delay in the model. We proved that when δ = 1, i.e., the infection does not alter the mean life span, Wolbachia can spread into the whole population as long as the infection frequency stays strictly above a threshold value for a period no less than the reproductive period τ . Interestingly, we also found that such a threshold value of infection frequency cannot be well defined for δ = 1. For the reduced system without maturation delay, we obtained sharp estimates on a threshold curve determined by the population sizes of both infected and uninfected mosquitoes. This is different from a differential equation for the fraction of Wolbachia-infected mosquitoes [29], with the assumptions of constant population size and perfect maternal transmission of Wolbachia, where the threshold infection frequency is unique.
The remaining part of this paper is divided into three sections. In Section 2, we prove invariance and boundedness of solutions to system (4)- (5), and present two useful identities together with additional preliminary results. These help us study the global dynamics in Section 3 and prove that the system can exhibit monomorphic, bistable, and polymorphic dynamics, and give sufficient and necessary conditions for each case. Imperfect maternal transmission could lead to infinitely many polymorphic states, see Theorem 3.1. Also, our discussions in Section 4 highlights some challenges associated with the release of w MelPop.
If none of (C1)-(C3) holds, then system (4)-(5) does not admit any interior equilibrium point. For convenience of later discussion, we classify the remaining cases as the following 6 cases. 2.3. The Jacobian matrix and two auxiliary functions. As f (x, y) > 0 and g(x, y) > 0 when both x > 0 and y > 0 are sufficiently small, E 0 (0, 0) is a repeller.
To determine the stability of other equilibria, we first compute the determinant and the trace of the Jacobian matrix [13] of system (4)- (5): At E 1 (β/δ, 0) and E 2 (0, 1), direct calculations give and at E * (x * , y * ) by using y * = kx * we find Proof. Write J(E * ) = (J ij ) 2×2 . By using the definition of x * we obtain By the definition of k in (11), we have the identity and hence J 11 can be further reduced to The other three entries J 12 , J 21 andJ 22 can be simplified similarly, and J(E * ) is reduced to WOLBACHIA SPREAD WITH IMPERFECT MATERNAL TRANSMISSION 529 (i) For notation simplicity, instead of calculating det J(E * ), we computē and then Note that (11) is equivalent to each of the following two identities which help us simplifyJ to the following (ii) Again, by using of (15), trJ * can be simplified to This completes the proof.
Then the derivative of z(t) satisfies Proof. By rewriting z(t) = (1/δ) ln x(t) − ln y(t) and taking derivatives, we find This completes the proof.
For fixed constant c > 0, we denote by L c the ray y = cx, x > 0. Also for a solution γ(t) = (x(t), y(t)) of (4)- (5) Proof. By interchanging cx with y when necessary, we obtain from (4)- (5) that By applying the relation (18) and the proof is completed.
3. Stability of equilibria. E 0 (0, 0) is always a repeller since f (x, y) > 0 and g(x, y) > 0 when both x > 0 and y > 0 are sufficiently small. In this section, we shall analyze the stability of equilibria and the global dynamics of system (4)- (5). The primary tool in our proof is (13) and (14) in Lemma 2.3, differential identities (18) in Lemma 2.4,and (19) in Lemma 2.5.
3.1. The degenerate case (C1). Assume that (C1) holds. Then both E 1 and E 2 are degenerate since J(E 1 ) and J(E 2 ) have zero as their eigenvalues. Furthermore, system (4)-(5) is reduced to and the two nullclines defined in (7) and (8) are identical. All points on this coincident nullcline in R 2 +0 , including E 1 and E 2 on the boundary, are equilibrium points.
Proof. The relation (21) is easily derived by solving (20) with the initial condition. The only non-trivial part in the proof of the remaining conclusions is to show that (21) intersects the nullcline Γ y (or Γ x as they are identical) exactly once at it holds that x + y < βµx + y = (x + y) 2 < βµx + βµy, 1 < x + y < βµ.
Now if βµ ≤ 2, by taking derivative with respect to x on both sides of (8), we have Hence y decreases in x along Γ y in R 2 +0 , and the increasing curve (21) intersects it exactly once for x + ∈ (0, βµ), see Figure 1-(A).

3.2.
Polymorphism and bistability. The following theorem shows that Condition (C2) leads to a polymorphic scenario.
Proof. Both E 1 and E 2 are saddle points since their Jacobian matrix have one positive and one negative eigenvalue. By (13), det J(E * ) > 0. Meanwhile, we have and hence, Thus E * is locally asymptotically stable. To prove its global stability, it suffices to show that there is no closed orbits. Take in R 2 +0 . The divergence of the vector field (Bf, Bg) satisfies ∇ · (Bf, Bg) which declares nonexistence of closed orbits in R 2 +0 by the Dulac's criterion [13]. The following theorem indicates that (C3) corresponds to a bistable scenario. Theorem 3.3. Assume that (C3) holds. Then the unique interior equilibrium point E * (x * , y * ) is a saddle point, and E 1 and E 2 are locally asymptotically stable. Moreover, denoting by y = h(x) the stable manifold of E * , h(x) is continuous and increasing, and is the separatrix that defines the basins of attraction of E 1 and E 2 in R 2 +0 , and it lies between h 0 (x) and h 1 (x), where Proof. Condition µδ < 1 (resp. β(1 − µ) < δ) in (C3) implies that J(E 1 ) (resp. J(E 2 )) has two negative eigenvalues, and hence both E 1 and E 2 are locally asymptotically stable under (C3). Moreover, det J(E * ) < 0 by (13), and hence E * is a saddle point. Noticing that along y = cx, (18) implies that We have By (19), we have There are three cases to consider.
(i) δ = 1. Then l k (t)| L k ≡ 0, implying that L k : y = kx consists of solution curves of system (4)- (5), which also serves as the stable manifold of E * .
(ii) δ > 1. Consider the domain We claim that D 1 is negatively invariant with respect to the dynamics of (4)- (5). In fact, on its upper boundary where y = h 1 (x), (27) implies that z (t) < 0. On its lower boundary y = h 0 (x), It follows that solutions of system (4)-(5) initiated within D 1 other than the origin and E * will leave this region in finite time and stay outside thereafter. Similarly, consider Then we have l k (t)| L k < 0 on its upper boundary y = h 0 (x), and z (t) > 0 on the lower boundary y = h 1 (x). Consequently, the stable manifold of E * must locate inside D 1 and D 2 and (2) is verified.
(iii) δ < 1. By using the same argument as above, the negative invariance of and D 4 = (x, y) : x > x * , h 0 (x) < y < h 1 (x) could be proved by dissecting the sign of l k (t)| L k and z (t) on the boundaries.

Monomorphism.
Theorem 3.4. Assume that one of the Conditions (C4), (C5) and (C6) holds. Then E 1 is globally asymptotically stable: for each point (x, y) ∈ R 2 +0 , ω(x, y) = E 1 . Proof. From (18), z (t) > 0 always holds for these 3 cases, showing that (4)- (5) does not have any non-trivival periodic solution. We claim that E 2 could not be the ω-limit set of the trajectory (x(t), y(t)) initiated from R 2 +0 . For cases (C4) and (C6), it is trivial since E 2 is a hyperbolic saddle point with the stable manifold on the y-axis. However, E 2 becomes non-hyperbolic for case (C5), whose local stability cannot be determined by the Jacobian matrix approach. In this case, since z (t) > 0, it holds that z(t) > z(0) for all t > 0, and so E 2 cannot be an ω-limit point. Hence, ω(x, y) = E 1 by the classical Poincaré-Bendixson theorem [13].
Theorem 3.5. Assume that one of the Conditions (C7), (C8) and (C9) holds. Then E 2 is globally asymptotically stable in the sense that for each point (x, y) ∈ R 2 +0 , ω(x, y) = E 2 . Proof. It is easily seen that z (t) < 0 for all cases, which again implies that (4)- (5) does not have non-trivial periodic solution. The equilibrium point E 1 could not be the ω-limit set of any trajectory (x(t), y(t)) in R 2 +0 for cases (C7) and (C8) since it is a hyperbolic saddle point with the stable manifold on the x-axis. For case (C9), since z (t) < 0, it holds that z(t) < z(0) for all t > 0, and so E 1 cannot be an ω-limit point. Hence, ω(x, y) = E 2 by the classical Poincaré-Bendixson theorem [13].

Discussion.
4.1. Polymorphism, bistability and monomorphism. The asymptotic behavior of the Wolbachia spreading dynamics in a mixed population of infected and uninfected mosquitoes depends strongly on the life parameters modified by the Wolbachia infection. We classify the system parameters into Cases (C1) to (C9) and identify their correspondence to three different scenarios: Polymorphism: Theorems 3.1 and 3.2 show that the mixed population remains polymorphic for which both infected and uninfected mosquitoes coexist when (C1) or (C2) holds. If the degenerate case (C1) occurs, then the ultimate coexistent state may change with the initial infection level. When (C2) occurs, the asymptotic population size is independent of the initial value and the infection frequency at the steady-state is x * /(x * + y * ).
Bistability: For case (C3), Theorem 3.3 offers a sharp threshold given by the separatrix y = h(x) on the initial values for successful Wolbachia invasion. If Wolbachia does not alter the longevity of the mosquito, i.e., δ = 1, then the stable manifold of E * coincides with h 0 (x) ≡ h 1 (x), see  Monomorphism: Theorem 3.4 implies that under either (C4) or (C5) or (C6), for any initial positive release, the Wolbachia infected mosquitoes will eventually replace the wild ones. Theorem 3.5 shows that under either (C7) or (C8) or (C9), Wolbachia infection will fail to establish, regardless of the number of initial release.
With perfect maternal transmission rate, i.e., µ 1 = µ 2 = 0, system (4)-(5) is reduced to for which only (C3), (C4) and (C5) can occur. When β/δ ≥ 1, the environment is more favorable (or at least equally favorable) for infected mosquitoes, and we say that Wolbachia infection has a fitness benefit. In this case, system (28)- (29) only admits the Wolbachia free equilibrium point E 2 and the Wolbachia replacement equilibrium point E 1 . Theorem 3.4 for Cases (C4) or (C5) shows that successful Wolbachia spreading is ensured from any initial release above 0 for the fitness benefit case. When β < δ, we say that the infection has a fitness cost, and besides E 1 and E 2 , there is a unique interior equilibrium point of system (28)- (29) which is a saddle point. Theorem 3.3 for Case (C3) displays bistable dynamics for the fitness cost case. These results are in line with the results in [41] which also include the corresponding results in [41] as special cases.
Comparing to perfect maternal transmission, the imperfect maternal transmission generates richer dynamics. To see this, we fix β and δ and consider separately the fitness cost case (i): β ≤ δ, and the benefit case (ii): β > δ. We explore the impact of µ > 0 in terms of its relation with two critical values The detailed outcomes in different cases are shown in Figure 3. For example, in the fitness cost case (i), small µ corresponds to the bistable case (C3). However, increasing µ leads to (C9) and (C7), for which all solutions are attracted to E 2 , and Wolbachia will be eventually wiped out, no matter how many infected mosquitoes are initially released. Thus, we see a transition from the bistable state to the monomorphic state E 2 as µ increases. In the fitness benefit case (ii) with ν 1 < ν 2 , small µ corresponds to (C4) and (C6), for which all solutions converge to E 1 : Wolbachia infected mosquitoes will eventually replace the wild ones, as long as µ < ν 1 . But when µ is large, (C4) is replaced by (C2), (C8) and (C7) consecutively, which again implies the failure of the release because all solutions converge to E 2 . Thus, there is also a transition of convergences to equilibria when µ > 0 is increased: transition from convergence to the replacement equilibrium E 1 to convergence to E * (the unique positive equilibrium, giving a unique polymorphic outcome), and finally to the failure equilibrium E 2 . This is in contrast to the transition in (i) which is a transition from the bistable state to the monomorphic state E 2 .

4.2.
Success or failure: Review on some recent releases. We consider two Wolbachia strains, the benign w Mel and the virulent w MelPop, of Aedes aegypti in Australia. It was found that w Mel caused only 10% longevity reduction and no significant fecundity cost [35], but w MelPop decreased the longevity of infected female mosquitoes by ∼50% and resulted in a fecundity cost of ∼56% [25]. By using the laboratory data in [25,35], we estimated the birth rates b I and b U from the oviposition rate, egg hatching rate and survival rate of larvae to adults, as well as the death rates δ I and δ U from the half-life of adults [41]. We found that b U = 0.3976 and δ U = 8.5034 × 10 −6 , and for the benign Wolbachia strain  In 2011, the Wolbachia strain w Mel was introduced into Aedes aegypti populations in some northern Australian field sites [17] and had almost completely replaced wild mosquitoes for more than 3 years [10,18]. Since September 2014, the w Mel mosquitoes were released in two neighborhoods in Brazil [5]: Tubiacanga in Rio de Janeiro (2014-2015) and Jurujuba in Niteroi (2015) to reduce dengue transmission. Successful or ongoing field trials also include the w Mel release in January 2014 in Indonesia, and in May 2015 in Colombia. Achievements of these field trials have been so encouraging that the strategy is becoming a worldwide plan to combat Zika and other mosquito-borne viruses.
In contrast to the effective and sustainable invasion of the w Mel release in Australia, the 2013 w MelPop trial in central Vietnam (Tri Nguyen, Hon Mieu Island) [26] was a failure: w MelPop could be established transiently, but not persistently, in any of the two sites. The failure was largely attributed to the virulent property of the w MelPop that causes the infected mosquitoes to have a shorter lifespan (δ = 1.6667 > 1) and produce fewer offspring than the wild ones (β = 0.5418 < 1).
We can claim that maternal transmission leakage even makes a transient establishment of Wolbachia in the wild mosquito population a big challenge. Consider again the w MelPop strain as an example. Given the uninfected mosquito's density y 0 and the leakage rate µ, if the maternal transmission leakage occurs, then (4)-(5) determines a unique Wolbachia release threshold x 0 = x 0 (y 0 , µ). We plot the ratio x 0 /y 0 against the leakage rate in the left panel in Figure 4. It displays a steep elevation of the threshold ratio x 0 /y 0 for a successful Wolbachia invasion as the leakage µ increases from 0.3 to 0.5, and reaches as high as 16 when µ = 0.5. Furthermore, should the leakage is accompanied by incomplete CI as shown in [28], the scenario is even more discouraging. It was shown in [28] that, among 1848 eggs from the crossing of uninfected females and w MelPop-infected males, 301 (16.31%) hatched. This leads to a reduction of the CI intensity s h from s h = 1 in (4)-(5) to s h = 0.8369. With the incorporation of incomplete CI, the model is modified to Again, we can determine the unique threshold value x 0 = x 0 (y 0 , µ, s h ). In this case, a small increase on the leakage rate could bring a catastrophic failure of Wolbachia release. As shown in the right panel in Figure 4, the threshold ratio is as large as 7 when µ = 0.15. Further increase of the leakage rate from 0.15 to 0.25 results in an augment of 6-fold on the threshold ratio. If µ = 0.3, we found that even if the ratio is as large as 200, Wolbachia will be wiped out. Our example highlights some challenges associated with releasing w MelPop, which is in line with the results in [26,28]. When the leakage rate is 0.15, the incomplete CI mechanism leads the threshold ratio increase from 1.95 to 7. Worse than that, when the leakage rate is 0.3, the ratio as high as 200 is still insufficient to ensure successful Wolbachia invasion.