Multiscale modelling the effects of CI genetic evolution in mosquito population on the control of dengue fever

Endosymbiotic Wolbachia bacteria are widely applied for the control of dengue fever by manipulating the reproductive mechanism of mosquitoes, including maternal inheritance and cytoplasmic incompatibility (CI). CI means that the offsprings from the matings between Wolbachia infected males and uninfected females can not be hatched. At present, CI effect is assumed as a constant in most of dynamic systems for the spread of Wolbachia. However, their spread may arouse the evolution of mosquitoes to resist CI. Thus, a multiscale model combining a birth-pulse model with a gene-induced discrete model for the frequencies of alleles is proposed to describe the spread of Wolbachia in mosquito population with resistance allele of CI. The main results indicate that the strategy of population eradication can not be realized, while the strategy of population replacement may be realized with the success of sensitive or resistance allele. If appropriate Wolbachia strains can not be selected, then there is a high probability of the failure of population replacement. Moreover, Wolbachia-induced parameters may arouse the catastrophic shifts among stable states of the model. In addition, the demographic parameters and Wolbachia-induced parameters may affect the level and the speed of population replacement and the density of uninfected mosquitoes.

Dengue fever as one type of vector-borne diseases, is caused by dengue virus mainly transmitted by mosquitoes 1,2 . This disease seriously threaten human health, especially in tropical and subtropical regions due to geographical distribution of mosquitoes, climatic characteristics and human movement patterns 3,4 . At present, endosymbiotic Wolbachia bacteria has been discovered to control dengue diseases as a new promising biological approach. Wolbachia is preferable to females and can be maternally transmitted from female mosquitoes to their offsprings 5,6 . Wolbachia exhibits a special phenomenon during the mating of their hosts, named as cytoplasmic incompatibility (CI) 7,8 . CI can be expressed as the death of the mosquito embryos when uninfected females mate with Wolbachia infected males whereas other types of matings are not affected 9,10 .
Owing to maternal inheritance and CI mechanism, Wolbachia can only vertically spread from females to their offsprings and terminate in infected males, which benefits to enhancing the fitness of infected hosts and reducing the proportion of uninfected hosts 11,12 . So when Wolbachia is utilized in the control of vector-borne diseases, we hope to realize two strategies, one is population suppression (eradication) and other is population replacement (Wolbachia infected insects established and replacing uninfected ones). Since Wolbachia bacteria is widely applied for the control of dengue fever, various mathematical models [13][14][15][16][17] have been established to explain the dynamics of Wolbachia invasion and the change on the size of host population. Especially, Zhang et al. 18,19 developed impulsive models to investigate the effects of different Wolbachia infected mosquito augmentations, including releases with different sex and release quantity, release moment and period, release sequence and times, on the success of different strategies of mosquito control.
In general, the spread of Wolbachia is not complete and often has a cost to offset host fitness because of intragenomic conflict and coevolution [20][21][22][23] . Martłnez-Rodrłguez et al. 20 proposed a matrix model to describe the variation of infection proportion in the host life cycle and predicted the evolution progress, and showed that the long-term infection was affected by the change of Wolbachia proportion. Turelli 21 proposed that the host genome should evolve to ameliorate action of Wolbachia, further he analyzed selection on Wolbachia and their host genes by infection frequency models. Cooper et al. 24 discovered that Wolbachia can cause intraspecific and interspecific CI based on the three hybridizing species of Drosophila, and both host genotype and Wolbachia variation modulate CI expression.
In addition, the dynamics of CI-expressing Wolbachia infections of insects have been extensively modeled [25][26][27] . Assuming random mating and discrete generations on host population, Caspari et al. 25,26 proposed difference equations to describe the frequency of infected adults at generation t + 1. They proved that there exist infected adult eradicated and established equilibria, and for Wolbachia spreading into the host population it needs to be introduced at an initial prevalence greater than the value of unstable equilibrium. In D. melanogaster, male development time can greatly influence CI expression discovered by Yamada et al. 28 . Thus, the infected males taking longer to undergo larval development quickly lose their ability to express the CI phenotype. Not only that, the level of CI expression in Drosophila melanogaster declines rapidly with male age, particularly when males are repeatedly mated 29 . Besides, Hurst et al. 30 improved the recursions to describe the dynamics of CI, and pointed that following the invasion of CI inducing Wolbachia into an uninfected population, not only will the sterilizing effect wane but the conditions become permissive for the spread of the uninfected cytotype.
Based on above works, it is interesting to investigate questions about the evolution of infections and their hosts, combining with their breeding process. However, a population-genetic model described the spread of insecticide resistance on mosquito population had been developed by Birget et al. 31 . In fact, the development of mosquito resistance to insecticide and resistance of CI are similar, we aim to explore the effect of resistance development of mosquito populations to CI mechanism on the spread of Wolbachia among mosquito populations by mathematical method. To do this, we assume that the evolution of CI mechanism is determined by a single gene with two alleles which conforms to Hardy-Weinberg Principle 32 . Thus the frequencies of alleles in infected and uninfected mosquitoes will change with the reproductive process of mosquito populations. Note that four factors may affect the reproductive progress, including maternal inheritance, resistance level, the dominance level of gene and the fertility cost. Although the frequencies of alleles are independent to the density of host population, while the latter may be impacted by the former. Therefore, we combine a birth-pulse model 33 of mosquito populations with a discrete model of the frequencies of alleles, i.e. the multiscale models have been developed in the present work.

