MODELING THE DYNAMICS OF WOLBACHIA -INFECTED AND UNINFECTED A EDES AEGYPTI POPULATIONS BY DELAY DIFFERENTIAL EQUATIONS

. Starting from an age structured partial diﬀerential model, constructed taking into account the mosquito life cycle and the main features of the Wolbachia -infection, we derived a delay diﬀerential model using the method of characteristics, to study the colonization and persistence of the Wolbachia - transinfected Aedes aegypti mosquito in an environment where the uninfected wild mosquito population is already established. Under some conditions, the model can be reduced to a Nicholson-type delay diﬀerential system; here, the delay represents the duration of mosquito immature phase that comprises egg, larva and pupa. In addition to mortality and oviposition rates characteristic of the life cycle of the mosquito, other biological features such as cytoplasmic incompatibility, bacterial inheritance


Introduction
Aedes aegypti is a widespread human blood-feeding mosquito responsible for the transmission of several arboviruses including Dengue, Yellow fever, Zika, Murray Valley, La Crosse, Chikungunya and Rift Valley fever.For most of these diseases an efficient vaccine is not available and the reduction of mosquito population still be the only way to prevent epidemics [17].The traditional approach to diminishing the mosquito population includes the reduction of breeding sites and the use of larvicides and pesticides for adults.In general, mechanical control and the application of larvicides are carried out before the period favorable to the proliferation of mosquitoes, while pesticides for adults are applied during epidemics when the number of infected humans is high [41].Environmentally-friendly techniques include the use of sterile males (SIT) [8,16] and Wolbachiainfected mosquitoes [7,15].While the first one focus on the reduction of mosquito population to halt disease transmission, the second one aims to replace the wild population for an infected one that is not able to transmit the virus.Both require the release of a large number of mosquitoes, hence a combination between traditional and new technologies are encouraged [33].
The intracellular bacteria of the genus Wolbachia manipulates host reproductive systems to increase its transmission by inducing parthenogenesis, feminization, male-killing or cytoplasmic incompatibility (CI) [2,37].Additionally, the bacteria is transmitted vertically from mother to its offspring.Put together, these characteristics confer a fitness advantage over uninfected population that can drive the Wolbachia-infected population to fixation [20,24].The same is not true for the SIT, since the mating with sterile mosquitoes does not produce viable offspring; they must therefore be introduced periodically.
In field, the Wolbachia strains that have been used in releases come from Drosophila melanogaster, wMel and wMelPop, and from Aedes albipicus, wAlbB.Artificial infections with new strains of bacteria still be done in laboratory in order to increase technique factibility.This is because recent research has shown that biotics and abiotics factors can influence Wolbachia densities and its distribution in mosquito tissue, and if the thresholds related to heritability and cytoplasmic incompatibility cannot be achieved the technique efficacy is lost [9,35,37].The thermal sensitivity of Wolbachia-infection is variable and can differ considerably between host species and strains.High temperatures might reduce its density in hosts, weaken the reproductive effects induced by Wolbachia-infection and even eradicate Wolbachia completely.In Ae. aegypti, the w AlbB infection type is more stable than w Mel and w MelPop at high temperatures [35].
Moreover, for poikilothermic species such as the Aedes aegypti mosquito, the body temperature depends on external factors and has strong effects over its entomological parameters and behavior.The lower and upper developmental threshold are 16 • C and 34 • C, being the development time shorter at higher temperatures [34].Also, the survival of immatures and adults may be negatively influenced by large diurnal temperature range since their mortality rates present U-shaped forms [42].The oviposition rate in turn increases quasi-linearly with temperature increasing [42].Both flight activity and mating rate were detected to increase with temperature ranging from 18 • C to 31 • C [11].Moreover, it was found that the length of the gonotrophic cycle was reduced with increasing mean temperatures.Besides, wing beat frequency, blood-feeding, biting activity, host seek among other behavior characters are significantly affected by temperature variation [34].
A lot of mathematical models have been addressing the use of Wolbachia-infected mosquitoes to control Dengue (and other viruses) transmission, because the presence of the bacteria reduces vector competence [28,31,32].In [14], a sex-structured model taking into account cytoplasmic incompatibility, male killing, incomplete maternal transmission, and different mortality rates for uninfected/infected population was developed.The boundedness of population was provided by considering competition among females for nesting places which give an upper limit for egg-laying rate.The ordinary differential model was studied analytically, and it was shown that the steady state where the Wolbachia-infected individuals dominate the population is possible when the maternal transmission is complete and cytoplasmic incompatibility is high.Coexistence of Wolbachia-infected and uninfected mosquito and Wolbachia-free equilibrium are found for a large set of relevant biological parameters.By considering that density-dependent death rate controls the exponential growth of populations [27], showed that only when the initial level of infection (given by the percentage of Wolbachia-infected population), breaks some critical thresholds that the infection takes off from the population (i.e. the threshold for invasion is achieved).In [31], the aquatic stage was also included and population boundedness was guaranteed by considering a logistic carrying capacity on this phase.It was shown that Wolbachia-infected mosquitoes always dominate the population provided they persist.The same approach was done in [29] and the existence of a minimum infection frequency above which Wolbachia could spread into the whole population of mosquitoes was explored.All of these mathematical models used ordinary differential equations to model the temporal dynamics of the mosquito population.
In turn, partial differential equations (PDE) are much less explored.Some studies such as [19,21] present reaction-diffusion models for the Wolbachia-infected and uninfected populations.They concluded that there is no spatial influence on the stability criteria for the steady states.Moreover [21], focused on determining the threshold for invasion of the wild population by the Wolbachia-infected one.Further [13], compares the stability results of the equilibriums obtained for an age-structured (PDE) with the one of an unstructured (ODE) model.For simplicity, two asexual population were considered, uninfected and Wolbachia-infected one.
Finally, rarer is the use of delay differential equations (DDE) for such problem.In [22], a DDE phasestructured model (larva and adult populations) evaluated the suppression of the wild population of Aedes mosquitoes by releasing a continuous constant number of Wolbachia-infected male.This is modeled by changing the growing rate of the population.The model considered two delays, one representing the average time from adult emergence to the hatching of the first larval stage (which determines larva population growth), and the other the average time from the first larval stage to adult emergence (which models adult population growth).Also, a strong density-dependent death rate was considered in the larval stage.They concluded that the delays do not impact population suppression.A modification of this model was proposed in [23] to compare two opposite phenomena which are the decrease on mating competitiveness of the released males relative to the wild males and the fitness advantage given by the cytoplasmic incompatibility probability to the Wolbachia-infected population over the wild one.The model considers only adult population and the delay gives the contribution of the last generation to the growth of the new population.They showed that CI plays a more important role in the suppression of Aedes population.
Here, starting from an age structured PDE model that considers the mosquito entomological parameters and also biological features associated to Wolbachia infection, a new two-population DDE model is obtained and carefully assessed.Analytical results such as positiveness, boundedness, and uniqueness of solutions are provided.Thresholds for existence and stability of the steady states were obtained and interpreted in the context of population fitness.The role played by the delay on the insect temporal dynamics can help to understanding the effect of changing on abiotic factors such as the temperature on the long-run behavior of this population.The model appears for the first time in [15] where numerical results concerning population dynamics were obtained.

