DYNAMICS OF DELAYED MOSQUITOES POPULATIONS MODELS WITH TWO DIFFERENT STRATEGIES OF RELEASING STERILE MOSQUITOES

. To prevent the transmissions of mosquito-borne diseases (e.g., malaria, dengue fever), recent works have considered the problem of using the sterile insect technique to reduce or eradicate the wild mosquito population. It is important to consider how reproductive advantage of the wild mosquito pop- ulation oﬀsets the success of population replacement. In this work, we explore the interactive dynamics of the wild and sterile mosquitoes by incorporating the delay in terms of the growth stage of the wild mosquitoes. We analyze (both analytically and numerically) the role of time delay in two diﬀerent ways of releasing sterile mosquitoes. Our results demonstrate that in the case of constant release rate, the delay does not aﬀect the dynamics of the system and every solution of the system approaches to an equilibrium point; while in the case of the release rate proportional to the wild mosquito populations, the delay has a large eﬀect on the dynamics of the system, namely, for some parameter ranges, when the delay is small, every solution of the system approaches to an equilibrium point; but as the delay increases, the solutions of the system exhibit oscillatory behavior via Hopf bifurcations. Numerical examples and bifurcation diagrams are also given to demonstrate rich dynamical features of the model in the latter release case.


(Communicated by Xingfu Zou)
Abstract. To prevent the transmissions of mosquito-borne diseases (e.g., malaria, dengue fever), recent works have considered the problem of using the sterile insect technique to reduce or eradicate the wild mosquito population. It is important to consider how reproductive advantage of the wild mosquito population offsets the success of population replacement. In this work, we explore the interactive dynamics of the wild and sterile mosquitoes by incorporating the delay in terms of the growth stage of the wild mosquitoes. We analyze (both analytically and numerically) the role of time delay in two different ways of releasing sterile mosquitoes. Our results demonstrate that in the case of constant release rate, the delay does not affect the dynamics of the system and every solution of the system approaches to an equilibrium point; while in the case of the release rate proportional to the wild mosquito populations, the delay has a large effect on the dynamics of the system, namely, for some parameter ranges, when the delay is small, every solution of the system approaches to an equilibrium point; but as the delay increases, the solutions of the system exhibit oscillatory behavior via Hopf bifurcations. Numerical examples and bifurcation diagrams are also given to demonstrate rich dynamical features of the model in the latter release case.
1. Introduction. For more than a century, human beings have attempted to control blood-feeding mosquitoes. This is because of the significant mortality and morbidity burden associated with mosquito-borne diseases (e.g., malaria, dengue fever, and West Nile virus), which are transmitted between humans via bloodfeeding mosquitoes [7,38]. Various control approaches have been explored, which include the development of more effective drug treatments, vaccines, and vector (mosquito) suppression. Vector control measures include the elimination or reduction of their nesting places by draining stagnant water deposits. This measure has been very effective in places where the prevalence of the disease was not very high. Indoor and outdoor insecticide spraying has been applied for many years for controlling mosquito populations. Despite some mosquito-borne diseases have been successfully controlled in many regions through vector-targeted intervention such as insecticide-treated bed nets (ITNs) and indoor residual sprays (IRS), massive and long time spraying of adulticide is not recommended. Since they have commonly been chemically-based, the effectiveness of this measure has been hampered by the appearance of insecticide resistant vector strains [1,2,8]. The genetically-altered or transgenic mosquitoes may provide a new and potentially effective weapon to fight these mosquito-borne diseases [14,20,29]. Sterile Insect Technique (SIT) is an environmentally friendly alternative strategy that is gaining renewed interest for the control of mosquito populations. The technique involves the massive release of male mosquitoes (sterilized through radiological or chemical means) into the environment to mate with wild mosquitoes that are present in the environment. Female mosquitoes mating successfully with a sterile male will either not reproduce, or produce eggs which will not hatch. Repetitive releases of sterile mosquitoes or releasing a significantly large number of sterile mosquitoes may eventually wipe out or suppress a wild mosquito population [5,37,3].
Mathematical models (see [11,13,17,22,18,34,25,9,26,40,41] and the references therein) have been formulated to explore the interactive dynamics of wild and sterile mosquito populations, and potentially estimate the effectiveness of the control strategy of mosquitoes. In particular, the models in [9,26] are ODE systems under the homogeneous assumption for the wild and sterile mosquitoes. However, these models hardly take into account the stage structured life-history of the mosquitoes, which can have significant effects on their dynamics (see [17,10,27,39]). In fact, all mosquito species go through four distinct stages (egg, pupa, larva, and adult) during a lifetime. The first three stages occur in water, but the adult is an active flying insect. Only the female mosquito bites and feeds on the blood of humans or other animals. Murdoch et al. [31] have shown short-period population oscillations in abundance occur from the developmental lags between mosquito life-history stages. In recent years, several interesting mathematical models have been developed to investigate the dynamics of mosquitoes with stage structure, for instance, the discrete models in [28], the continuous time models in [16,27] and the delayed models in [15,24].
Based on the models in [9,26,32,15], in this paper, we assume that the wild mosquito population growth is stage-structured. Let w(a, t) denote the density of wild mosquitoes at time t of age a. Let τ be the maturation time for all mosquitoes (the total time from egg to adult). During the larvae maturation, they are subject only to the possibility of natural death. Let µ 0 be the per-capita natural mortality. We have ∂w(a, t) ∂t + ∂w(a, t) ∂a = −µ 0 w(a, t), 0 < a < τ, t > 0.
Let W (t) = ∞ τ w(a, t)da be the wild adult mosquitoes population at time t ≥ 0, µ 1 the per-capita natural death rate of adult mosquitoes, and ξ 1 the density-dependent death rate of adult mosquitoes. Let g(t) be the sterile mosquitoes population at time t ≥ 0, µ 2 and ξ 2 the density-independent and dependent death rates of sterile mosquitoes. The interactive dynamics of wild adult mosquitoes and sterile DYNAMICS OF DELAYED MOSQUITOES POPULATIONS MODELS 1183 mosquitoes are governed by the following equations: where B(·) is the release rate of the sterile mosquitoes. Assume that w(∞, t) = 0. Integrating the first equation of the above system over the interval [τ, ∞) with respect to a gives Use the method of characteristics to solve (1) and get w(a, t) = w(0, t − a)e −µ0a for 0 ≤ a ≤ τ and t ≥ a. In particular, w(τ, t) = w(0, t − τ )e −µ0τ for t ≥ τ . Supposed that the interactions between the two types of mosquitoes lead to egg-laying rate given by where a 0 is the number of wild offspring produced by per mating between the wild mosquitoes. Then it follows that Substituting this formula into (2) and replacing W (t) by w(t) (for easing the notation) leads to the delayed model of the adult wild and sterile mosquitoes: The factor e −µ0τ in (4) stands for the survival rate of the immature mosquitoes who were born at time t − τ and still remain alive at the time t. The initial conditions for system (4) take the following form where t 0 ≥ τ , φ(θ) and ψ(θ) are positive continuous functions for θ ∈ [−τ, 0]. The main aim of this paper is to investigate how the delay τ affects the dynamics of system (4) with different strategies of releasing rate of sterile mosquitoes. In Section 2, we first explore the model (4) with the constant release rate of sterile mosquitoes and analyze the effect of the delay τ on the dynamics of mosquito populations. Our results show that the delay does not affect the dynamics of the model. In Section 3, we investigate (4) with the release rate of sterile mosquitoes proportional to the wild mosquito population, and obtain the conditions on the stability of the positive equilibria and the occurrence of Hopf bifurcations for certain values of the delay. In Section 4, numerical examinations are provided to demonstrate the complexity of the model dynamics in the latter case. In Section 5, we provide a brief discussions on our findings, in particular on the impact of the time delay on the mosquito control measures.

