Wolbachia infection dynamics in mosquito population with the CI effect suffering by uninfected ova produced by infected females

Mosquito-borne diseases pose a great threat to humans' health. Wolbachia is a promising biological weapon to control the mosquito population, and is not harmful to humans' health, environment, and ecology. In this work, we present a stage-structure model to investigate the effective releasing strategies for Wolbachia-infected mosquitoes. Besides some key factors for Wolbachia infection, the CI effect suffering by uninfected ova produced by infected females, which is often neglected, is also incorporated. We analyze the conditions under which Wolbachia infection still can be established even if the basic reproduction number is less than unity. Numerical simulations manifest that the threshold value of infected mosquitoes required to be released at the beginning can be evaluated by the stable manifold of a saddle equilibrium, and low levels of MK effect, fitness costs, as well as high levels of CI effect and maternal inheritance all contribute to the establishment of Wolbachia-infection. Moreover, our results suggest that ignoring the CI effect suffering by uninfected ova produced by infected females may result in the overestimation of the threshold infection level for the Wolbachia invasion.


1.
Introduction. Mosquito-borne diseases pose a great threat to humans' health, including dengue fever, chikungunya, and Zika [12]. Wolbachia is a maternally transmitted endosymbiotic bacteria. About 65% of insect species and 28% of mosquito species have been surveyed [32]. It may sometimes be artificially introduced to assist in the control of vector-borne diseases. There are two main reasons why Wolbachia may reduce disease transmission. First, Wolbachia may shorten the lifespan of mosquitoes [27,29,30,42]. It will have a disproportionate effect on disease transmission through a reduction in the number of mosquitoes that live long enough to go through the cycle of pathogen infection, maturation, and transmission [27]. Second, Wolbachia affects the ability of mosquitoes to transmit diseases, probably through interactions with the pathogen mediated via the insect immune system [3,16].
Most common effects of Wolbachia infections include cytoplasmic incompatibility (CI), feminization of genetic males and male killing (MK) [23]. CI occurs when Wolbachia in the bodies of infected males somehow modify the sperms of their hosts, such that the arrest of embryonic development occurs unless the egg also carries the bacterium [44]. Therefore, uninfected females that mate with infected males rarely produce viable eggs, while infected females are not affected. This gives infected females an advantage and helps the bacteria to spread quickly through a mosquito population [6,22,31,47,49]. Besides, compared with uninfected ones, Wolbachia-infected mosquitoes may suffer fitness cost which is not inevitable as there may be some fitness advantages, which allows them to spread more easily and establish themselves in field trials [20,49]. Hence, in this work, we consider that the mortality rates of adult infected mosquitoes can be greater, equal, or less than the mortality rates of adult uninfected mosquitoes.
Many mathematical models have been proposed for the spread of Wolbachia ( [1,7,13,14,16,18,21,23,24,29,30,36,37,40,43,45,46,48,49,50,51], etc.) In [14], Farkas and Hinow proposed a single-sex model for Wolbachia infection with both age-structured and unstructured models to study the stability and equilibrium based on the assumption that Wolbachia infection leads to increased mortality or reduced birth rate. In [21], Keeling et al. proposed fixed sex ratio ODE models to study the competition and coexistence between multiple strains of Wolbachia in a well-mixed population. 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 ( [30]). Xue et al. incorporated 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 ( [46]). Zheng et al. showed that the successful population replacement would depend on the strains of Wolbachia used and require a careful release design ( [52]). Qu et al. extended the model in [46] to include female mosquitoes mating once ( [32]). Li and Liu proposed an impulsive mosquito population model with general birth and death rate functions and showed that strategies may be different due to different birth and death rate functions, the type of Wolbachia strains and the initial number of Wolbachia-infected mosquitoes ( [26]). Recently, some Wolbachia invasion model with optimal control have been proposed (See [4,5], etc.) Meanwhile, some literatures concentrated on the spatial spread of Wolbachia ( [19,40], etc.) In [9], the authors reviewed how mathematical models provide insights into the dynamics of the spread of Wolbachia, the potential impact of Wolbachia on dengue transmission, and discussed the remaining challenges in evaluation and development. Most recently, [53] showed that combining incompatible and sterile insect techniques (IIT-SIT) enables near elimination of field populations of the world's most invasive mosquito species, Aedes albopictus.
Most of these models consider either complete CI effect or incomplete CI effect suffering by eggs produced by uninfected females. Actually, if maternal transmission is imperfect, whether an embryo experiences CI is determined solely by the infection status of the male and the ovum, so that the uninfected ova produced by infected females are as susceptible to the sperms from infected males as are ova from uninfected females [10,33,38,15]. In other words, uninfected ova produced by infected females due to incomplete vertical transmission also suffer CI effect when they bind with the sperms from infected males. However, this critical factor is barely incorporated in Wolbachia invasion models to the best of our knowledge.
Our objective is to formulate a new stage-structured mosquito population model incorporating the key factors of Wolbachia infection including cytoplasmic incompatibility, male killing, incomplete maternal transmission, reduced reproductive output and different mortality rates, in particular, the CI effect suffering by uninfected ova produced by infected females, which is often neglected in other models, and to understand how Wolbachia infection can be established in a wild population of mosquitoes.
This work is organized as follows. In section 2, we propose a novel Wolbachiainfection considering the CI effect suffering by uninfected ova produced by infected females and other key effects of Wolbachia infection. In section 3, we analyze the existence and stability of equilibria, as well as backward bifurcation. In section 4, we present numerical analysis and sensitivity analysis to analyze the stability of coexistence equilibrium and the impact of some critical factors of Wolbachia infection. Besides, we propose some different releasing strategies to establish the Wolbachia infection effectively. In section 5, we provide some discussions and conclusions.