Age structured partial differential model
Let w and u be Wolbachia-infected and uninfected mosquito status.We denote f j := f j (t, a), with j ∈ {w, u} the female population density of mosquito, a ∈ [0, τ ) the physiological age of the immature phase including egg, larva and pupa, a ≥ τ the physiological age of the mature phase (fertile adults), t the calendar time, µ w and µ u the adult mortality rates, and µ the immature mortality rate.We assume that the parameters τ and µ are the same for infected and uninfected mosquitoes [31].The temporal evolution of mosquito population satisfies the following age-structured Lotka-McKendrick system where The Wolbachia bacteria is transmitted from mother to its offspring with probability ξ w ∈ (0, 1).Thus, with probability an infected female can produce uninfected offspring.The mating probability between an uninfected female and an infected male is denoted by ν ∈ (0, 1) and the probability of cytoplasmic incompatibility occurrence is q ∈ (0, 1).This means that the fraction of matings between uninfected females and infected males that produce viable eggs is given by 1 − qν.The average birth rates are φ w > 0, φ u > 0 with the average percentage of female births r w , r u ∈ (0, 1).We denote by the total population of infected and uninfected adult females, respectively.Therefore, the newborn individuals introduced into the population are given at a = 0 by where G(X, Y ) = e −α(EwX+EuY ) η measures competition among individuals.More precisely, to take into account the competition between mosquitoes for oviposition sites, the number of eggs per female is multiplied by the density-dependent factor G (F w (t), F u (t)).The parameters α > 0 and η > 0 are, respectively, the environmental carrying capacity and the measurement of how rapidly it is achieved [12].The parameters E w and E u take into account the different behaviors between infected and uninfected females.In addition, the populations are assumed to satisfy lim a→+∞ f j (t, a) = 0, j ∈ {w, u}, t > 0. (2.5) The initial age-distribution f j (0, a), j ∈ {w, u}, is assumed to be known.Finally, as we are assuming that the mating between male and female is given by a constant parameter ν, we do not need to explicitly model the male population, and we will omit it from the analysis.