2.
Constant release rate of the sterile mosquitoes. In this section, we investigate the model (4) with constant release rate of the sterile mosquitoes, i.e., we take B(·) = b with b > 0. Then (4) becomes The initial conditions of system (6) satisfy (5). We first establish the following result: Lemma 2.1. For any given initial data in (5) that are positive on [t 0 − τ, t 0 ] for some t 0 ≥ 0, system (6) has a unique solution (w(t), g(t)), which are defined, positive and bounded for all t ≥ t 0 − τ .
Proof. The local existence and uniqueness of the solution (w(t), g(t)) of (6) subject to a prescribed initial condition (5) follows from the standard theory for delay differential equations (e.g., Theorem 3.1 in [19]). The global existence of this solution in [t 0 − τ, +∞) follows from the positiveness and boundedness of this solution that we prove below. Let [t 0 − τ, T ) be the maximal interval of existence of (w(t), g(t)), and let We claim that (w(t), g(t)) ∈ Ω 1 := {(w, g) : 0 < w < K 1 , 0 < g < K 2 } for [t 0 − τ, T ). First, since w(t 0 ) > 0 and g(t 0 ) > 0, and for t ∈ [t 0 , T ) it follows that w(t) > 0 and g(t) > 0 for t ∈ [t 0 , T ). Now we prove that w(t) < K 1 . If not, there would exist the smallest t 1 > t 0 such that w(t 1 ) = K 1 , w (t 1 ) ≥ 0 and w(t) < K 1 for any time t 0 − τ ≤ t < t 1 . However, from the first equation of system (6) we have a contradiction. This shows that w(t) < K 1 for t ∈ [t 0 , T ). Similarly using the second equation of (6) we show that g(t) < K 2 for t ∈ [t 0 , T ).
We thus conclude the above claim. Since Ω 1 is bounded, it follows that T = ∞ and (w(t), g(t)) ∈ Ω 1 for t ≥ t 0 .
We note that the above proof implies that for any K 1 > a − µ 1 ξ 1 and K 2 > b µ 2 , the set Ω defined in the above proof is positively invariant for the flows of (6).
Since the presence of the delay does not change the number of the equilibrium solutions in system (6), the existence of the equilibria of (6) follows from the same argument as for ODE system in [9] . Let N =N (τ ) : and define the threshold release value of the sterile mosquitoes as b 0 (τ ) := 1 2a We have the following result similar to the one in the reference [9].
Theorem 2.2. Let τ ≥ 0 and b > 0. Then the system (6) has a unique boundary nonnegative equilibrium E 0 (0, g 0 ) where g 0 is given in (7). Furthermore, regarding the existence of positive equilibria, we have: does not have any positive equilibrium; given by: We now investigate the stability of the equilibria of system (6). Upon linearizing system (6) at the equilibrium E 0 , E * 1 , or E * 2 , we obtain the characteristic equation where
At the positive equilibria E * 1 and E * 2 (when 0 < b < b 0 ), by direct calculations we reduce the characteristic equation (8) as: Recall that for τ = 0, the system (6) becomes the system (2.2) in [9]. From Theorem 2.2 in [9], we know that the positive equilibrium E * 1 is a saddle, and E * 2 is a stable node or focus. The following results show that the stabilities of these equilibria remain the same for τ > 0.
is locally asymptotically stable. Proof. The stability of E 0 is already proved above. So we only need to prove (ii). We shall complete the proof by the following two steps.
Step 2. E * 2 is stable when τ > 0. To do so, we need the following Proposition (for the clearness of the proof we write L τ (λ) for L(λ) for τ ≥ 0): Proposition. For each τ > 0, L τ (λ) = 0 has countably many roots, all lying to the left of a vertical line Reλ = β for some real number β; moreover, for any α < β there are at most finitely many of these roots lying in the vertical strip α < Reλ < β.
Using this proposition, the Rouché's theorem [12], and the fact that for τ = 0, L 0 (λ) = 0 has only two roots with negative real parts, we conclude that for some sufficiently δ > 0 and any 0 < τ < δ, all the roots of L τ (λ) = 0 has negative real parts (NRP). For simplicity, we say L τ has the NRP-property if all its roots have negative real parts. Recall we are trying to show that for any τ > 0, L τ has the NRP-property. Assume by the contradiction that this is false. Then, letting we have that τ 0 ≥ δ > 0, and L τ0 has at least one root, say, λ(τ 0 ), such that Reλ(τ 0 ) ≥ 0. Using the above property, and the Rouché's theorem we can exclude the case that Reλ(τ 0 ) > 0, yielding Reλ(τ 0 ) = 0. Hence we may assume that λ(τ 0 ) = iω for some ω ≥ 0 and then insert it into (9) with τ = τ 0 to get Taking the absolute values of both sides yields From (10), we have c(τ 0 ) − d(τ 0 ) > 0. Note that by the same argument used To simplify the expressions in the following calculations, we let A 1 = µ 1 + ξ 1 N and A 2 = µ 2 + ξ 2 N. It follows from direct calculations that and so From (14) that B 1 = b(τ 0 ) 2 + 2c(τ 0 ) − a 2 (τ 0 ) < 0. Now using the above estimates for B 1 and c 2 (τ 0 ) − d 2 (τ 0 ), we conclude that (13) cannot hold for any real ω. This contradiction shows τ 0 = ∞, thereby completing the proof of Step 2 and the proof of Theorem 2.3.
In Figure 1 and Figure 2, we perform simulations to show that in system (6), the bistable phenomenon still occur and the time delay τ of stage-structure growth has only effect on the level of the positive equilibria with τ being varied.
3. Release rate proportional to wild mosquito population. In this section, we consider the system (4) with the release rate B(·) proportional to the wild mosquito population. As in the reference [9], we also incorporate the Allee effect [35] to account for the difficulty and stochasticity of finding mates when the populations of mosquitoes are small. Then the system (4) becomes: The initial conditions of system (15) satisfy (5). Note that when τ = 0, the system reduces to the system (3.1) in [9]. As in the previous section, we investigate how the delay τ affects the dynamics of (15) when varying τ . First, by a similar proof for Theorem 2.1, we have the following result: Theorem 3.1. For any given initial data in (5) that are positive on [t 0 − τ, t 0 ] for some t 0 ≥ 0, system (15) has a unique solution (w(t), g(t)), which are defined, positive and bounded for all t ≥ t 0 − τ .
Again use the fact that the number of the steady state solutions of the system (15) with τ > 0 are the same as that with τ = 0. Hence we can apply Theorem 3.1 in [9] (with a = a 0 e −µ0τ ) to obtain Theorem 3.2 below. Let and LetN be the point in (N 1 , N 2 ) such that G (N ) = 0. We then define the threshold release value of the sterile mosquitoes as b 0 (τ ) := G(N ).
(ii) If (H 1 ) and (H 2 ) hold for some τ * ∈ (0, τ 0 ), then the system (15) undergoes a Hopf bifurcations at E * 2 when τ = τ * . In theorem 3.1 (ii), the conditions (H 0 ) − (H 2 ) guarantee the occurrence of the Hopf bifurcation at some τ * . However, since we do not have the explicit formulas of E * 2 = (w * 2 , g * 2 ) for the general parameters in (15), these conditions are hard to be verified analytically. Below we show that the Hopf bifurcations do occur for some parameter ranges.
To this end, we consider the system (15) with the parameters satisfying:
Proof. We verify the assumptions in (28). Under the either assumption of the corollary we have ∆ 0 2 > 0. Since lim ξ2→0 ∆ 0 1 = 3 16 − b − 1 2 ξ 1 , we see that ∆ 0 2 > 0 under the either assumption. Finally, if the first assumption holds, then if the second assumption holds, then This shows that ∆ 0 > 0 under the either assumption. Applying Theorem 3.2 completes the proof of the corollary.

