Characterization of Wolbachia enhancing domain in mosquitoes with imperfect maternal transmission

ABSTRACT A novel method to reduce the burden of dengue is to seed wild mosquitoes with Wolbachia-infected mosquitoes in dengue-endemic areas. Concerns in current mathematical models are to locate the Wolbachia introduction threshold. Our recent findings manifest that the threshold is highly dependent on the initial population size once Wolbachia infection alters the logistic control death rate of infected females. However, counting mosquitoes is beyond the realms of possibility. A plausible method is to monitor the infection frequency. We propose the concept of Wolbachia enhancing domain in which the infection frequency keeps increasing. A detailed description of the domain is presented. Our results suggest that both the initial population size and the infection frequency should be taken into account for optimal release strategies. Both Wolbachia fixation and extinction permit the oscillation of the infection frequency.


Introduction
The ability of Wolbachia to spread through cytoplasmic incompatibility (CI) [33,36] and grant mosquitoes resistant to dengue virus [2,21] has triggered the development of Wolbachia-based strategies for population replacement of Aedes mosquitoes, the primary vectors of dengue. Since 2009, small-scale trials have been conducted in dengueaffected communities, including Australia, Vietnam, Brazil and Colombia [10,11]. Before the implement of large-scale trials, one major concern is the introduction threshold of Wolbachia-infected individuals that must be surpassed for infection to spread and become established in the population.
Several key parameters are critical for determining the introduction threshold. (i) The intensity of CI s h , which is the probability of embryo death from the crossing of infected males with uninfected females; (ii) The maternal transmission leakage μ ∈ [0, 1], which is the percentage of uninfected progeny produced by an infected mother. (iii) The fitness cost/benefit, including changes in fecundity, hatching, pupation, eclosing and adult longevity, of infected females caused by Wolbachia infection. The classical finding of Turelli-Hoffman [29,31,32] on Wolbachia invasion in Drosophila simulans, based on discrete models parameterized using laboratory and field data, showed that the dynamics of Wolbachia spread can be determined completely by the frequency of infected among all individuals p(t). They showed that there is a unique threshold of the infection frequency p * ∈ (0, 1): If the initial infection frequency p(0) > p * , then p(t) increases for all t > 0 and p(t) → 1 as t → ∞ for Wolbachia fixation; if p(0) < p * , then p(t) decreases for all t < 0 with p(t) → 0 as t → ∞ to drive Wolbachia extinction. Such a simple correspondence between the infection frequency and Wolbachia spreading dynamics has been justified repeatedly for laboratory experiments of mosquitoes with separated generations [6,10,36].
However, for Wobachia spread in a natural setting where mosquitoes are present in overlapping generations, the situation is much more complicated. Various mathematical models have been developed to characterize the thresholds, including models of ordinary differential equations [16,22,[38][39][40], delay differential equations [15,39], impulsive differential equations [37], stochastic equations [12], and reaction-diffusion equations [13,14,30]. In most of the above theoretical studies on characterization of the introduction thresholds, it was assumed that the maternal transmission is perfect (μ = 0) and CI is complete (s h = 1), which are supported by Wolbachia strain WB1 in Aedes aegypti [36], wAlbA and wAlbB in Aedes albopictus [35]. Although imperfect maternal transmission (μ > 0) has been repeatedly documented with respect to the wRi, wRu infection in D. simulans [5,17,18,32] and the wMel in Drosophila melanogaster [9], It was predominantly believed that the maternal transmission of Wolbachia in mosquitoes is perfect. However, a very recent study [26] demonstrates that high temperatures (26-37 • C) can induce both imperfect maternal transmission of Wolbachia in A. 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 [3].
Motivated by these findings, we developed models of ordinary differential equations in [40,41] to describe Wolbachia spreading dynamics, aiming to understand the effect of imperfect maternal transmission. To the end, let I(t) and U(t) be respectively the population size of infected and uninfected adults at time t, both equally distributed in sex. Let b I (or b U ) be the natural birth rate of infected (or uninfected) individuals. As infected offsprings are born only from infected mothers, by taking the maternal transmission leakage into account, the birth of infected offsprings reads as b I (1 − μ)I/2. The birth function of uninfected offsprings is in which the first term accounts for the leakage of infected mothers, and the second term is from uninfected mothers, with the CI intensity s h = 1 for matings between uninfected females and infected males. On the loss of mosquitoes, we follow the conventional approach by replacing the natural death of adults with a logistic-like density dependence term, which increases in I (or U) and the sum I+U. Let δ I (or δ U ) be the logistic control coefficient for infected (or uninfected) mosquitoes. The loss of reproduction due to density dependent death is modelled by δ I I(I + U) (or δ U U(I + U)) for infected (or uninfected) mosquitoes. In summary, our model takes the form Upon the rescaling and rewriting d · /ds as d · /dt, system (1) and (2) is reduced into where β (or δ) is the birth (or death) rate constant of infected females relative to uninfected females, and x, y are respectively the rescaled population size of infected and uninfected individuals.
Our recent studies [40,41] showed that system (4) and (5) exhibits monomorphic, bistable, and polymorphic dynamics. A detailed description of the threshold curve is offered in terms of β, μ and δ, as we have done for perfect maternal transmission case in [39]. The results suggest that the largest maternal transmission leakage rate supporting possible Wolbachia spreading does not necessarily increase with the fitness of infected mosquitoes. By exploring the analytical property of the threshold curve, we find that with the presence of imperfect maternal transmission rate, Wolbachia in a completely infected population could be wiped out ultimately if the initial population size is small. All these findings point to the fact that even when μ = 0 which corresponds to perfect maternal transmission, there is no more single threshold infection frequency unless δ = 1, and the classical result of Turelli-Hoffman needs to be interpreted differently. Indeed, we showed that there is a threshold curve, depending on the initial population sizes of both infected and uninfected mosquitoes when δ = 1. Although we offered a complete classification of the Wolbachia infection dynamics by the proof of a unique threshold curve, the separatrix of an unstable saddle point, the implication of results requires a knowledge of the initial population sizes x(0) and y(0). This is, practically, a formidable task in large areas with complex landscape structures, in particular, in residential areas. On the other hand, the infection frequency may be estimated by undertaking PCR assay on samples or progeny tests on subsamples from various locations [32]. To make our theoretical results be more helpful in surveillance and adjustment of release strategies, here, we proposed the concept of Enhancing domain in terms of the growth rate of the Wolbachia infection frequency which consists of all points (x, y) at which the infection frequency increases. We call E I the Wolbachia infection enhancing domain in which infected mosquitoes is more favoured comparing to uninfected ones. The aim of this paper is to characterize E I and elucidate the relation or most importantly, the difference between the enhancing domain and the attraction domain of Wolbachia fixation equilibrium, denoted by A I . Our results show that, when δ = 1 or μ = 0, the relative complement of E I in A I is not empty, which implies that the infection frequency could go to Wolbachia fixation even though it does not always increase. Meanwhile, the relative complement of A I in E I is not empty either, in which initially, the infection frequency increases, but eventually, Wolbachia will go to extinction. These deceptive phenomena should be taken into account when designing release strategies to guarantee the success.

Preliminaries
When the maternal transmission of Wolbachia is perfect and the infection does not affect the lifespan of infected females, i.e. μ = 0, and δ = 1. From (4) and (5), we have dp dt = d dt If β > 1, then dp/dt > 0 always holds, and hence E I = E o . If β < 1, then the infection frequency presents an Allee effect: p(t) goes to 1 for Wolbachia fixation when p(0) > p * := 1 − β, and to 0 for Wolbachia extinction when p(0) < p * . In this case, Noticing that system (4)and (5) admits no interior equilibrium when β > 1, and a unique interior equilibrium E * (x * , y * ) = (β − β 2 , β 2 ) when β < 1, combining the characterization of A I in [39] (Theorem 2.1), we can conclude that the enhancing domain is identical to the attracting basin of Wolbachia fixation equilibrium (β, 0) when β < 1. In summary, we have However, things are quite different for μ = 0 or δ = 1. Since and g(x, y) is independent of δ, we have the following simple monotone properties In order to identify the domain E I in a systematic way, we introduce, for fixed constants c > 0, l c (t) := cx(t) − y(t). We denote by L c the ray y = cx, x > 0. If γ (t) = (x(t), y(t)) ∈ L c for some t ≥ 0, then l c (t) = 0. If it is additionally known that l c (t) > 0, then l c (s) > 0 for all s sufficiently close to t with s > t, implying that A similar argument shows that p(s) < p(t) for all s sufficiently close to t and s < t. As a result, For (x, y) ∈ E o , the line connecting it with the origin cuts the x−nullcline exactly once. We denote by (x s , y s ) the unique cutting point, and call it the shadow point of (x, y) on x . If (x, y) is itself on the nullcline, then it coincides with its shadow. In general, (x, y) and its shadow point are colinear with the origin, and (4) and (5) All points in L c share the same shadow point on the x−nullcline c given by The following two results will be frequently used to characterize the enhancing domain, which have been proved in [40].
When δ = 1, z(t) and p(t) do not relate in any obvious fashion and (11) is no longer powerful for the study of E I . Lemma 2.1 provides a primary tool for us to study E I . For instance, assume that δ > 1 and the system has two interior equilibria E * 1 and E * 2 . We see from (9) that l c (t)| L c > 0 in the triangle connecting E * 1 , E * 2 , and the origin, including the three edges but not the vertices. Hence the triangle is contained in E I . It is also clear from (9) that E I must contain some neighbouring areas of the triangle, and so E I is not completely contained within the whole sector bounded by the two rays L c passing through E * 1 and E * 2 . Our study shows that when δ > 1, E I is contained inside a larger sector, whose boundary, to our surprise, consists of the rays connecting the origin with each of the interior equilibrium points corresponding to the accompanied system with δ = 1. Let δ > 1, β > 0, and μ ∈ (0, 1). Then we have

Theorem 3.2:
where (x c , y c ) is the shadow point of (x, y) on the x−nullcline x + y = κ. Then (3) If β(1 − μ) < 1, β < 2, and μ < β/4, then system (4) and (5) with δ = 1 has two interior equilibria E 1 (x 1 ,ȳ 1 ) and E 2 (x 2 ,ȳ 2 ) withȳ 2 <ȳ 1 . Letc 1 =ȳ 1 /x 1 andc 2 =ȳ 2 /x 2 . Then Proof: From Lemma 2.1 we see that E I consists of the segments in L c such that Let c > 0 be fixed, and (x, y) ∈ L c . If (x, y) ∈ E I , then Recall from (8) that x c = κ/(1 + c) and y c = cx c , we find This indicates that the segment of L c with 0 < x < x c belongs to E I , provided that x c > 0. If x c ≤ 0, then the whole ray L c stays outside of E I . As y c = cx c , x c > 0 if and only if By recalling the definition of G in (10), we obtain its equivalent conditions as As κδ + βμ = β(1 − μ) + βμ = β, we find Now we denote by (x c ,ỹ c ) the shadow point of (x, y) ∈ L c on x + y = β (1 − μ), the x−nullcline for system (4) and (5) with δ = 1. Thenx c = δx c ,ỹ c = δy c . With δ = 1, the G function introduced in (12) takes the form and therefore We finally derive from (17) the following useful relation (1) Since we asserted in Theorems 3.1 that E I is empty if μ ≥ β/4, or β(1 − μ) ≤ 1 and β ≥ 2 for the case δ = 1, it can be deduced from the monotone property (7) that E I is empty for δ > 1.
We remark that the conditions on β and μ are the same in each of Parts (1), (2), and (3) in Theorems 3.1 and 3.2. Under the conditions of Parts (2) or (3) in Theorem 3.2, E I is a bounded domain whose boundary is parameterized as x = x c , y = cx c , where c ≥c in Case (2) andc 2 ≤ c ≤c 1 in Case (3). The boundary passes through the origin and the interior equilibrium points should they exist. It has two tangent lines at the origin, which are given by the y− axis and Lc in Part (2), and Lc 1 and Lc 2 in Part (3). Let δ < 1, β > 0, and μ ∈ (0, 1). Then the infection enhancing domain E I is an unbounded set (20) and (x c , y c ) is the shadow point of (x, y) on the x−nullcline x + y = κ. Furthermore,

Theorem 3.3:
(2) If β(1 − μ) > 1, then the system defined by (4) and (5)  Proof: Let c > 0 be fixed, and (x, y) ∈ L c . By repeating the discussion in the proof of Theorem 3.2, from the beginning to (16), we see that (x, y) ∈ E I if and only if, as opposed to (16), .
Let (x c ,ỹ c ) be the shadow point of (x, y) ∈ L c on x + y = β (1 − μ), the x−nullcline for the system defined by (4) and (5) with δ = 1. By repeating the arguments from (16) to (19), we find that x c ≤ 0 ⇔ G 1 (ỹ c ) ≤ 0. The rest of the proof can be done by exactly the same arguments as in the proof of Theorem 3.2.
It is interesting to see from Theorem 3.3 that, when δ < 1, E I is not only non-empty but also unbounded, a sharp contrast to the case δ > 1 for which E I is either empty or a bounded set. In Case 1), the y−nullcline y stays above 1 , but y may intersect δ . If it happens that y stays above δ , then E I locates entirely above y . Indeed, the triangle formed by δ and the two axes does not belong to E I since x c > x c , and the region located between δ and y does not belong to E I either as x < 0 and y > 0 in the region. In Case 2), the boundary of E I is tangent to Lc at the origin. In case 3), the boundary is tangent to Lc 1 and Lc 2 at the origin.