Reduction to a delay differential system
Henceforth, we reduce the system (2.1)-(2.5) to delay differential equations.We denote by the total population of immature females, Wolbachia-infected and uninfected, respectively.By integrating the system (2.1) over the age variable from 0 to τ and from τ to +∞, respectively, we get for j ∈ {w, u}, where F j is given by (2.3) and F i j by (3.1).On the other hand, the method of characteristics (see [36]) implies that As we are interested on the asymptotic behavior of the population, we can assume that t is large enough such that t > τ .Then, By adding the boundary conditions (2.4), we get the following delay differential system where F i j , F j , j ∈ {w, u}, are the total population of immature and adult females, Wolbachia-infected and uninfected, respectively.We can see that the equations of mature population F j are independent on the equations of immature one F i j .Then, we will omit the system of F i j and concentrate only on F j .Remembering that the nonlinear function G is given by we carry out the transformations and we define the new parameters respectively, the number of infected eggs per time per infected individual that will hatch, the number of uninfected eggs per time per uninfected individual that will hatch, and the number of uninfected eggs per time per infected individual that will hatch.Then, the model can be reduced, for t > τ , to with initial conditions given by We make a translation in time so as to define the system (3.3) on the interval [0, +∞) and the initial conditions where N = w + u and p = P u e −µτ .The equation (3.5) has been extensively studied in the literature [5,6,10].
The main results on equation (3.5) deal with the global attractivity of the positive steady state and the existence of oscillatory solutions (see [6]).