Models and Methods
Birth-pulse model. In order to study the effect of the evolution of Wolbachia-induced CI on the dynamic behaviors of mosquito populations, we first introduce the following birth-pulse model [33][34][35] where n(n = 1, 2, …) denotes the sequence of birth pulse and T is the breading period. N(t) is the total density of mosquito population, which includes the infected I(t) and uninfected U(t). The densities of infected and uninfected mosquitoes are decreased owing to density-dependent death whenever t ≠ nT, while birth pulses with the CI mechanism and matrilineal inheritance occur at t = nT. The items I(nT + ) and U(nT + ) denote the mosquito densities of the infected and uninfected ones at time nT after birth pulse. The description of parameters are shown in Table 1. For more details of above model, please refer to the work in 33 .
Birth-pulse model with CI genetic evolution. Suppose that the spread of Wolbachia and the evolution of CI in mosquito population are determined by a single gene with two alleles, including dominant (here resistance) allele R and recessive (here sensitive) allele S. It gives rise to three different genotypes in mosquito population, including resistant homozygote (RR), heterozygote (RS) and sensitive homozygote (SS).
The frequencies of three genotypes. The genotype frequencies for both males and females are assumed to be equal. The frequencies of alleles R and S in Wolbachia infected (uninfected) mosquitoes at generation n are denoted as p I,n and q I,n (p U,n and q U,n ), respectively, then we have p I + q I = I/(I + U) and p U + q U = U/(I + U). Thus, if the density of mosquitoes is zero, then p I,n + q I,n + p U,n + q U,n = 0, otherwise p I,n + q I,n + p U,n + q U,n = 1. According to Hardy-Weinberg Principle and maternal inheritance, the frequencies of alleles satisfy the equation ((p I,n + p U,n ) + (q I,n + q U,n )) 2 = 1. The fitnesses of three genotypes are different at generation n owing to CI effect in the progress of genetic evolution, which causes the frequencies of alleles in infected and uninfected mosquitoes change constantly. Thus, it is necessary to investigate the frequency dynamics of alleles. The frequencies of different genotypes with different mating types P i (i = 1, 2, 3) at generation n are shown in Table 2.
The fitnesses of three genotypes. Note that Wolbachia infected females never undergo CI mechanism, and all the female mosquitoes mating with uninfected males can not be affected by CI. For convenience, these matings without CI mechanism are renamed as group one. However, uninfected females mating with infected males may arise the occurrence of CI which is renamed as group two. The reproduction of mosquito, which is determined by the natural birth rate b(N(t)) and their resistance of CI effect, called fitness 31 . In the group one, the effect of resistant homozygote individuals RR on the fertility cost of resistance is assumed as z, then the cost in heterozygote individuals RS is the product of z and the level of dominance h, homozygote sensitives is not affected by the resistance and reproduce normally. For the group two, we assume that all the sensitive mosquitoes SS can not reproduce and not be affected by the resistance, that is complete CI mechanism will occur. There is a survival of resistant homozygote RR because of the resistance ρ and the cost of fertility is reduced by z. The fitness in heterozygote individuals RS is the product of ρ and the level of dominant allele h, meanwhile, the cost of fertility is also reduced by hz. Therefore, the fitnesses of different genotypes in different groups at generation n are calculated as shown in Table 3, with all the parameter values of ρ, h and z belonging to [0,1]. The detailed flow diagram of genetic evolution is shown in Fig. 1.
The genotype frequencies of new born offsprings. The genotype frequencies of new born offsprings at generation n + 1 are determined by the frequencies and the fitnesses of three genotypes RR, RS and SS. Infected females occur maternal inheritance with a probability τ (also called as transmission rate τ for convenience), or in other words, infected females can produce a proportion 1 − τ of uninfected offsprings. So the ratios of new born infected and uninfected mosquitoes from three different genotypes to the total mosquito population at generation n are shown in Table 4.
The frequencies of two alleles. Since the development of mosquito populations is a process of overlapping generation, the survival mosquitoes from generation n may affect the frequencies of alleles R and S in infected and uninfected mosquitoes at generation n + 1. So the frequencies of alleles at generation n + 1 consist of the survival mosquitoes from generation n and the new born offsprings at generation n + 1. Therefore, the process of the frequencies of alleles in infected and uninfected mosquitoes from generation n to n + 1 is depicted by Fig. 2    where Q is the total birth rate of Wolbachia infected and uninfected mosquitoes at generation n with genotypes RR, RS and SS, while Q I,RR , Q I,RS and Q I,SS (Q U,RR , Q U,RS and Q U,SS ) denote the frequencies of genotypes RR, RS and SS for new born Wolbachia infected (uninfected) offsprings at generation n + 1, respectively. Based on (2.1), (2.2) and Table 2, a multiscale model combining the birth-pulse model with the genetic evolution to resist CI mechanism in mosquito populations is proposed as follows , ,