Discussion
As a safe and novel strategy for controlling mosquito-borne diseases, releasing mosquitoes carrying Wolbachia or mosquitoes with lethal gene to suppress or replace the wild mosquito population has been implemented in areas where mosquito-borne diseases such as Zika, dengue and chikungunya are endemic, which has made the spread dynamics of Wolbachia a hot topic [1,4,19,20,24,25,27,28]. For Wolbachia fixation in wild areas, it is crucial to locate the introduction threshold, i.e. the minimum number of infected mosquitoes released to guarantee the spread and establishment of Wolbachia in the mosquito population. Our recent findings [39][40][41] shows that, with or without the maternal transmission leakage, there is no unique introduction threshold as stated in earlier systematic studies [9,10,29,31,32] once Wolbachia infection alters the lifespan of infected females, i.e. δ = 1. Instead, the introduction threshold depends on the initial population sizes x(0) and y(0), which is an arduous task in wide areas with complex landscape structure, especially in residential areas. An alternative method with extensive application is to test the infection frequency by PCR essays on samples, and then deduce the changes in the frequency of Wolbachia infected mosquitoes [7,8,10,11,23]. This makes the surveillance of temporal profile on the infection frequency as the basis of design and adjustment of release strategies, which is also the motivation of introducing the concept of 'enhancing domain' in this paper.
To see the difference between the enhancing domain E I and the attraction domain A I of Wolbachia fixation equilibrium, we take the case when κ < 1, β < 2, and μ < β/4 as an example due to the reason that Wolbachia infection usually causes fitness damage (κ < 1) to the infected mosquitoes. In this case, system (4) and (5) admits two interior equilibria E * 1 (x * 1 , y * 1 ) and E * 2 (x * 2 , y * 2 ). Let c * 1 = y * 1 /x * 1 and c * 2 = y * 2 /x * 2 . When δ = 1, the infection enhancing domain reads as from Theorem 3.1, which is Domain II in Figure 1(a). However, the attracting domain of E * [40]. The difference between E I and A I is the Domain I, in which initially, p(t) decreases, but eventually, it goes to the stable frequency at E * 2 .
When δ > 1, the enhancing domain is Domain II ∪ IV, which is bounded, while the attracting domain of E * 2 is Domain (I+II), see Figure 1(b). Hence, any solution to (4) and (5) initiated from Domain I approaches to E * 2 , but p(t) does not always increase. And any solution initiated from Domain IV approaches to Wolbachia-free equilibrium point E 1 (0, 1), but p(t) does not always decrease.
When δ < 1, the enhancing domain is Domain (II+III), which is unbounded, while the attracting domain of E * 2 is Domain (I+II+V), see Figure 1(c). Hence, any solution to (4) and (5) initiated from Domain I or Domain V approaches to E * 2 , but p(t) does not always increase. And any solution initiated from Domain IV approaches to Wolbachia-free equilibrium point E 1 (0, 1), but p(t) does not always decrease. The gap between E I and A I suggests that the design or adjustment of releasing strategies could not be based solely on the surveillance of the infection frequency. To see what does the conclusion imply biologically that may be instructive to the design of release strategies, we take the benign Wolbachia strain, wMel, of A. aegypti as an example since wMel infected mosquitoes are currently being deployed in several countries for the control of arboviruses [34]. The findings in [33] shows that wMel brings no significant fecundity cost to the host, and causes only 10% longevity reduction. Based on the laboratory data in [33], from the oviposition rate, egg hatching rate, the survival rate of larvae to adults we estimated the birth rates b I and b U as The decay rates δ I and δ U are estimated from the half-life of adults as δ I = 9.4482 × 10 −6 , δ U = 8.5034 × 10 −6 .
Hence, from the rescaling (3)  We take μ = 0.05 to account for the maternal leakage of infected mothers [26]. Under these parameters, system (4) and (5) admits two interior equilibria E * 1 is a saddle point, and E * 2 is an asymptotically stable polymorphic point [40], which denotes the coexistence of infected and uninfected mosquitoes with infection frequency 0.8045 0.8045 + 0.0505 ≈ 0.9409 := p stable .
To proceed, we take three special initial values: (0. 25 Compared to estimating the absolute population sizes of infected and uninfected mosquitoes, a more direct and plausible method is to monitor whether the more measurable quantity p(t) is going up or down. However, our observations imply that when