2.
Model. Mosquitoes undergo complete metamorphosis going through four distinct stages of development during a lifetime, including egg, larva, pupa, and adult [2]. We group the three aquatic stages of mosquitoes into one stage, i.e. aquatic stage. We denote the numbers of uninfected and infected immature mosquitoes by U A and W A , respectively. Adult mosquito population is subdivided into four groups, namely, Uninfected males (U M ), Infected males (W M ), Uninfected females (U F ), and Infected females (W F ). The state variables are listed in Table 1.
We assume that both infected and uninfected eggs have the same natural birth rate, θ and natural death rate, d a . The intraspecific competition rate between immature mosquitoes in aquatic stage is κ, while the intraspecific competition between adult mosquitoes is omitted. The offspring of Wolbachia-free and Wolbachia-infected mosquitoes have a proportion f u and f w of females, respectively. We denote ϕ as the development rate of immature mosquitoes, i.e., 1/ϕ is the time of aquatic stage for a mosquito. Mortality rates for uninfected/infected male/female individuals may be different, which are denoted as d m , d mw , d f , d f w , respectively. Wolbachia is mostly passed from infected females to their offspring. However, the transmission is usually not perfect in nature, occurring with a probability v ∈ (0, 1]. The effect of the CI mechanism results in zygotic death of eggs with a probability q ∈ [0, 1] when an infected male mates with an uninfected ovum. We introduce a parameter λ to show different scenarios of incomplete CI. The parameter λ = 0 means that the CI effect suffering by uninfected ova produced by infected females is neglected, while λ = q means that this effect has been considered. The parameter, α ∈ [0, 1) represents the reduced reproductive output of Wolbachia infected females, τ ∈ [0, 1) measures MK in the sense that it is the probability that a Wolbachia infected male larva dies during its development. A complete list of the parameters involved in model (1) is given in Table 2.
Based on the above assumptions, we formulate the following system to describe the dynamics of Wolbachia infection.
(1)  Table 2. Definition of the parameters 3. Dynamical analysis. In this section, we derive the basic reproduction number, R 0 , and analyze the existence and stability of the equilibria, as well as backward bifurcation.