Genotype
Group one (W 1 ) Group two (W 2 ) Table 3. The fitnesses of three genotypes in different groups.

Figure 1.
Flow diagram of the fitnesses from different mating types of Wolbachia infected and uninfected mosquitoes at generation n. Total density of mosquito population is divided into four classes, including infected females, infected males, uninfected females and uninfected males, with three genotypes RR, RS and SS. Assume that there is the same sex ratio in infected and uninfected mosquitoes. The fitnesses W i (i = 1, 2) from three mating types P 1 , P 2 and P 3 with three genotypes are shown in Table 2. When infected females mate with males (P 1 ), their offsprings may be infected owing to maternal inheritance; when uninfected females mate with uninfected males (P 2 ), their offsprings are uninfected; while when uninfected females mate with infected males (P 3 ), their offsprings may not be hatched owing to CI mechanism, so the three mating types P 1 , P 2 and P 3 represent maternal inheritance, normal and CI effect groups, respectively. The densities of mosquitoes at generation n + 1 are composed of the survivals of mosquitoes at generation n and new born offsprings from generation n.
where one scale is based on the dynamic system of mosquito populations, and another is based on the frequencies of alleles in mosquito populations. The parameter interpretations and notations in the model are given in Table 1.

Stroboscopic map of system (2.3). The stroboscopic map which plays an important role in investigating
the dynamic behavior of system (2.3) is determined by the densities of infected and uninfected mosquitoes after each birth pulse at discrete time nT. Note that we aim to study the effect of CI genetic evolution on the density of mosquito population. For simplicity, birth rates of both infected and uninfected mosquitoes are assumed to be equal and as a constant b. In addition, we mainly study the case when extra death rate of infected mosquitoes D = 0, which only affects the change rates for the densities of mosquito populations, but does not change the qualitative behavior of system (2.3) 17 .
To obtain stroboscopic map of system (2.3), we first solve the two differential equations of system (2.3) in any birth pulse interval (nT, (n + 1)T] as follows , and deform some equations, we obtain the corresponding stroboscopic map of system (2.3) as follows   where Based on system (2.5), the densities of infected and uninfected mosquitoes depend on variables p i and q i (i = I, U), while the frequencies of alleles R and S in infected and uninfected mosquitoes are independent to the densities of infected and uninfected mosquitoes. Therefore, system (2.5) is divided into two subsystems, and the first two equations are renamed as subsystem (I) and the other four equations as subsystem (II). For convenience, denote the right hand functions of subsystem (I) as f i (i = 1, 2), while the right hand functions of subsystem (II) as f i (i = 3, 4, 5, 6). Since subsystem (II) is independent to subsystem (I), we first analyze it in the coming section.