4.
Numerical simulation of the model (15). In this section, we will focus on the numerical simulation of the model (15). First, let the parameter values a 0 = 30, µ 0 = 0.1, µ 1 = 0.5, µ 2 = 1.5, b = 2, ξ 1 = 4, ξ 2 = 0.51 in system (15). In Figure  3, we show that the time delay τ have effect on stability of the positive equilibrium E * 2 and there is a stable periodic solution for τ = 4.9. Then, we choose to use the same set of parameters as in Figure 3 except τ and b. We let τ vary and use it as a bifurcation parameter. Other parameter values are fixed except parameter b in Figure 9, where b is varying as a secondary free parameter to generate the two dimensional bifurcation diagram. Using Matlab, we obtain the bifurcation diagram (see Figure 4).  From the bifurcation diagram in Figure 4, we can see that there is a stable steady state for delay τ from 0 up to around 4.5. There exists a stable periodic solution for delay between 4.5 to 5.2. Interestingly, there is a discontinuity at τ ≈ 5.2. Actually after further exploration, we find the discontinuity occurs since there exist two stable periodic solutions for delay τ close to 5.2. One of them has a smaller amplitude and the other one has a larger amplitude. Due to the change of their attraction basin as delay τ varies, a solution which originally approaches to the periodic solution with a smaller amplitude, switches to the periodic solution with a larger amplitude. When this switch occurs, we see a jump in our bifurcation diagram in Figure 4. There is another interesting observation on the bifurcation diagram around τ ≈ 6.8 where the bifurcation diagram stop abruptly. We did further exploration there too and found the amplitude of the periodic solutions shrink abruptly within a very small interval of delay and our Matlab code fails to work appropriately due to a precision issue.
To demonstrate the existence of two stable periodic solutions, we choose delay τ = 5.18 and choose two different sets of constant initial values (w, g) = (10, 5) and (w, g) = (0.5, 1.5). The solution associated with the initial values (w, g) = (10, 5), approaches to the periodic solution with a larger amplitude and is plotted in dotted line in Figure 5. The solution associated with the initial values (w, g) = (0.5, 1.5), approaches to the periodic solution with a smaller amplitude and is plotted in solid line in Figure 5.
Similarly there is a secondary discontinuity at delay τ ≈ 6.5. There also exist two stable periodic solutions for delay close to 6.5. To illustrate that, we choose (w, g) = (0.5, 1.5). Accordingly we obtain two stable periodic solutions which are plotted in Figure 6.
To have a better understanding about system (15), we use DDEbiftool to investigate the system further and obtain some bifurcation diagram. First we obtain a one dimensional bifurcation diagram using delay τ as a bifurcating parameter (see, Figure 7). In this figure we can see, the system has two positive equilibria which merge and disappear through a limit point bifurcation at around τ = 6.2. For τ between 4.7 and 7.2, there are three periodic solutions. Some of them is stable and some of them is unstable. We calculate the stability of periodic solutions along these curves and place the result in Figure 8.
In Figure 8, it is easy to see that there exists two stable periodic solutions for delay between τ ∈ (4.8, 5.8). This observation is consistent with our simulation result in Figure 4. The first branch of periodic solution is stable for τ ∈ (4.55, 5.5). The second branch periodic solution is stable for delay τ ∈ (4.8, 5.7), and the third branch of periodic solution is stable from delay τ ∈ (5.55, 7.2). This bifurcation diagram confirms the existence of two stable periodic solutions and clearly indicate delay values where one can expect the two stable periodic solutions coexist.
In addition, via tracing the periodic branches, we generate a two dimensional bifurcation diagram (see, Figure 9 ) in the parameter space of (τ, b). From the bifurcation diagram, we can see the system undergoes supercritical Hopf bifurcation (the red curve), subcritical Hopf bifurcation (the cyan curve which is very close to red curves and very hard to see), the fold-Hopf bifurcation (blue curve) and torus bifurcation (the black curve). In addition, along Hopf bifurcation curves, other more degenerated points are detected including Hopf-Hopf bifurcation (six black circles) and Generalized Hopf bifurcation (one black square). At those degenerated points, Hopf bifurcation curves, Fold-Hopf bifurcation curves or torus bifurcation curves meet each other. We believe more interesting dynamics can be expected around  these Hopf-Hopf points or Generalized Hopf point. For torus bifurcation, we also plot the profile of orbits on the torus at the bifurcation point and placed them in Figure 10. In summary, Figure 9 shows that very complicated dynamics can be expected from the model. In another follow up work, we will investigate the model further around the Hopf-Hopf bifurcation, Generalized-Hopf bifurcation points.
5. Discussion. In this work, we studied (both analytically and numerically) the dynamics of interactive systems of the wild and sterile mosquitoes with different releasing sterile mosquitoes strategies by incorporating the delay in the growth stage of wild mosquito populations. Our obtained results in Theorem 2.3 have shown that the growth time delay of the mosquito population has no destabilizing effect on the solution's behavior in the constant release rate of sterile mosquitoes. That is, in this case, for any size of the growth delay, the solutions of the system approach to either the boundary equilibrium or the stable positive equilibrium, depending on the size of the release rate. This implies that the control of the wild mosquito population is highly dependent on the rate at which the sterile mosquito are released, with only high release rates giving sufficient control. However, when the number of the sterile releases is s proportional to the wild mosquito population, analytical results in Theorems 3.1 and 3.2 and numerical results in Section 4 suggest that the solutions of the model exhibit complex dynamical behavior. Some phenomena are observed on fluctuating interactions between the wild and sterile mosquitoes. As the time  On the graph, the torus bifurcation curve is very close to Fold-Hopf bifurcation curve. To have a better view, we include a zoomed figure. delay varies, we have found Hopf bifurcations, Bautin bifurcations and Hopf-Hopf bifurcations in the system. The existence and stability of periodic solutions created in these three types of bifurcations are investigated by numerical analysis. These obtained results describe the equilibrium of system process. In particular, when a stable periodic orbit exists, it can be understood that the wild and sterile mosquitoes system can coexist for a long term although the wild mosquito is not eliminated. The conditions for the existence of the bifurcations indicate that the parameters of the system are important in controlling the development and progression of wild mosquitoes.
In summary, by investigating two different delay differential equations models, we have provided convincing evidences that a simple delay that accounts for such processes as the growth stages of mosquito populations is not alone responsible for sustained oscillations between the interactive dynamics of the wild and sterile mosquitoes. If such oscillations occur, then they must be the consequence of delayed dependence on other processes, or of a more complex dependence on the past population density and Allee effect [4,21]. Therefore, our research in this paper shows that the choice of the releasing sterile mosquito strategies is an important determinant of overall mosquito population dynamics.
Finally, we would like to stress that for many mosquito populations, seasonal environmental effects are an important feature of their life-history, giving rise to seasonal cyclic dynamics [36,30]. For example, seasonal rainfall can increase the abundance of mosquitoes, where the reproduction depends on the availability of suitable breeding sites such as water-filled containers. Moreover, many species of mosquitoes in temperate zones overwinter in a diapausal state [33]. Hence, in order to fully understand the influences of the releasing sterile mosquitoes on control wild mosquitoes, further investigation on seasonal environmental effects must be addressed in the modeling, which is planned in our future research.