State variables Uninfected Infected
3.1. Wolbachia-free equilibrium and R 0 . A straightforward calculation leads to a unique Wolbachia-free boundary equilibrium (WFE) .
It should be pointed out that R 0u can be biologically interpreted as the basic reproduction number of wild mosquito population, which can be explained by the means of biology. Following [41], the basic reproduction number can be computed using the next generation matrix approach. The matrices F (the new infection terms) and V (the transition terms) are given respectively, by Since κU 0 A + d a + ϕ = fuϕθ d f at E 01 , it follows that the basic reproduction number, denoted by R 0 , is given by Here, R 0 is the product of the birth rate, θ(1 − α) per unit time, the maternal transmission rate, v, the average life time of Wolbachia-infected female mosquitoes, 1/d f w , and the fraction f w ϕ/(κU 0 A + d a + ϕ) of immature mosquitoes surviving through the aquatic stages and becoming female mosquitoes. Thus, R 0 can be interpreted as the average number of secondary infections due to an infective individual during the infectious period, when everyone else in the population is susceptible [17].
Therefore, by Theorem 2 of [41], we establish the following results.
Theorem 3.1. Suppose that R 0u > 1. Then the Wolbachia-free equilibrium E 01 exists, and it is locally asymptotically stable if R 0 < 1, and unstable if R 0 > 1.
The condition R 0 > 1 is sufficient for the persistence of Wolbachia. It is easy to see that R 0 > 1 is equivalent to θ(1 − α)vf w /d f w > θf u /d f , the biological implication is that infected females produce more infected daughters than those produced by uninfected females. This is in line with some previous studies (see, for instance, [10,23]). On the other hand, we show in the next section that Wolbachia may still persist even if R 0 < 1 due to the backward bifurcation.