Results
The analysis of subsystem (II). The existence and stabilities of equilibria for subsystem (II) in different cases are shown in Table 5, while for their detailed calculation and analysis, see in Supplementary Information (SI). The expressions of equilibria for subsystem (II) are listed as follows In the following, we mainly analyze subsystems (I) and (II) without fertility cost, and the effects of fertility cost and dominance level on the speed of population replacement are considered in subsection 3.5. In order to understand the existence of equilibria of subsystem (II) in more detail, we take τ and ρ as bifurcation parameters, and four key lines are defined as follows L L L L : 1, : 0 75, : 4 3, : 1 Then the four lines (L i , i = 1, 2, 3, 4) divide the τ and ρ parameter space into eight regions as shown in Fig. 3.
i coexist in regions Ω 1 , Ω 2 and Ω 3 , respectively. Especially, equilibria coincide in lines L 2 and L 3 , respectively. Next we state the biology significance for the five stable equilibria E 0 Only in the neighbourhood of E 0 ⁎ , the solutions of subsystem (II) tend to E 0 ⁎ , so it is difficult to make the frequencies of alleles in infected and uninfected mosquitoes tend to zero at the same time, that is to say, it is difficult to eradicate infected and uninfected mosquitoes simultaneously. Thus in the following, we focus on the other four equilibria. When the solutions of system (II) stabilize at equilibrium E 1 ⁎ , then all alleles to depict CI effect are in uninfected mosquitoes, which indicates that some uninfected mosquitoes have resistance allele, while the rest of uninfected mosquitoes only have recessive allele. The stability of equilibrium ⁎ E 2 indicates that the success of recessive allele in mosquito population, and it means that all mosquitoes don't have resistance allele. However when the solutions of system

