A DELAYED DIFFERENTIAL EQUATION MODEL FOR MOSQUITO POPULATION SUPPRESSION WITH STERILE MOSQUITOES

. The technique of sterile mosquitoes plays an important role in the control of mosquito-borne diseases such as malaria, dengue, yellow fever, west Nile, and Zika. To explore the interactive dynamics between the wild and sterile mosquitoes, we formulate a delayed mosquito population suppression model with constant releases of sterile mosquitoes. Through the analysis of global dynamics of solutions of the model, we determine a threshold value of the release rate such that if the release threshold is exceeded, then the wild mosquito population will be eventually suppressed, whereas when the release rate is less than the threshold, the wild and sterile mosquitoes coexist and the model exhibits a complicated feature. We also obtain theoretical results including a suﬃcient and necessary condition for the global asymptotic stability of the zero solution. We provide numerical examples to demonstrate our results and give brief discussions about our ﬁndings.


1.
Introduction. Sterile insect technique (SIT) is a promising technique that has been first introduced by E. F. Knipling and his collaborators and experimented successfully in the early 50's by nearly eradicating screwworm population in Florida. Thereafter, SIT has been applied to control of different pests and vectors such as fruit flies and mosquitoes [10]. In recent years, the use of sterile mosquitoes has provided a bio-safety control method for eradicating or reducing mosquito populations and thus to control the mosquito-borne diseases, such as dengue, malaria, yellow fever, west Nile, and Zika. It is different from the traditional methods of killing mosquitoes by spraying insecticides. Utilizing radical or other chemical or physical techniques, sterile male mosquitoes are produced who are incapable of producing offsprings despite being still sexually active. These sterile male mosquitoes are then released into the field to mate with wild mosquitoes. A wild female mosquito that mates with a sterile male mosquito will either not reproduce, or produce eggs but the eggs will not hatch. Repeated releases of a significantly large number of sterile mosquitoes then can possibly wipe out a wild mosquito population eventually, or more realistically reduce the density of the wild mosquito population [1,30].
Mathematical models have been formulated to study the interactive dynamics of wild and sterile mosquito populations or the control of mosquitoes [6,7,9,11,13,14,15,16,17,18,19,23,24,34,35,36,38,39], among which the models in [6,23] involve three different strategies in the releases of sterile mosquitoes with constant release rate, the rate proportional to the wild mosquitoes, and the rate proportional to wild mosquitoes with saturation. It is assumed in [6] that the mosquito population follows the logistic growth in the absence of interactions, and considered homogeneous mosquito populations without distinguishing the four distinct stage of development during a lifetime: egg, pupa, larva, and adult. Continuously in this direction and considering the density dependence of the aquatic stages (egg, pupa, larva) but ignoring the density dependence of the adult, new revised simple models were formulated in [23], which are more reasonable for such studies.
While the larvae stages are implicitly included, the models in [23] and those models in [6] are based on ordinary differential equations (ODE) without considering the growth stage of the larvae mosquitoes. Following the line in [23], we explicitly include the time duration of the aquatic stages and formulate a delayed differential equation (DDE) model for the interactive wild and sterile mosquitoes but only consider constant releases of sterile mosquitoes in this paper. Detailed descriptions for the model formulation are given in Section 2. To fully analyze the model dynamics of the models, we begin with the case of absence of sterile mosquitoes in Section 3. We investigate the existence and the local and global stability for all equilibria of the model system. We then include the sterile mosquitoes and formulate a DDE model system to study the interactive dynamics of the wild and sterile mosquitoes in Section 4. With the assumption of constant releases of sterile mosquitoes, the dynamics of the sterile mosquitoes can be determined and then we investigate the coupled DDE model system. We derive thresholds of releases based on which the wild mosquito population will be eliminated or suppressed, or the two types of mosquitoes coexist. To testify our main theoretical results and exhibit the complexity of the model dynamics, we provide numerical examples in Section 5. Finally, we briefly discuss our main findings in this study in Section 6. 2. The model formulation. In this section, we give the motivation to formulate a delayed sterile mosquito model with a constant release rate.
Suppose that the wild mosquito population follows the logistic growth in the absence of interactions, and consider homogeneous mosquito population without distinguishing the metamorphic stages. Let S ν (t) be the number of wild mosquitoes and S g (t) the number of sterile mosquitoes at time t. The following continuous-time mathematical model for the interactive dynamics of the wild and sterile mosquitoes was formulated in [6]: where a is the number of offsprings produced per individual, per unit of time, with all matings; S ν (t)/ (S ν (t) + S g (t)) is the fraction of wild mosquitoes over all mosquitoes; b is the constant release rate of sterile mosquitoes; µ i and ξ i , i = 1, 2, are the density independent and dependent death rates of the wild and sterile mosquitoes, respectively. Notice that the intraspecific competition of mosquitoes basically takes place in water, and is relatively weak between adult mosquitoes. Then the death rates for adult wild and sterile mosquitoes are considered as density independent. Assume that the death rates of wild and sterile mosquitoes are the same. The revised models were formulated in [23], where the model dynamics are governed by the following system: Here ξ ν is the carrying capacity parameter such that 1 − ξ ν S ν describes the effect of density dependent survival probability of the larvae. Models (1)-(2) are both ODE models, and the progression of the larvae stages of mosquitoes is not explicitly presented. Let the duration of the larvae progression or the average maturation time be τ . Based on model (2), we formulate the following delayed model for the interactive wild and sterile mosquitoes where e −µ0τ is the survival rate of the immature mosquitoes who were born at time t − τ and are still alive at the time t. Suppose the initial value condition is given by Then for a given t , 3. Delayed model without sterile mosquitoes . We assume that, before sterile mosquitoes are released, the dynamics of the wild mosquitoes are governed by the following equation and that the population size in the beginning is less than the carrying capacity 1/ξ ν so that the initial condition for (6) is given by For given t 0 ≥ 0, S ν (t) = S ν (t, t 0 , φ ν ) is said to be a solution of (6) through Remark 1. Models similar to (6) have been used to describe the dynamics of stagestructured population or the population with adult recruitment in several previous studies [20,21,3,29], and model dynamics have been particularly investigated in [28,31] with more details. We briefly apply and extend existing results to model (6) as follows.
Define the intrinsic growth rate of the population by r 0 := ae −µ0τ /µ 1 . We have the following basic properties of the solutions of equation (6).
Proof. The local existence and uniqueness of the solution of (6)-(7) follows from the standard results of the theory of delay differential equations (See, e.g. [33]). The global existence of the solution in [t 0 , ∞) follows from the boundedness we prove below.
We now study the stability of the equilibria of equation (6) with the initial condition (7). We can easily obtain that (6) has a trivial equilibrium S  ν is globally asymptotically stable, and thus the population goes extinct, if 0 < r 0 < 1 and is unstable if 1 < r 0 < 4.
Proof. (1) The characteristic equation for the linearization of (6) at S It is easy to see that h(λ) = 0 cannot have any root with zero or positive real part and thus S (0) ν is locally asymptotically stable if r 0 < 1. We then prove the global attractiveness of S We then claim that every solution of (6) eventually enters the interval (0, C 0 ). In fact, if it does not occur, there exist sufficiently large t 4 > t 0 and a solution y(t) of (6) such that y(t) ≥ C 0 for all t > t 4 , but from (6) we have y (t) < −µ 1 C 0 +g 0 (C 0 ) < 0 for all t > t 4 + τ , which indicates that y(t) → −∞ as t → ∞. This leads to a contradiction to y(t) > C 0 for all t > t 4 .
We further investigate the stability of S ν through Hopf bifurcations for r 0 > 3. Denoteτ Then r 0 = 3 when τ = τ , and r 0 = 4 when τ =τ . The characteristic equation for the linearized equation of (6) around S (9) can be also written in the general form with is asymptotically stable when τ = 0. In the following, we consider the stability of S (1) ν when τ ∈ [0, τ ] by using the method introduced by Beretta and Kuang [4], where a geometrical criterion which gives the existence of purely imaginary roots of a characteristic equation with delaydependent coefficients was established. More specifically, we investigate the existence of purely imaginary roots λ = iω, (ω > 0), to equation (11), with the coefficient Q(λ, τ ) depending on τ .
In order to employ the criterion in [4], we first need to verify the following properties for all τ ∈ [0, τ ].
, as a function of τ , has a finite number of zeros; (v) Each positive root ω(τ ) of F (ω, τ ) = 0 is continuous and differentiable in τ whenever it exists.
In addition, by applying Corollary 2.4 in [32] and the Hopf bifurcation theorem of functional differentia equations [12], we can make the following conclusion.
We give an example below to demonstrate the bifurcation results in Theorem 3.3.
4. Delayed interactive model with sterile mosquitoes. Now suppose that sterile mosquitoes are released into the field of wild mosquitoes. Then the interactive dynamics are governed by system (3). For any equilibrium (S ν , S g ) of (3), we can easily find that it is globally asymptotically stable with respect to the second variable S g = b/µ 2 . Therefore, we only need to study its stability with respect to the first variable.
First, we study the property of the birth function of wild mosquitoes in (3).
(2) Suppose that S + ν < C * 1 . We divide the proof into two parts as follows.
Next, we claim that the first variable S ν (t) of (3) eventually enters the interval (0, C * 1 ). In fact, if it does not occur, there exists t 10 > t 0 large enough and a solution (y(t), S g (t)) of (3) such that y(t) ≥ C * 1 for all t > t 10 . However, from (3) we have y (t) ≤ −µ 1 C * 1 + g(C * 1 , b/µ 2 ) < 0 for all t > t 10 + τ , which indicates y(t) → −∞ as t → ∞ and leads to a contradiction to y(t) ≥ C * 1 for all t > t 10 . Finally, since every solution of (3) eventually enters the positively invariant set for all t > t 11 + τ . Hence the lower limit S ν ≥ S + ν . Let {s n } be a divergent sequence along which S ν (s n ) →S ν and S ν (s n ) → 0 as n → ∞. Taking the limit of the first equation of (3) along this sequence leads to Otherwise, there exists t 12 > t 0 such that S ν (t 12 ) = S + ν − ε and S ν (t) > S + ν − ε for all t ∈ [t 0 − τ, t 12 ). Hence S ν (t 12 ) ≤ 0. However, it follows from the first equation of (3) that (25), we have S + ν − ε < S ν (t) < S + ν . Hence the lower limit satisfies S ν ≥ S + ν − ε > S − ν , and the upper limit satisfies S ν ≤ S + ν . Let {w n } be a divergent sequence along which S ν (w n ) → S ν and S ν (w n ) → 0 as n → ∞. Taking the limit of the first equation of (3) along this sequence leads to (1) and (2), we immediately conclude that N − is unstable. We complete the proof.  (1): N 0 is stable, and lim t→∞ (S ν (t), Proof. (1). We skip the proof since it is virtually identical to that of Theorem 4.4 (1).

Model dynamics in case
In this case, f (S ν , b) > 0 for all S ν ≥ 0. Then N 0 is the only non-negative steady-state. We show that it is globally asymptotically stable in (3) below.
Remark 3. From Theorem 4.6, we determine a release threshold of sterile mosquitoes, above which the wild mosquito population will eventually be suppressed or eliminated. From Theorem 4.4, Theorem 4.5, and Theorem 4.6 above, we further conclude that b > b * is a sufficient and necessary condition for the global asymptotic stability of the zero solution of the delayed model (3).
5. Numerical simulations. We give numerical simulations, in this section, to demonstrate the results we obtain for model (3).

Estimation of parameters.
Before a simulation, it is critical to obtain a detailed estimation of the parameters for our model. We first provide important parameter values for Aedes albopictus and sterile mosquitoes, taken from earlier measurements in Guangzhou [2,8,25,26,27,37], in Table 1. Note that parameters related to environmental changes such as temperature and rainfall are not considered.
We let the survival probability of the average maturation time be 0.05 as in [26], and assume that the sterile mosquitoes die in about 7 days after releasing to the wild so that the death rate of sterile mosquitoes is 1/7. We then choose the model parameter values for our simulation as a = 6, τ = 30, e −µ0τ = 0.05, µ 1 = 0.08, µ 2 = 1/7, ξ ν = 0.0025.
The values of parameters a, τ , e −µ0τ , µ 1 , and µ 2 are within the range of our estimations given in Table 1.  (3) . We present simulations to illustrate the dynamics of system (3) as follows in Example 5.1. In FIGURE 3, it shows a bi-stable phenomenon when the release rate of sterile mosquitoes is smaller than the threshold of release rate b * , which is consistent with the results of Theorem 4.4. FIGURE 4 shows a clear description of the global dynamics of (3) when the release rate of sterile mosquitoes is larger than the threshold of release rate b * , which results in the global asymptotic stability of the zero solution of system (3).  : Simulations for the case of b > b * . Parameters are given in (28) such that the release threshold b * ≈ 28.81 and we set b = 35 greater than b * . Then the zero equilibrium N0 is globally asymptotically stable, which indicates that the wild mosquito population is eventually eradicated.

5.3.
The effects of the releases on suppression efficiency. In this section, we investigate the effects of different release rates of sterile mosquitoes on suppression efficiency when the release rate b is greater than the threshold b * . The suppression efficiency is measured by the least time needed to have 95% of the initial wild mosquitoes suppressed, which is usually used in a field trial or a laboratory experiment to suppress a wild mosquito population.
Example 5.2. Let the parameters of system (3) be given in (28), and φ ν = 100, φ g = 200. Then the threshold of release rate is again b * ≈ 28.81. We take b = 50, 70, 90, and 130, respectively. While our simulation shows that the wild mosquito population will be eventually suppressed in 300 days, comparisons of the four cases indicate that the efficiency is improved as b increased such that if we want to suppress the wild mosquito population successfully, that is to reduce the wild mosquito population to 95% of its initial size, in 200 days, we should take a release rate b larger than 70 as shown in FIGURE 5.  6. Concluding remarks. The use of sterile mosquitoes has provided to be an effective and safe method to biologically control mosquito populations for prevention of mosquito-borne diseases, such as dengue, Zika, yellow fever, and west Nile. There are many works in the literature on mathematical modeling of releases of sterile mosquitoes and their interactive dynamics with wild mosquitoes. Comparing with the existing models (1)-(2) in [6,23], we take the maturation delay of wild mosquitoes into consideration and formulated our dynamical models with time delay in this paper. We first gave complete descriptions in Section 3 for the model without sterile mosquitoes, where the model dynamics are characterized by the intrinsic growth rate of the population r 0 . A numerical example to confirm the results was also provided in Section 3.
We then included the sterile mosquitoes with a constant release rate in the model in Section 4 and fully investigated the interactive dynamics for model system (3). We determined a release threshold of sterile mosquitoes b * such that the model system has no, unique, or two positive equilibria when the release rate is greater than, equal to, or less than the threshold. For the three possible cases, we further studied the stability of each of the equilibria. In addition, we established a sufficient and necessary condition for the global asymptotic stability of the zero solution.
Time delay is included in the model formulation in this study, which is needed and more realistic from the biological perspective. The mathematical analysis becomes more difficult and the thorough investigations of the model features and detailed techniques developed in this paper help us to gain deeper understanding of those features. In the meantime, nevertheless, we realize that the model dynamics such as the bi-stability revealed in this study are similar to those of the homogeneous models and stage-structured models in [6,23]. This seems to imply that the dynamical features for constant releases of sterile mosquitoes are possibly model independent and characterize the biological essentials of such a biological control measure. Moreover, this might be also another evidence that simplified models may not necessarily lose key features that more complicated models exhibit, as pointed out in [23].