3.2.
Wolbachia-infected equilibria. We next calculate the Wolbachia-infected equilibrium (WE) at which the sizes of Wolbachia-infected mosquito population classes are not zeros.
We define the basic reproduction number of infected mosquito population as It is the product of the fecundity rate of infected females, (1 − α)θ, the maternal transmission rate, v, the fraction f w ϕ/(d a + ϕ) of immature mosquitoes surviving through the aquatic stages and becoming female mosquitoes, and the average life time of infected females, 1/d f w .
If the maternal transmission of Wolbachia is complete, i.e. v = 1, by a straightforward caculation, model (1) has another boundary equilibriun, i.e. all-Wolbachiainfected equilibrium (AWE) When the variables are ordered as ( and all the eigenvalues of D have negative real parts. Moreover, it is easy to prove that all the eigenvalues of A have negative real parts if R 0 > 1 − q and A has a positive eigenvalue if R 0 < 1 − q. Since B = 0, then the eigenvalues of J are eigenvalues of A and D. Therefore, the following theorem is proven.
Theorem 3.2. Suppose that R 0w > 1 and v = 1. Then the boundary equilibrium E 02 exists, and it is locally asymptotically stable if R 0 > 1 − q, and unstable if R 0 < 1 − q.
Since R 0w = R 0u R 0 , if R 0u > 1 and R 0 > 1 hold simultaneously, then R 0w > 1 and R 0 > 1 − q. Therefore, if the Wolbachia-free equilibrium E 01 is unstable, another boundary equilibrium E 02 must exist and is locally asymptotically stable when v = 1.
Although O(0, 0, 0, 0, 0, 0) is not an equilibrium of Model (1) due to the singularity, we prove that both kinds of mosquitoes die out under certain conditions in the following. Theorem 3.3. We assume that R 0u < 1 and R 0w < 1. Then the solution of system (1) tends to O as t → ∞ if all of the six variables are sufficiently small initially.
Proof. Linearisation of the last three equations of system (1) at O yields the following partially decoupled system. Obviously The Jacobian matrix of (5) at (0, 0, 0) is It is clear that all eigenvalues have negative real parts when f u ϕθ < d f (d a + ϕ), i.e., R 0u < 1, which implies that (0, 0, 0) is locally asymptotically stable for system (5). Thus, we conclude that the solution of the whole system (1) tends to O(0, 0, 0, 0, 0, 0) as t → ∞ if all of the six variables are sufficiently small initially.
In the following, we explore the existence of coexistence equilibria (CE) (1), i.e., the steady states where each component is strictly positive.
In terms of (1), we obtain Substituting (6) into the fourth equation of (1) yields which is equivalent to Substituting (6) into the first equation of (1) yields Let From (8), we obtain and which is equivalent to where From (8) and (12), we find that the number of coexistence equilibria of Model (1) is the same as the number of positive roots of equation (12) under the condition of R 0w > 1, which is necessary for the existence of a coexistence equilibrium.
We now examine the possible existence of positive roots of (12) in two cases.
For complete maternal transmission of Wolbachia, the condition for the existence of coexistence equilibrium is likely to hold when the coefficient of CI is sufficiently close to 1. The biological implication is that Wolbachia infection may induce the increase of mortality rate due to fitness cost if the effect of CI is large enough. When the unique coexistence equilibrium exists, i.e., R 0w > 1 and 1 − q < R 0 < 1, by Theorem 3.1 and Theorem 3.2, two boundary equilibria E 01 and E 02 exist and they are both locally asymptotically stable. If an unstable Wolbachia-free equilibrium E 01 exists, i.e., R 0u > 1 and R 0 > 1, according to Theorem 3.2 and Theorem 3.4, a stable all-Wolbachia-infected equilibrium E 02 exists and there is no coexistence equilibria, which is consistent with elementary competition theory. Our results are similar to the results of models studied in [13,23].
When R 0 < 1, then C 3 > 0. (12) has at least one positive root provided If ∆ > 0, there exist two roots which implies the existence of two coexistence equilibria. If ∆ = 0, the two roots coalesce and the unique root is of multiplicity 2. The condition (14) is abstract, we will examine condition (14) in terms of model parameters in more details later. We first summarize the above results in the following theorem. 1. If R 0 > 1 and R 0w > 1, there exists a unique coexistence equilibrium E * 1 , 2. If R 0 = 1 and R 0w > 1, there exists a unique coexistence equilibrium E * 3. If R 0 < 1 and R 0w > 1, two coexistence equilibria E * 1 and E * 2 exist, if and only if C 2 < 0 and ∆ > 0, And these two coexistence equilibria coalesce if and only if C 2 < 0 and ∆ = 0, 4. Otherwise, there is no coexistence equilibria, where E * 1 and E * 2 are corresponding to χ * 1 and χ * 2 , respectively.
3.3. Backward bifurcation. According to (12), the ratio of the natural death rate of Wolbachia-free mosquitoes d f to the natural death rate of Wolbachiainfected mosquitoes d f w , and the ratio of the fraction of births that are female mosquitoes for Wolbachia-infected mosquitoes f w to the fraction of births that are female mosquitoes for Wolbachia-free mosquitoes f u are crucial factors for the existence of coexistence equilibria. Let Then R 0 = (1 − α)vξ and the coefficients of (12) become Next we address how the four critical parameters τ, v, q, and ξ affect the occurrence of backward bifurcation when 0 < v < 1 and v = 1, respectively. ( then C 2 > 0 and there exists a positive root of (12) iff R 0 > 1, i.e. forward bifurcation occurs at R 0 = 1. In other words, is a necessary condition for backward bifurcation. In order to describe the backward bifurcation, we assume that and choose ξ and q as parameters. We will describe the conditions (14) under above assumption and then present the bifurcation diagram in the positive quadrant of the (ξ, q) plane where 0 ≤ q ≤ 1. First, we note that C 2 = 0 defines a line, denoted by L 1 (see Fig. 1): This line intersects with the ξ − axis at the point K 1 (ξ 1 , 0), where .
Moreover, it intersects with the line q = 1 at the point K 0 (0, 1). The basic reproduction number R 0 = 1 defines a vertical line L 2 (see Fig. 1), It is easy to check that L 2 has an intersection with the ξ−axis at the point K 2 (ξ 2 , 0), where Since then the point K 1 is always located to the right hand side of We can verify that L 2 also intersects with the line L 1 and q = 1 at K 3 (ξ 2 , q 3 ) and K 4 (ξ 2 , 1), respectively, where Let the curve defined by ∆ = 0 be L 3 (Fig. 1). A straightforward calculation gives where and It is easy to check that L 3 passes through K 3 . A straightforward calculation can also verify that L 3 intersects with the line q = 1 at the point K 0 (0, 1).
Based on the above discussion and Theorem 6, then we have the following theorem regarding the backward bifurcation and distribution of WE.
Then Model (1) undergoes backward bifurcation on the line segment L 2 between K 3 and K 4 . Moreover, in the region D 1 where R 0 > 1, there exists a unique coexistence equilibrium; in the region D 2 , there are two coexistence equilibria; and on the segment of L 2 between K 3 and K 4 , i.e. D 3 , and the segments of the curve L 3 between K 0 and K 3 , i.e. D 4 , there also exists a unique coexistence equilibrium, where The corresponding bifurcation diagram is given in Fig. 1.
From the above analysis, we know that if ξ = f ∆ (q) and q 3 < q ≤ 1, i.e., for (ξ, q) ∈ D 4 , two coexistence equilibria of Model (1) coalesce. LetR 0 = R 0 | (ξ,q)∈D4 , then on L 3 , we havê We summarize the backward bifurcation in terms of the basic reproduction number in the following theorem. Fig. 1, it does not have any coexistence equilibrium point if R 0 <R 0 ; there are two coexistence equilibria ifR 0 < R 0 < 1; there is a unique endemic equilibrium if R 0 > 1. A backward bifurcation occurs at R 0 = 1 (see Fig.  2).
(II): v = 1 When v = 1, we also choose ξ and q as parameters to describe the distribution of WE and then present the bifurcation diagram in the positive quadrant of the (ξ, q) plane where 0 ≤ q ≤ 1.
There are two coexistence equilibria in D 2 whereR 0 < R 0 < 1, and the two equilibria coalesce when the parameters approach the curve segment of L 3 between K 0 and K 3 where R 0 =R 0 . There is a unique coexistence equilibrium on the segment of L 2 between K 3 and K 4 , where R 0 = 1. In the other regions of the first quadrant, no WE exists.  Table 4 except that d f w varies.
According to Theorem 3 and Theorem 5, we can get the following theorem regarding the distribution of WE when v = 1.
Theorem 3.8. Suppose that v = 1 and R 0w > 1. In the region D 11 , where R 0 < 1 − q, there exists a unique unstable Wolbachia-infected equilibrium (WE), in the region D 22 where 1 − q < R 0 < 1 there are two WE. One is coexistence equilibria and the other is the boundary equilibrium AWE, which is stable. Besides, in the region D 12 , there also exists a unique stable WE, where The corresponding bifurcation diagram is given in Fig. 3. Figure 3. The distribution of Wolbachia-infected equilibria (WE) in the plane of (ξ, q) when v = 1 and R 0w > 1. In the region D 11 , where R 0 < 1 − q, there exists an unstable WE E 02 . In the region D 22 , where 1−q < R 0 < 1, There are two WE. One is a coexistence equilibria E * and the other is the boundary equilibrium, E 02 , which is stable. In the region D 12 , there also exists a unique stable WE E 02 .
Then, we can summarize the results regarding existence and stability of equilibria in Table 3. 4. Numerical analysis. In order to study the stability of the coexistence equilibria (CE) and the impact of some critical factors of Wolbachia-infection, we present some numerical analysis using the software Matlab. 4.1. Stability of coexistence equilibria and different cases. As mentioned above, the Wolbachia-free equilibrium E 01 is unstable when R 0 > 1 and Wolbachia infection can be established (Table 3: Case 1, Case 2, and Case 5). If R 0 < 1, the essential condition for Wolbachia establishment is the existence of Wolbachiainfected equilibria E * or E 02 (Table 3: Case 3 and Case 7) when ignoring the singularity cases.
To analyze local stability of these coexistence equilibria (CE) E * , one may apply a standard technique based on calculation of eigenvalues of the Jacobian matrix v Threshold condition   evaluated at the steady states. However, this approach looks rather knotty and impracticable from the computational standpoint [5]. Alternatively, one can use Monte Carlo method (see [5] and references therein) to repeatedly calculate the eigenvalues of the Jacobian evaluated at each steady state E * 1 , E * 2 , and E * c . According to (6), (10), (13), and (15) the coordinates of all three steady states can be expressed in terms of fourteen parameters (θ, q, v, d a , d f , d f w , d m , d mw , ϕ, α, f u , f w , τ , κ) whose baseline values and ranges are given in Table 4. The sampling pool + was defined by a Cartesian product of fourteen intervals P i , i = 1, ..., 14 listed in Table 4. Our sampling comprised 10 6 confounding scenarios S = (s 1 , ..., s 14 ) ∈ S, where each s i ∈ P i , i = 1, ..., 14, was randomly chosen under uniform distribution with no correlation between parameters. Simulation results are presented in Fig. 4 and Fig. 5.
Excluding all singularity cases (real parts of some eigenvalues are zero), as is shown in Fig. 4, all the real parts of eigenvalues corresponding to E * 1 are negative, which implies that E * 1 is stable. Fig. 5 shows that E * 2 and E * c are unstable because the products of all six eigenvalues are always negative except singular cases.
Therefore, for the two scenarios in Case 3 and 7, Model (1) exhibits bistability phenomenon. Whether Wolbachia persists or not depends on the initial values of the state variables involved. Complete population replacement is possible when v = 1, as long as the initial values lies in the basin of attraction of E 02 . We focus on them in the following discussion.

4.2.
The threshold for initial condition. The stable manifold of the saddle E * 2 and E * c are shown in the U A U F W M -space using parameter values given in Table 4 (see Fig. 6). In fact, this manifold can be depicted in any three-dimensional space of the model. We can also calculate the threshold for initial values of Wolbachiainfected males and females for Wolbachia invasion if other initial values are fixed. Fig. 7 manifests that the threshold line in the W F W M -plane isolates the plane into two parts. The initial values lying in the basin of attraction of the stable Wolbachiainfected equilibrium is essential to guarantee that Wolbachia can be established in a mosquito population. Moreover, the smaller the initial number of Wolbachiainfected female W F (0) is, the lager the initial number of Wolbachia-infected male W M (0) is needed. It provides an approach for evaluating the threshold number of infected mosquitoes needed to be released for one-time augmentation strategy.  Table 4. In order to demonstrate the significance of augmentation of infected females, we compared three different cases in Fig. 8 (I) No Wolbachia-infected mosquitoes are released.
(II) 1,500,000 Wolbachia-infected male mosquitoes are released for 300 days and no female mosquitoes are released.
(III) One Wolbachia-infected female and 150,000 infected males are released once at the beginning.
In (III), as in [23], we give an extreme case that only one infected female mosquito is released and we assume that the infected female mosquito can survive before laying eggs. It is our aim to evaluate the significance of the augmentation of infected females even though the quantity is tiny. Fig. 8 (a) and (c) imply that although releasing infected males can reduce the number of wild females at the early stage, the number of wild females goes up to the level of case (I) gradually once the augmentation is terminated. Nevertheless, even if only one infected female is released at the beginning in case (III), Wolbachia establishment and population replacement can be realized as long as enough number of males are released simultaneously. In order to dispel the doubt of the release of infected females, the time courses of the number of wild females and the total number of females are simulated in case of imperfect and perfect maternal transmission, correspondingly. Fig. 8 (b) and (d) demonstrate that the total number of females is less than that of case (I) in the whole process, so that people do not need to worry about female mosquito augmentation as long as the number of females released is small enough. Thus, releasing as few infected females as possible (at least one female) and a large enough quantity of infected males is a desirable strategy for the perspective of economy and disease control [23]. This conclusion is in line with the conclusion obtained in [23]. The threshold value of infected males required to be released can be obtained by Fig. 7. 4.3. Sensitivity analysis. Sensitivity analysis (SA) is defined as the study of how uncertainty in the quantity of interest(QoI) is attributed to different sources of uncertainties [34,35,39] etc. It helps us understand how the parameters of a model affect the quantity of interest(QoI). SA is often employed to rank the model parameters in the order of their significance so that we can determine the most significant parameters. We performed global sensitivity analysis to identify the key factors that determine the magnitudes of R 0 and W * F /U * F , which manifest the level of population replacement (see Fig. 9).
By the expression of R 0 , W * F and U * F , we can find that where χ * denotes the positive roots of Equation 12. Hence R 0 and W * F /U * F are independent of d a , ϕ, θ and κ, which are not considered in Fig. 9.
We performed PRCC analysis with 10 6 random sample uniformly distributed in the range of parameter values from Table 4. Fig. 9 demonstrates that the basic reproduction number R 0 is the most sensitive to the death rates of wild female mosquitoes (d f ) and wild female mosquitoes (d f w ), whose PRCC indices are 0.9387       Table 4. There are four critical parameters for Wolbachia infection, that is, MK effect parameter τ , CI effect parameter q, maternal transmission rate v, and fertility cost coefficient α. As shown in Fig. 10, there are threshold conditions of these four parameters for the existence of Wolbachia-infected equilibria. If one of the values of τ and α is too large, or one of the values of v and q is too small, Wolbachia establishment will not be achieved. Moreover, the number of Wolbachia-infected females at the stable coexistence equilibrium will decrease with the increase of τ and α or the decrease of v and q. Besides, Fig. 11 depicts the threshold values of Wolbachia-infected males with respect to the four parameters for Wolbachia establishment. The results manifest that the Wolbachia strains with smaller τ and α or larger v and q may lead to less infected males to be released.

4.4.
The impact of CI effect suffering by uninfected ova produced by infected females. In order to study the the impact of the CI effect suffering by uninfected ova produced by infected females, we compare the threshold lines of initial values for Wolbachia esbabilishment for λ = 0 and λ = q. Fig. 12 shows that, if maternal transmission is imperfect, neglecting the CI effect suffering by uninfected ova produced by infected females, i.e. λ = 0 in model (1), may result in the overestimation of the number of infected mosquitoes needed to be released. 5. Discussion. In this paper, we developed a novel Wolbachia infection model in a two-sex mosquito population with stage-structure. The model captures many of the key effects of Wolbachia infection, including incomplete cytoplasmic incompatibility (CI), male killing (MK) effect, maternal transmission, fertility cost of the Wolbachia infection to reproductive output. Particularly, the CI effect suffering by uninfected ova produced by infected females, which is often neglected, is incorporated.
We provided the existence and stability of the equilibrium points associated with the proposed model for both complete and incomplete maternal transmissions. For both cases, there is a Wolbachia-free equilibrium (WFE), which is stable when 0 < R 0 < 1. When the maternal transmission is complete, there is a stable all-Wolbachia-infected equilibrium (AWE), and a unstable coexistence equilibrium (CE). Meanwhile, when the maternal transmission is incomplete, there is no AWE and only two CE, one is stable and the other is unstable. Under this circumstance, bistability occurs for both cases leading to backward bifurcation. We can find the stable manifold of the unstable saddle that isolates the state space into two parts, which are the basins of attraction of the two stable states, respectively. One is Wolbachia-infected state and the other is Wolbachia-free state. Whether Wolbachia persists or not depends on the initial values of all the state variables involved. By the stable manifold, we provided the method to evaluate the threshold number of infected mosquitoes required to be released at the beginning. Besides, considering the perspective of economy and disease control, a desirable strategy is keeping the number of infected female mosquitoes released to be a necessary minimum by replying on higher number of male mosquitoes released. Besides, people do not need to worry about female mosquito augmentation as long as the number of females released is as few as possible. This result is in line with that provided in [23].
Global sensitivity analysis and numerical simulation have been performed to study the impact of model parameters to Wolbachia establishment and population replacement. The results manifest that different parameter values with different Wolbachia strains may lead to different Wolbachia establishment results. There are four critical parameters for Wolbachia infection, including MK effect τ , CI effect q, maternal transmission rate v, and fertility cost coefficient α. If one of the values of τ and α is too large, or one of the values of v and q is too small, Wolbachia establishment will not be achieved. Moreover, the number of Wolbachia-infected females at the stable coexistence equilibrium decreases with the increase of τ and α or the decrease of v and q. Besides, the threshold values of Wolbachia-infected males with respect to the four parameters for Wolbachia establishment, smaller τ and α or larger v and q may lead to less infected males to be released. Therefore, low levels of MK effect and fitness costs as well as high levels of CI effect and maternal inheritance are in favor of Wolbachia establishment. Suitable strains of Wolbachia should be chosen deliberately for the success of population suppression or replacement as well as the control of some mosquito-borne disease transmission.
It is worth emphasizing that ignoring the CI effect suffering by uninfected ova produced by infected females, i.e. λ = 0 in model (1), may result in the overestimation of the number of infected mosquitoes needed to be released. Therefore, considering the CI effect suffering by uninfected ova produced by infected females is very essential.
As a preliminary work, we are trying to model the impact of all the critical factors on Wolbachia establishment in wild mosquito population, especially the CI effect suffering from uninfected ova produced by infected females. Nevertheless, mosquito population may be affected by some other environmental and climatic factors, such as humidity and precipitation. Moreover, our results are lack of detailed comparison of how different Wolbachia strains lead to different Wolbachia establishment, which will be incorporated in future work.