Basins of attraction of subsystem (II).
Based on the above subsection, the stabilities of equilibria of subsystem (II) have been determined by analytical and numerical techniques. In order to study the effect of different initial frequencies of alleles R and S in infected and uninfected mosquitoes on the final distribution of the allele frequencies, it is necessary to further study the effects of different initial values on the solutions of subsystem (II).
To do this, all parameter values are fixed, except ρ, and a wide ranges of frequencies of alleles R and S in infected and uninfected mosquitoes are chosen, then we obtain the basins of attraction of equilibria for subsystem (II), as shown in Fig. 4. When parameter values are chosen from region Ω 2 , then it follows from Figs 3 and 4(a) that the two stable equilibria ⁎ E 1 and E 2 ⁎ coexist. Consequently the solutions of subsystem (II) from the initial frequencies of alleles lying in deep red or green region will stabilize at equilibrium E 1 ⁎ or E 2 ⁎ , respectively, which indicates that the two regions are the basins of attraction of the two equilibria, respectively. Especially, in this case, although there exist initial frequencies of alleles in infected mosquitoes, most solutions of subsystem (II) will tend to equilibrium ⁎ E 1 such that the frequencies of alleles tending to zero in infected mosquitoes has a high probability 99.81%, which is obtained by calculating the ratio of the number of points in deep red region to that of in region Ω. All those confirm that a high probability for the failure of population replacement.
Similarly, Figs 3 and 4(b) show that the tri-stable equilibria E 1 ⁎ , E 2 ⁎ and ⁎ E 4 coexist when parameter values are chosen from region Ω 3 , then the solutions of subsystem (II) from the initial frequencies of alleles lying in deep red, green or blue region will stabilize at equilibrium ⁎ E 1 , ⁎ E 2 or E 4 ⁎ , respectively, which indicates that the three regions are the basins of attraction of the three equilibria, respectively. In this case, most of solutions of subsystem (II) will tend to equilibrium ⁎ E 1 or E 4 ⁎ such that there is a small probability (9.87%) for the frequencies of alleles tending zero in infected mosquitoes or a significant probability (89.94%) for the success of resistance allele in mosquito populations. All those confirm that we will fail to realize population replacement but can realize partial population replacement with resistance allele. Figure 4(a,b) show that there is a low probability for the solutions of subsystem (II) stabilizing at equilibrium ⁎ E 2 , that is to say, the probability for the establishment of populations with sensitive allele owing to their low fitness to CI is low. Fig. 3, the existence of stable equilibrium ⁎ E 1 , the coexistence of two stable equilibria ⁎ E 1 and ⁎ E 2 , and the coexistence of three stable equilibria E 1 ⁎ , E 2 ⁎ and ⁎ E 4 are in regions Ω 1 , Ω 2 and Ω 3 , respectively. In order to investigate the distribution of allele frequencies, it is important to study the effects of different parameters and initial frequencies on the coexistence, catastrophic shifts and the values of three stable equilibria. Note that for the special case z = 0, the level of dominance h does not appear in the expressions of equilibria, so these equilibria are mainly determined by parameters τ and ρ. Therefore, we separately fix the two parameters, and plot one or two-parameter bifurcation diagrams from three different initial values, as shown in Figs 5, 6, 7 and 8.