Positivity and boundedness of solutions
The positivity and boundedness of solutions are important in biological models.We first establish an existence and uniqueness theorem about the positive solution for the nonlinear delay differential system (3.3)-(3.4).
Then, from the system (3.3),we can write the following estimation .
Hence, y is bounded on the interval [−τ, t 0 ).This gives a contradiction and proves that the problem (3. .

Existence of steady states
Let (w * , u * ) be a steady state of the system (3.3).Then, (5.1) Defining and the system (5.1) can be rewritten as δ w e −µτ w * e −(w * +u * ) η − w * = 0, (δ u u * + δ wu w * ) e −µτ e −(w * +u * ) η − u * = 0, and we obtain three solutions S 0 , S u and S wu of this system: (i) Extinction of both populations (trivial equilibrium) (ii) Extinction of infected population and persistence of uninfected one , with R u = δ u e −µτ ; (5.5) (iii) Persistence of both populations (coexistence of infected and uninfected mosquitoes) with R w = δ w e −µτ and β wu = δ wu δ w − δ u + δ wu . (5.6) By examining the components of the steady states, we can deduce the following existence conditions.In terms of the original parameters, we have The two dimensionless parameters R w and R u are, respectively, the mean number of female infected offspring produced by a Wolbachia-infected female mosquito during her whole life, and the mean number of female uninfected offspring produced by an uninfected female mosquito during her whole life.
As we are interested on the relationship between temperature variation (that affects strongly the maturation time τ ) and population dynamics of both Wolbachia-infected and uninfected mosquitoes, we have to study the existence of the steady states in terms of the delay τ .We consider the following thresholds of the maturation time In fact, τ j ∈ R and we have: 1. τ j ≥ 0 if and only if δ j ≥ 1, and 2. τ j < 0 if and only if 0 < δ j < 1.
In terms of the delay, we obtain the following result.
1. S u exists if and only if 2. S wu exists if and only if 0 ≤ τ < τ w and τ u < τ w .
We remark that the steady states S u and S wu exist in the same time if and only if 0 ≤ τ < τ u and τ u < τ w .
In summary, we can distinguish four situations.
Then, for all τ ≥ 0, S 0 is the only steady state.(ii) Assume that τ u > 0 and τ w ≤ τ u .Then, (a) if 0 ≤ τ < τ u , there are two steady states S 0 and S u , with Components X and Y of the steady states S j = (X, Y ), with j = {0, u, wu} versus the delay τ .The parameter sets used in the simulations are given in Table 1, and in all cases η = 3.0.In (a), we plotted case (ii) in which τ u > 0 and τ w ≤ τ u .Then, if 0 ≤ τ < τ u , there are two steady states S 0 and S u , and if τ ≥ τ u , S 0 is the only steady state.In (b), we plotted case (iii) in which τ u ≤ 0 < τ w .Then, if 0 ≤ τ < τ w there are two steady states S 0 and S wu , and if τ ≥ τ w , S 0 is the only steady state.In (c), we plotted case (iv) in which 0 < τ u < τ w .Then, if 0 ≤ τ < τ u , there are three steady states S 0 , S u and S wu , if τ u ≤ τ < τ w , there are two steady states S 0 and S wu , and if τ ≥ τ w , S 0 is the only steady state.In all panels S 0 component Y appears as a dotted line, S u component Y appears as a solid, and S wu component X and Y appear, respectively, as a dot-dashed and dashed line.(c) if τ ≥ τ w , S 0 is the only steady state.
Figure 1 summarizes the results obtained in Proposition 5.3.The different scenarios correspond to the three parameter sets shown in Table 1.For each τ the corresponding steady state was obtained from (5.5) or (5.6).The panel (a) corresponds to case (ii), the panel (b) to case (iii), and the panel (c) to case (iv).We prove in Theorem 6.2 that the local asymptotic stability of the trivial steady state S 0 is given by the condition (6.1).We will need the following useful lemma.

Stability analysis of the steady states
Lemma 6.1.[Theorem A.5 in [26], p. 416, see also [18]] All roots of the algebraic equation λ + a + be −λτ = 0, have negative real parts if and only if Theorem 6.2.
1.If (6.1) is satisfied, then the trivial steady state S 0 is locally asymptotically stable.2. If (6.1) is not satisfied, then the trivial steady state S 0 is unstable.
Proof.We linearize the system (3.3)around S 0 .Then, we obtain The objective is to find conditions such that the solutions of (6.2) have negative real parts.We have to analyze separately the solutions of each factor of the product given by (6.2).For the equation λ + µ j − P j e −µτ e −λτ = 0, j ∈ {u, w}, the statement (i) aτ = µ j τ > −1 of Lemma 6.1, is always satisfied.On the other hand, the condition (ii) a + b = µ j − P j e −µτ > 0 of Lemma 6.1, is equivalent to τ > τ j := 1 µ ln (δ j ) , j ∈ {u, w}.
To prove our result (Thm.6.4) on the global asymptotic stability of the trivial steady state, we use the following lemma.
Lemma 6.3.( [6,25]).Suppose that there exists γ > 0 such that, for all y ∈ R 2 + , where . is a norm in R 2 and θ(A) is the matrix measure of A, θ(A) = lim Then, the trivial equilibrium of (6.3) is globally asymptotically stable.
Then, the trivial steady state S 0 of the system (3.3) is globally asymptotically stable.
We also have Finally, provided that 0 < e −µτ max{P w , P u + P wu } < min{µ w , µ u }, (6.4) Lemma 6.3 guarantees the global asymptotic stability of the trivial steady state S 0 .In fact, the condition (6.4) is equivalent to This finishes the proof of the theorem.
The local asymptotic stability of the other steady states will be analyzed by increasing the delay τ from zero with the possibility of eigenvalues to cross on the imaginary axis and the appearance of Hopf bifurcation.

Local asymptotic stability and Hopf bifurcation of Wolbachia-free steady state S u
The Wolbachia-free steady state S u := 0, (ln R u ) The linearization of the system (3.3)around the steady state S u is given by (6.5) Note that P u = µ u δ u , P w = µ w δ w and P wu = µ u δ wu .Then, the characteristic equation associate to (6.5) is given by The roots of the first term of the characteristic equation (see Lem. The statements (i) and (iii) are always satisfied and the statement (ii) is equivalent to Then, we can immediately conclude the following result.
Proposition 6.6.Suppose that Then, the steady state S u is unstable.

Now suppose that
Then, the local asymptotic stability of the steady state S u is given by the roots of the second term of the characteristic equation (6.6): Thanks to Lemma 6.1, the roots of (6.7) have negative real parts if and only if The statement (i) is always satisfied and the statement (ii) is equivalent to the condition that gives the existence of the steady state S u .Suppose that Then, the statement (iii) is satisfied and we have the local asymptotic stability of the steady state S u .In fact, the condition (6.8) is equivalent to max 0, τ u − 1 ηµ < τ < τ u .(6.9) We directly conclude the following result.
Proposition 6.7.If the condition (6.9) is satisfied then, the steady state S u is locally asymptotically stable.In particular, if then, for all 0 ≤ τ < τ u , S u is locally asymptotically stable.
This inequality is equivalent to We proved that S u is locally asymptotically stable for τ u < τ < τ u .Suppose that 0 ≤ τ < τ u .
When τ = 0, the characteristic equation (6.7) reads It has only one root As δ u > 1, then λ 0 < 0. We conclude that the steady state S u is locally asymptotically stable for τ = 0.By using a continuity argument, it is straightforward that there exists ∈ (0, τ u ), such that S u is locally asymptotically stable for τ ∈ [0, ).Consequently, when τ ∈ [0, τ u ) increases, the stability of S u can only be lost if characteristic roots cross on the imaginary axis.We look for purely imaginary roots ±iω, ω ∈ R. Remark that if λ is a characteristic root then its conjugate λ is also a characteristic root.Then, we look for purely imaginary roots iω with ω > 0. By separating real and imaginary parts in the characteristic equation (6.7), we get Adding the squares of both hand sides of the last system and using the fact that cos 2 (τ ω) + sin 2 (τ ω) = 1, it follows that For the existence of ω > 0, it is necessary to have Then, it is immediate to conclude the following result.
then, for all 0 ≤ τ < τ u , S u is locally asymptotically stable.

Now suppose that
and consider the function : [0, τ u ) → (0, +∞) defined by In fact, we have Then, for each τ ∈ [0, τ u ), there is a unique solution Θ(τ ) ∈ [0, 2π) of the system Then, Θ(τ ) ∈ (π/2, π) and it is given by We conclude that the system (6.10) is equivalent to find τ ∈ [0, τ u ) solution of with (τ ) given by (6.11) and Θ(τ ) by (6.12).We remark here that in all this study, the set N includes 0. This is equivalent to solve More precisely, we have to solve for k ∈ N and τ ∈ [0, τ u ), The functions Z k (τ ) are given explicitly.However, we cannot determine explicitly their roots.The roots can be found numerically.The following lemma states some properties of the functions Z k , k ∈ N. Lemma 6.9.For all k ∈ N and τ ∈ [0, τ u ), Therefore, provided that no root of Z k is a local extremum, the number of positive roots of Z k , k ∈ N, on the interval [0, τ u ) is even.
This lemma implies, in particular, that, if Z k has no root on [0, τ u ), then no function Z j , with j > k, has roots on [0, τ u ).The next proposition is a direct consequence of Lemma 6.9.Proposition 6.10.If the function Z 0 defined on the interval [0, τ u ), by has no root, then the steady state S u is locally asymptotically stable for all τ ∈ [0, τ u ).
We now suppose that Z 0 , under the condition η > 2 ln(δ u ) , has at least one positive root on the interval [0, τ u ).Let τ * u ∈ (0, τ u ) be the smallest root of Z 0 .Then, S u is locally asymptotically stable for τ ∈ [0, τ * u ), and loses its stability when τ = τ * u .A finite number of stability switch may occurs as τ increases and passes through roots of the Z k functions.
Our next objective is to prove that S u can be destabilized through a Hopf bifurcation as τ ∈ [0, τ u ) increases.We start by proving that if an imaginary characteristic root iω exists then, it is simple.Suppose, by contradiction, that λ = iω is not a simple characteristic root.Then, λ is a solution of ∆(τ, λ) = 0 and This is equivalent to The two equations of the system (6.16)lead to This a contradiction with the fact that λ = iω.As τ * u is the smallest root of Z 0 then, from the definition of Z 0 , the characteristic equation (6.15) has purely imaginary roots ±i (τ * u ), where is defined by (6.11).The stability of the positive steady state switches from stable to unstable as τ passes through τ * u .Other stability switch occur when τ passes through roots of the Z k functions (see [4]).Now, we rewrite the characteristic equation (6.15) in the following form ∆(τ, λ) := A(τ, λ) + B(τ )e −λτ = 0.
We define, for λ = iω, the polynomial function Then, Let λ(τ ) be a branch of roots of (6.15) such that λ(τ * u ) = i (τ * u ).The Hopf bifurcation theorem says that a Hopf bifurcation occurs at S u when τ = τ * u if sign d (λ(τ )) dτ We know from [4] that That is to say It is clear that It follows The following proposition states the existence of a Hopf bifurcation at τ = τ * u that destabilizes the positive steady state S u .Proposition 6.11.If Z 0 (τ ) has at least one positive root on the interval (0, τ u ), then the positive steady state S u is locally asymptotically stable for τ ∈ [0, τ * u ), where τ * u is the smallest root of Z 0 (τ ) on (0, τ u ), and S u loses its stability when τ = τ * u .A finite number of stability switch may occur as τ passes through roots of the Z k functions.Moreover, if then a Hopf bifurcation occurs at S u for τ = τ * u .Figure 2 shows, for two set of parameters, the existence or non-existence of roots for the functions Z 0 and Z 1 , given by the equation (6.13).In each case, τ ∈ (0, τ u ).On the left, we have η = 2.0 and no root for Z k .Then, the equilibrium S u stays locally asymptotically stable on the interval (0, τ u ).On the right, we have η = 3.0 and two roots for Z 0 (no root for Z 1 ).A Hopf bifurcation occurs at τ = τ * u and periodic oscillations around the equilibrium S u are observed until the threshold τ = τ + u .The equilibrium S u corresponds to the extinction of Wolbachia-infected mosquito and the persistence of uninfected one.
For the set of parameters given in the case (ii), Table 1, where the equilibrium S u exists for τ ∈ [0, τ u ), we can see in Figure 3, for each value of η, the roots of Z 0 .For η less than a threshold, η min (η min ≈ 2), Z 0 has no Figure 2. The functions Z k , k = 0, 1, given by the equation (6.13), versus τ .On the left, we can see that Z 0 and Z 1 have no roots.The vertical line shows the right end of Z k domain, τ u = 13.65.On the right, we can see that Z 0 has two roots τ * u = 4.01 and τ + u = 12.67, and Z 1 has no root.The vertical line shows the right end of Z k domain τ u = 15.98.At τ = τ * u a Hopf bifurcation occurs and periodic oscillations around the equilibrium are observed until τ = τ + u ; outside the interval (τ * u , τ + u ), the equilibrium S u is locally asymptotically stable.In both panels, we use the parameter set from case (ii) (Tab.1).On the left, we set η = 2.0, and on the right, we set η = 3.0.Figure 3. η versus τ .The continuous curve provides the roots of Z 0 (τ ) given η (see 6.14).For the set of parameters chosen, case (ii) (Tab.1), the horizontal dotted line highlights the lower threshold to have Hopf bifurcation.The domain of Z 0 is given by (0, τ u ) with τ u = 17.76.As an example, the dashed line is set for η = 3.0.For this value of η, the corresponding roots of Z 0 are τ * u = 3.85 and τ + u = 14.67.As dZ 0 (τ * u )/dτ > 0 a Hopf bifurcation occurs at τ = τ * u by Proposition 6.11.At this point, the equilibrium S u looses his stability and periodic oscillations can be seen until τ + u .At τ = τ + u the stability of S u is recovered.This corresponds to the extinction of Wolbachia-infected mosquito and persistence of uninfected one.root which implies that the steady state S u is stable.For η greater than η min , Z 0 has two roots τ * u and τ + u .In this case, S u looses stability at τ = τ * u and periodic oscillations can be seen.For τ = τ + u the stability of S u is recovered.The equilibrium S u corresponds to the extinction of Wolbachia-infected mosquito and persistence of uninfected one.u , τ + u ); and (c) the temporal evolution of u(t) (solid line) and w(t) (dot-dash line) given by the system (3.3) with τ = 10.In all panels, the parameters were taken from Table 1 case (iv) and η = 3.The thresholds given by τ * u = 3.85, τ + u = 14.67, τ u = 17.76 and τ u = 22.95 are, respectively, the first and second roots of Z 0 , the right end of the domain of the function Z 0 , and the right end of τ that allows the existence of S u .
Using the set of parameters from case (ii), Table 1, and η = 3.0, in Figure 4, we can see in the panel (a) the minimum and the maximum values of the periodic solutions u(t) of the system (3.3)plotted for τ * u ≤ τ < τ + u (where S u is unstable).This corresponds to the amplitude of u(t) for τ * u ≤ τ < τ + u .In the panel (b), we plotted the period of these periodic oscillations which is an increasing function of τ ∈ (τ * u , τ + u ).Finally in the panel (c), we can see an example of temporal evolution of the system (3.3) for η = 3.0 and τ = 10.The component u(t) oscillates around the steady state and the component w(t) tends to zero with positive damped oscillations.The other parameters were taken from the case (iv) of Table 1.A Hopf bifurcation occurs at τ ≈ 3.5 and the stability of S wu is restored at τ ≈ 8.2; τ ∈ (0, 22.33).
and interesting questions, such as the equilibriums and their local and global stability, can be analytically addressed while keeping all the biological assumptions behind the model.Thus, the positivity, boundedness and uniqueness of solutions were proved.The model admits three steady states: extinction of both populations (S 0 ), extinction of infected population and persistence of uninfected population (S u ), and persistence of both populations (S wu ).The conditions of existence of each equilibrium were established as a combination of the entomological parameters that describe mosquito life's cycle (such as mortality, oviposition and development rates) and the effects of the Wolbachia presence in the host (vertical transmission of the bacteria, cytoplasmic incompatibility, and sex-ratio-distorting).The two thresholds R w and R u can be interpreted as the net reproductive rates which are defined as the average numbers of female offspring that a female produces during her lifetime.Therefore, R w measured the number of infected offspring produced by an infected female and R u measured the number of uninfected offspring produced by an uninfected female.Moreover, we could prove that S 0 always exists, S u exists if and only if R u > 1 and S wu exists if and only if R w > max{1, R u }.
As we were interested in analyzing the effect of variation on developmental time in population dynamics, we rewrote the condition of existence of each equilibrium in terms of the delay τ which measures the time spent from egg to adult.Four scenarios can be drawn (Fig. 1): (i) assume that τ u , τ w ≤ 0 (these thresholds are equivalent to the ones on R u and R w , respectively).Then, for all τ ≥ 0, S 0 is the only steady state; (ii) assume that τ u > 0 and τ w ≤ τ u .Then, if 0 ≤ τ < τ u , there are two steady states S 0 and S u , and if τ ≥ τ u , S 0 is the only steady state; assume that τ u ≤ 0 < τ w .Then, if 0 ≤ τ < τ w , there are two steady states S 0 and S wu , if

Proposition 5 . 1 .
(a) S 0 always exists; (b) S u exists if and only if R u > 1; (c) S wu exists if and only if R w > max{1, R u }.

Figure 4 .
Figure 4. From the left to the right we have: (a) the minimum and the maximum values of the component u(t) of the system (3.3)plotted against τ ∈ (τ * u , τ + u ); (b) the period of the oscillations of u(t) versus τ ∈ (τ *u , τ + u ); and (c) the temporal evolution of u(t) (solid line) and w(t) (dot-dash line) given by the system (3.3) with τ = 10.In all panels, the parameters were taken from Table1case (iv) and η = 3.The thresholds given by τ * u = 3.85, τ + u = 14.67, τ u = 17.76 and τ u = 22.95 are, respectively, the first and second roots of Z 0 , the right end of the domain of the function Z 0 , and the right end of τ that allows the existence of S u .