Parameter bifurcation analyses of subsystem (II). According to
From Fig. 5, when the transmission rate is low (τ ∈ (0, 3/4]), then the frequencies of alleles in infected mosquitoes tend to zero for all three initial values. With the increase of transmission rate (τ ∈ (3/4, (3 + ρ)/4]), it occurs the coexistence of two attractors ⁎ E 1 and ⁎ E 2 . As the transmission rate further increase (τ ∈ ((3 + ρ)/4, 1]), the three attractors E 1 ⁎ , ⁎ E 2 and E 4 ⁎ can coexist. Especially, with the increase of transmission rate τ, the solution of subsystem (II) from initial frequencies of alleles (p I,1 , q I,1 , p U,1 , q U,1 ) = (0, 0, 0.1, 0.9) can stabilize at ⁎ E 1 , and the catastrophic shift of equilibria does not occur at all (see magenta curves); the stable states of subsystem (II) from initial frequencies of alleles (p I,1 , q I,1 , p U,1 , q U,1 ) = (0, 0.8, 0,   Similarly, it follows from Fig. 6 that when the resistance level is low (ρ ∈ (0, 4τ − 3]), then it occurs the existence of three attractors ⁎ E 1 , ⁎ E 2 and ⁎ E 4 . With the increase of resistance level (ρ ∈ (4τ − 3, 1]), it occurs the coexistence of two attractors E 1 ⁎ and ⁎ E 2 . Especially, with the increase of resistance level ρ, the solutions of subsystem (II) from the same initial frequencies of alleles as those of the deep magenta and blue curves in Fig. 5 stabilize at E 1 ⁎ and E 2 ⁎ , respectively, and there is no switch of equilibria (see magenta and blue curves); the stable states of subsystem (II) from the same initial frequencies of alleles as those of the red curves switch from E 4 ⁎ to ⁎ E 1 at ρ = 4τ − 3 (see red curves).
In the following, it is necessary to further study how initial frequencies of alleles affect the catastrophic shift and the values of stable equilibria of subsystem (II), as shown in Figs 7 and 8. It follows from Fig. 7, with the increase of transmission rate τ, the solutions of subsystem (II) from initial frequencies of alleles (p I,1 , q I,1 , p U,1 , q U,1 ) = (0, 0, 0.9, 0.1) stabilize at E 1 ⁎ , and the catastrophic shift of equilibria does not occur at all (see magenta curves); the stable states of subsystem (II) from initial frequencies of alleles (p I,1 , q I,1 , p U,1 , q U,1 ) = (0, 0.25, 0, 0.75) switch from ⁎ E 1 to E 2 ⁎ at τ ≈ 0.77 (see blue curves); the stable states of subsystem (II) from initial frequencies of alleles (p I,1 , q I,1 , p U,1 , q U,1 ) = (0.05, 0.05, 0.4, 0.5) switch from ⁎ E 1 to ⁎ E 4 at τ ≈ 0.93 (see red curves). Similarly, the bifurcation diagrams of the frequencies of alleles at equilibria of subsystem (II) from the three initial frequencies of alleles with respect to ρ are shown in Fig. 8. The stable states of subsystem (II) from some initial values may occur catastrophic change for resistance level ρ ∈ (0, 4τ − 3]. Based on the above bifurcation analyses, both parameters τ and ρ can significantly affect the switching behavior among the stable equilibria. Therefore, to further show how the two parameters τ and ρ affect the switching behavior, we take the two-parameter bifurcation analyses, as shown in Fig. S2. From which we can obtain the switching curves of stable equilibria from E 1 ⁎ to E 2 ⁎ (see from deep red region to green region) and from E 1 ⁎ to ⁎ E 4 (see from deep red to blue region). Note that the switching curves will alter with the change of initial frequencies of alleles, by comparing Fig. S2(a) with Fig. S2(b). By the same way, the τ − ρ plane can be drawn to obtain the switching curves of stable equilibria, corresponding to Figs 5, 6, 7 and 8, so here we omit it. Therefore, for any initial frequencies of alleles in mosquito populations, we can numerically obtain the corresponding switching curves of stable equilibria of subsystem (II), which is important for the control of mosquito populations.
Replacement level and the density of uninfected mosquitoes. In the above subsections, the dynamic behaviors of subsystem (II) have been investigated by theoretical and numerical methods. Moreover, the two variables of subsystem (I) dependent upon the four variables of subsystem (II). In fact, in order to fight against dengue disease, we aim to realize the strategy of population replacement with high level and high speed. So in the following, we will focus on studying the effects of subsystem (II) and different parameter values on the dynamic behavior of subsystem (I). The detailed analysis of the existence and stabilities of equilibria for subsystem (I) is given in SI. And the expressions and stabilities of equilibria for subsystem (I) are shown in Table 6. Especially, the stabilities of E 1 ⁎ and ⁎ E 6 imply the success of complete population replacement or the failure of population replacement, respectively, which is determined by initial frequencies of alleles.
In practice, Wolbachia infected mosquitoes should be released to realize the strategy of population suppression or population replacement. Based on Table 6, on the one hand, equilibrium E 0 ⁎ is locally stable, but not asymptotically stable. Thus, the success of population suppression occurs provided only the initial frequencies of alleles in the neighbourhood of E 0 ⁎ , which indicates that the probability of the failure of population suppression may be high. On the other hand, E 1 ⁎ and E i ( 2, 4) i = ⁎ coexist and are locally stable provided Δ 2 > 0, which respectively means either the failure of population replacement or the success of partial population replacement determined by initial frequencies of alleles. Therefore, we aim to realize partial population replacement with high level, meanwhile, to keep the density of uninfected mosquitoes as low as possible. Particularly, the complete population replacement may be succeeded, provided perfect transmission rate (τ = 1).
Denote ξ = I/(I + U) to exactly depict replacement level. Then corresponding replacement levels of equilibria i are shown in Table 6. Here we just focus on the replacement levels of stable equilibria = ⁎ and the density of uninfected mosquitoes of E 2 ⁎ is less than that of E 4 ⁎ . For convenience, replacement levels ξ ⁎ 2 and ξ ⁎ 4 are defined as high level with sensitive allele and low level with resistance allele, respectively.
In order to discuss the effect of each parameter on replacement level and the density of uninfected mosquitoes, all the parameters in system (2.3) are divided into two distinct types, one is Wolbachia-induced parameters, including τ, ρ, h, z, another is demographic parameters, including b, δ, T. Note that no matter subsystem (I) with or without fertility cost, the values of equilibria is not affected by parameters z and h except equilibrium E 1 (1) . Thus we do not consider the effects of the two parameters on the level of population replacement. According to the expressions of replacement levels 2 Wolbachia-induced parameters and demographic parameters. High replacement level 2 ⁎ ξ is increased with the increase of transmission rate. So the two parameter values are chosen and others are fixed, then we draw the contour plots of low replacement level 4 ⁎ ξ and the densities of uninfected mosquitoes ⁎ U 2 and U 4 ⁎ , as shown in Fig. 9. It follows from Fig. 9(a) that low replacement level ⁎ ξ 4 is increased with the increase of transmission rate τ and the decrease of resistance level ρ. Which indicates that high transmission rate and low resistance level contribute to improve low replacement level. Figure 9(b) shows that the density of uninfected mosquitoes is reduced with the decrease of transmission rate τ and the increase of resistance level ρ, which indicates that low transmission rate and high resistance level contribute to reduce the density of uninfected mosquitoes. So transmission rate and resistance level have reverse effects on replacement level and the density of uninfected mosquitoes. Therefore, in practice, suitable Wolbachia strains should be selected to balance the reverse effects so that population replacement can be realized with a relative high level and low density of uninfected mosquitoes. Figure 9(c,d) show that the densities of uninfected mosquitoes ⁎ U 2 (see red lines) and ⁎ U 4 (see black lines) are reduced with the decrease of birth rate b and the increase of death rate δ or reproductive period T.  Table 6. The expressions and stabilities of equilibria for subsystem (I). S: Stable, U: Unstable, Replacement speed. In order to fight against dengue disease, it is important to study the effects of initial densities of mosquito populations and parameter values on the speed of population replacement. Therefore we numerically obtain the effects of these factors on the speed of the solutions of subsystem (I) stabilizing at stable states, as shown in Fig. 10. It follows from Fig. 10(a), when the density of uninfected mosquitoes is assumed as a constant, then the solutions of subsystem (I) from the more initial infected mosquitoes can stabilize at stable equilibria with less time. Figure 10(b) reflects that short reproduction period T contributes to the high speed of population replacement. Similarly, Fig. 10(c,d) indicate that high transmission rate τ and low resistance level ρ contribute to the high speed of population replacement. Although fertility cost z and dominance level h have no effect on the values and stabilities of equilibria for subsystem (II), the two parameters may affect replacement level and speed, as shown in Fig. 10. Figure 10

Discussion
In this paper, based on Mendelian inheritance, we propose a multiscale model to study the effects of resistance allele of CI mechanism on dynamic behaviors of mosquito populations so that the strategies of population replacement with high level and speed can be succeeded. To establish the model, we first study the frequencies of three genotypes, the fitnesses of three genotypes, the genotype frequencies of new born offsprings and the frequencies of two alleles step by step. Then we formulate a six-dimension discrete system by the stroboscopic map   In practice, in order to control the spread of dengue virus, the strategy of population replacement should be succeeded at first. On the one hand, note that different parameter spaces of τ and ρ may lead to different cases for the coexistence of attractors. Therefore, based on special mosquito species, suitable Wolbachia strains need to be selected to ensure Wolbachia-induced parameters τ and ρ lying in the best region Ω 2 or the second region Ω 3 , so that high level of population replacement may be realized (see Fig. 3). However, suppose Wolbachia strain fixed, there is a low probability for the success of high level replacement with sensitive allele (see deep green attractive region of E 2 ⁎ in Fig. 4). In this case, Wolbachia strain should be reselected so that parameters τ and ρ lying in region Ω 3 , then corresponding mosquito control strategies should be carried out to drive initial frequencies of alleles in mosquito population lying in deep blue attractive region of E 4 ⁎ in Fig. 4(b), which indicates the success of low level replacement with resistance allele. Especially, if some Wolbachia can induce perfect transmission rate, then the success of complete population replacement or the failure of population replacement is determined by the initial frequencies of alleles. Similarly, corresponding mosquito strategies should be carried out to drive initial frequencies of alleles in mosquito population lying in the objective basins of attraction for the success of complete population replacement. On the other hand, the switching curves for the stable equilibria of subsystem (II) are determined by initial frequencies of alleles. Therefore, suppose mosquito population confirmed, that means the initial density of mosquito population and initial frequencies of alleles in mosquitoes are fixed. Then based on special mosquito species, suitable Wolbachia strains should be selected to ensure Wolbachia-induced parameters τ and ρ lying in deep green regions (or second in deep blue regions in Fig. S2, so that high level replacement with sensitive allele (or low level replacement with resistance allele) can be realized.
Further, we aim to realize high level of population replacement with high speed and low density of uninfected mosquitoes which are determined by Wolbachia-induced parameters, demographic parameters of mosquitoes or initial frequencies of alleles in mosquitoes. According to Fig. 9, high transmission rate τ and low resistance level ρ contribute to improve the level of population replacement, while low transmission rate τ, high resistance level ρ, low birth rate b, high death rate δ and long reproduction period T benefit to reducing the density of uninfected mosquitoes. It follows from Fig. 10(a-d), high initial density of infected mosquitoes I 0 , high transmission rate τ, low resistance level ρ and low reproduction period T benefit to improving the speed of population replacement. It is worth noting that we should balance the inverse effects of parameters τ,ρ and T on control objectives. Especially, although fertility cost and dominance level has no effect on the level of replacement, the high fertility cost and low dominance level contribute to the success of population replacement with high speed.
In most empirical and experimental works, the evolution of CI is estimated by the intensity of CI expression 24,28,36 . Most researchers try to study the effects of some phenotypic factors on the evolution of CI and the distribution of host population. In addition, most mathematical models about CI expression focus on its macroscopic effects on the dynamics of host population 18,19,25,26,33 . Weeks et al. 37 said that many factors can influence the expression of Wolbachia effects which is rarely considered in laboratory studies, and it is necessary to find a comprehensive approach to understand the host evolution. Thus, our hypothesis as well as multiscale model combine the host population and the microgenetic inheritance, which should be reasonable by theoretical support and actual possibility.
Note that the evolution of CI can greatly affect the spread and distribution of Wolbachia in some host populations. In this work, as far as we know, it is the first time to try exploring the relationship between mosquito population dynamics and CI genetic evolutionary dynamics. We assume that the introduction of Wolbachia infected mosquitoes in practice may lead to the evolution of mosquito population to resist Wolbachia-induced CI mechanism, which could affect the spread of Wolbachia in mosquitoes by the genetic development of the alleles related with CI mechanism. As it is the first step, we only focus on describing the intensity of CI expression by natural birth rate b, maternal transmission rate τ and the genetics parameters (here z, h and ρ) in different genotypes and different groups. In fact, the effects of the phenotypic factors on CI expression can be depicted by changing the values of these parameters in some extent.
This study shows that the evolution may weaken the ability of Wolbachia spread in mosquito popuation, and give some interpretation for the incomplete CI by genetic evolution, which may lead to the failure of Wolbachia-based control in some field trials. Consequently, the CI evolution in mosquito population should be aroused sustained attention. Tortosa et al. 38 showed that incompatible matings effectively lower the fertility of infected males, so that selection acting on the host genome should tend to reduce the expression of CI in males, for example, by reducing the density of Wolbachia infection in males before sexual maturation. Thus, in the further work, we will explore the evolution of mosquito population with sex structure and some phenotypic factors to resist Wolbachia-induced CI mechanism, and its effect on the the spread of Wolbachia in mosquito populations. In addition, we will enhance the cooperation with the department of disease control to collect relative data of mosquito population with Wolbachia infection, and explore and develop the model to control other vector-borne diseases.