Dynamics of interactive wild and sterile mosquitoes with time delay

ABSTRACT We develop a delay differential equation model for the interactive wild and sterile mosquitoes. Different from the existing modelling studies, we assume that only those sexually active sterile mosquitoes play a role for the interactive dynamics. We consider the cases where the release amount is either constant or described by a given function of time. For the constant releases, we establish a threshold of releases to determine whether the wild mosquito suppression succeeds or fails. We study the existence and stability of the model equilibria. When the releases are described by given functions, the trivial equilibrium is no longer globally but locally uniformly asymptotically stable if the amount of releases is below the threshold whereas it is still globally uniformly asymptotically stable if the release amount is above the threshold. Numerical examples demonstrating the model dynamical features and brief discussions of our findings are also provided.


Introduction
In the face of outbreaks of mosquito-borne diseases, such as dengue fever, malaria and Zika, plaguing the world, traditional prevention and control programs are to eradicate or suppress wild mosquitoes, involving the use of insecticides. The spraying insecticides bas been an effective method to control these disease outbreaks in the short term in some areas. However, the method of spraying insecticide does not typically last long enough to keep the mosquito population density low enough to avoid the risk of epidemic. Regular spraying of insecticides can also lead to serious environmental problems and insecticide resistance [22]. These problems lead to the urgent need of finding environmentally-safe ways to control the outbreaks of the mosquito-borne diseases. Through the continuous efforts of scientists, effective biological control methods have been developed and used. These include the genetic approaches [9,15,21,24], the Wolbachia driven mosquito control technique [4,25,27], and the sterile insect technique (SIT) [26].
Mathematical models have been developed and analysed to investigate the effects of these biological control measures. Modelling of transgenesis and paratransgenesis has helped us to gain insight into better strategies in genetically engineering mosquitoes [16,17,19,20]. Models based on difference or differential equations to study the suppression effects on the control of mosquitoes by releasing Wolbachia-infected male mosquitoes have been formulated with the earliest difference equation model, incorporating laboratory data, by Caspari and Watson [7] in 1959. The work then has been followed one after another in [11][12][13][14][29][30][31][32][33] and the references therein. To investigate the impacts of utilizing SIT on mosquito control, mathematical models have also been formulated and analysed in [1][2][3]6,[16][17][18]. Several strategies of releases of sterile mosquitoes are proposed, qualitatively analysed, and compared in these studies. Different from the methods of transgenesis and paratransgenesis, releasing sterile or Wolbachia-infected male mosquitoes has a common basis in mathematical modelling process.
For these biological control measures, the major difficulty is to estimate how many among the released engineered mosquitoes are sexually active. Existing models involving sterile mosquitoes assume the entire population of sterile mosquitoes is sexually active and its population dynamics are governed by a different differential equation coupled with an equation for the dynamics of the wild mosquitoes. However, many sterile mosquitoes have a short sexual lifespan and lose their ability to mate. As a result, their presence has no influence on the dynamics of wild mosquitoes after they lose their mating ability.
We, in this paper, focus on further developing models for the interactive wild and sterile mosquitoes with the idea that only sexually active mosquitoes are considered and their death is ignored. We describe our model formulation in detail in Section 2. We consider two different cases of releases of sterile mosquitoes. Model dynamics for the sterile mosquitoes staying at a constant level and varying as a given function are investigated in Sections 3 and 4, respectively. Numerical examples are also provided to confirm our analytic results. Brief discussions are given in Section 5.

Model formulation
To explore the interactive dynamics for mosquitoes after sterile mosquitoes are released into the field of wild mosquitoes, the following ordinary differential equations model was introduced in [6]: (1) Here w and g are the numbers of wild and sterile mosquitoes at time t, respectively, a is the total number of offspring per wild mosquito, per unit of time, w/(w + g) is the fraction of mates with wild mosquitoes, μ i and ξ i (w + g), i = 1, 2, are the density-independent and density-dependent death rates of wild and sterile mosquitoes, respectively, and B is the rate of releases of sterile mosquitoes.
Since mosquitoes undergo complete metamorphosis, going through four distinct stages of development during a lifetime [8], inclusion of the larvae maturation to adults makes more realistic models. Instead of explicitly formulating stage-structured models, by using time delay τ for the larvae maturation of wild mosquitoes, the following delay differential equations model was studied in [5]: where e −μ 0 τ is the survival rate of lava mosquitoes, and the initial conditions for the model system are with t 0 ≥ 0, and φ(θ) and ψ(θ) both positive and continuous in [−τ , 0]. We note that in both of models (1) and (2), the sterile mosquitoes are assumed to have their own dynamic equation. However, among the sterile mosquitoes determined by the second equation of both (1) and (2), some are sexually active, which can play a role in the interactive dynamics, but some have lost their abilities to mate. The interactive dynamics of wild and sterile mosquitoes after sterile mosquitoes are released are mostly affected only by the mating between the two types of mosquitoes. Thus, the modelling focus should be only on those sexually active sterile mosquitoes. Recently the author in [28] introduced a non-autonomous mosquito population suppression model, in which it is assumed that the number of Wolbachia-infected mosquitoes released is a given positive function. It is also assumed that the released mosquitoes have no offspring and their losses can only be supplemented by new releases.
Employing the modelling ideas in [28] and assuming that the death and the dynamics of the sterile mosquitoes with respect to the time variable are negligible, we use the following model equation for the wild mosquitoes where g(t) is the number of sterile mosquitoes released. It is a non-negative valued function of t ≥ 0 that satisfies g(t) ≡ 0 for t < 0. Comparing (3) with both (1) and (2), we see that g(t) is no longer a variable whose dynamics are governed by an independent equation, and Equation (3) is no longer an autonomous, but a non-autonomous or a time-varying equation. Many classical methods and techniques, such as qualitative analysis and the characteristic root method, are then not applicable.
For given t 0 ≥ 0, with an initial function φ ∈ C([t 0 − τ , t 0 ], (0, ∞)), and a given nonnegative function g(t) specified separately on . It is easy to see that, since φ(t) > 0 for t 0 − τ ≤ t ≤ t 0 , we have w(t) > 0 for all t ≥ t 0 . We assume that the release starts at time t = 0. Then g(t) = 0 for t ≤ 0 and g(0 + ) denotes the number of the released sterile mosquitoes for the first time.
Before we proceed with further analysis of model dynamics, a few definitions are given below. Definition 2.1: [28] A solution w 0 (t) of (3) is said to be • Stable if for any t 0 ≥ 0 and ε > 0, there exists δ > 0 such that and If δ is independent of t 0 , then w 0 (t) is uniformly stable. • Stable from the right-side (or left-side) if for any t 0 ≥ 0 and ε > 0, there exists • Semi-stable if it is stable from the right-side (or left-side) but unstable from the left-side (or right-side). • Globally asymptotically stable if it is uniformly stable, and for any t

Sterile mosquitoes staying at a constant level
In reality, sterile mosquitoes cannot be released continuously, but only at discrete times. We now assume that the sterile mosquitoes are released impulsively with a constant amount c at time kT, k = 0, 1, . . ., where T > 0. The lifespan of male mosquitoes is relatively short, usually less than 10 days. While mating behaviour, one of the critical behaviours that characterize the mosquito life strategy, is the least understood and most understudied [10,23], it is well known that males can mate several times and their sexual lifespan is even much shorter than their ages. During the period when the sterile mosquitoes are sexually active, it is reasonable to neglect their death, and more importantly the efficacy of the releases of sterile mosquitoes is closely related to how often they are released and how the waiting time for the next release is correlated with the sexually active period of the sterile mosquitoes.
Let the sexual lifespan of sterile male mosquitoes beT. If the new sterile mosquitoes are released exactly at the end of the sexual lifespan of sterile mosquitoes such that T =T, then the sterile mosquitoes in the field are almost kept at a constant level. In this case, we let g(t) = c for all t ≥ 0. Then model equation (3) becomes An equilibrium of (4) satisfies the following quadratic equation in w: Define the threshold value of releases for (4) by We have the following existence results for the equilibria of Equation (4). (4) has a trivial equilibrium E 0 := 0 and, in addition,
For the case when c ∈ (0, g * ), Equation (4) exhibits more complex dynamics such as the bistable dynamics. We give relatively complete stability results below.
On the other hand, from (4), we have Because {w(t n − τ )} is bounded, it has a convergent subsequence. For convenience, we assume lim n→∞ w(t n − τ ) = w 1 . Then w 1 ≤w. Taking the limit on both sides of (10) yields and then . This contradicts (8), and as a result, we get w = 0 and lim t→∞ w(t) = 0. Hence (9) is proved.
In fact, if (11) is not true, then there iss > t 0 as the least time such that w(s) = E + (c) + ε or w(s) = E + (c) − ε. If w(s) = E + (c) + ε, then w (s) ≥ 0. It follows from (4) that This is a contradiction to the fact that Q(w, c) < 0 for x ∈ (E − (c), E + (c)). This proves the uniform stability of E + (c).
Finally, we show By using an argument similar to that in the proof of (11), we can show w(t) > E − (c) for t ≥ t 0 and hence the lower limit w = lim inf t→∞ w(t) ≥ E − (c). Let {s n } be a monotonic time sequence such that s n → ∞, w(s n ) → w, as n → ∞, and w (s n ) ≤ 0. Then we have, from (4), that Since {w(s n − τ )} is bounded, we let lim n→∞ w(s n − τ ) = w 2 . Then w 2 ≥ w. Taking the limit on both sides of (12) yields Thus w ≥ E + (c), and furtherw ≥ E + (c). The proof will then be completed if we can provē w ≤ E + (c).
To this end, we let {t n } be an infinite divergent sequence along which w(t n ) →w, as n → ∞, and w (t n ) ≥ 0. Then from (4), we obtain (10). By using an argument similar to that in the proof above in (1), we have , and hence the proof is complete.
If the number of sterile mosquitoes released is the same as the threshold value such that g(t) ≡ g * , for all t ≥ 0, then Equation (3) is reduced to (4) with c replaced by g * where g * is given in (6); that is, Equation (3) becomes with a unique positive equilibrium The dynamics of (13) can be summarized as follows. (1) Equilibrium E 0 = 0 is uniformly asymptotically stable and lim t→∞ w(t, g * , t 0 , φ) = 0, if 0 < φ(t) < E * for t ∈ [t 0 − τ , t 0 ]; (2) Equilibrium E * is globally uniformly asymptotically stable from the right-side that

Proof: (1)
The proof is the same as that in Theorem 3.2 (1) and is omitted.
(2) To show the uniform stability of E * from the right-side, it suffices to verify that 0 < φ(t) − E * < ε on [t 0 − τ , t 0 ] for any ε > 0 implies By an argument similar to that in the proof of Theorem 3.2, we can easily prove that (14) is not true, then there is t > t 0 such that w(t ) = E * + ε, and w(t) < E * + ε for t ∈ [t 0 − τ , t ) and w (t ) ≥ 0. From (13), we get which yields Q(E * + ε, g * ) < 0, a contradiction to the non-negativity of Q(x, g * ). This shows that (14) is true and hence E * is uniformly stable from the right-side. We then prove that lim t→∞ w(t) = lim t→∞ w(t, g . To this end, we just need to provew = lim sup t→∞ w(t) = E * . Otherwise, if w > E * , then there must be an increasing sequence {t n } such that t n → ∞, w(t n ) →w, as n → ∞, and w (t n ) ≥ 0.
From (13), we have Since {w(t n − τ )} is bounded, it has a convergent subsequence. For simplicity, we let lim n→∞ w(t n − τ ) = w 1 , and then w 1 ≤w. Taking the limit on both sides of (15), we have which leads to Q(w, g * ) ≤ 0, a contradiction to the fact that Q(w, g * ) > 0 for all nonnegative x = E * . This forces Q(w, g * ) = 0 and hencew = E * . The proof is complete.
Remark 3.1: Theorem 3.3 tells us that E * is semi-stable, that is, stable from the rightside and unstable from the left-side. Such a phenomenon of semi-stability is due to the coincidence of unstable equilibrium E − (c) and stable equilibrium E + (c) when c = g * .
We give an example below to demonstrate the results from Theorems 3.2 and 3.3.
we have threshold g * = 7.4304. With the release c = 6.6874 < g * , there exist two positive equilibria E − = 4.1171 and E + = 14.1108. The origin E 0 is uniformly asymptotically stable, E − is unstable, and E + is uniformly asymptotically stable. Solutions with initial values Figure 1. With parameters given in (16) and t 0 = 0, the threshold of releases is g * = 7.4304. When the amount of releases is c = 6.6874 less than the threshold, there exist two positive equilibria E ∓ . The origin E 0 is uniformly asymptotically stable, E − is unstable, and E + is uniformly asymptotically stable. Solutions approach either E 0 or E + , depending on their initial values, as shown in the left figure. When the amount of releases is equal to the threshold g * , there exists a unique positive equilibrium E * = 8.3709 which is globally uniformly asymptotically stable from the right-side such that lim t→∞ w(t, g * , 0, φ) = 0, for φ(t) < E * , t ∈ [−9, 0], and lim t→∞ w(t, g * , 0, φ) = E * , for φ(t) > E * , t ∈ [−9, 0], as shown in the right figure. φ(t) < 4.1171, for t ∈ [−9, 0], approach E 0 and solutions with initial values φ > 4.1171, for t ∈ [−9, 0], approach the positive equilibrium E + , as shown in the left figure in Figure 1.
In the case of c > g * , equilibrium E 0 = 0 is the unique equilibrium. Its dynamics are described as follows.

Theorem 3.4:
Assume that c > g * . Then the origin E 0 of (13) is globally uniformly asymptotically stable.
The proof follows directly from that of Theorem 4.2 and is omitted. From Theorems 3.2, 3.3, and 3.4, it is easy to obtain the following useful conclusion.

Corollary 3.5:
The origin E 0 of (13) is globally uniformly asymptotically stable if and only if c > g * . Remark 3.2: Corollary 3.5 provides a necessary and sufficient condition for the global asymptotic stability of E 0 . In comparison, the conditions for the global asymptotic stability of E 0 given in [5] (for the case when the dynamics of g(t) is governed by a differential equation) are only sufficient but not necessary.

Sterile mosquitoes varying as a function of time
In this section we no longer assume T =T and assume that g(t) is not a constant but a nonnegative given function of time for all t > 0. Then (3) has no positive constant equilibrium, and we have the following results. Theorem 4.1: Assume that there exist two positive constants g 1 < g 2 ≤ g * such that g(t) ∈ [g 1 , g 2 ] for all t > 0. Then the trivial equilibrium E 0 of (3) is uniformly asymptotically stable. Moreover, lim t→∞ w(t, g, t 0 , φ) . Here E − (g i ), i = 1, 2, are given in (7) with c replaced by g i .

Proof:
We first show that E 0 is locally uniformly stable.
We next show that It suffices to prove that Assume, otherwise, w > 0. It is easy to prove w < E − (g 1 ). Then we may choose a strictly monotone increasing sequence {t n } such that w (t n ) ≥ 0 and w(t n ) → w as n → ∞. It follows from Equation (3) that which yields, by substituting g(t) ≥ g 1 , for all t > 0, into (22), Since w(t n − τ ) is bounded, from (18), sequence {w(t n − τ )} contains a convergent subsequence, such that lim n→∞ w(t n − τ ) = w 1 ≤ w for the elements in the subsequence. Then, taking the limit on both sides of (23) as n → ∞, we have which leads to This implies w ∈ [E − (g 1 ), E + (g 1 )], a contradiction to w < E − (g 1 ). Thus (21) holds and therefore E 0 is locally uniformly asymptotically stable. Finally, we show that To this end, it suffices to show that if We prove this by contradiction again as follows. Otherwise, we can chooset > t 0 such that w(t) = E − (g 2 ) and w(t) > E − (g 2 ) for t 0 − τ < t <t. Hence w (t) ≤ 0 and it follows from (3) that which leads to a contradiction. The proof thus is complete.
For the case of g(t) ≡ c > g * in Section 3, Equation (3) becomes (4) and has only the unique trivial equilibrium E 0 = 0 and it is globally uniformly asymptotically stable. Now we show that the dynamics are similar to those for the case of variable releases as follows.

Proof:
The uniform stability of E * can be proved by an argument similar to that in the proof of Theorem 3.2. To complete the proof, it suffices to show that for any In fact, if we setw = lim sup t→∞ w(t), then we only need to provew = 0. If the limit lim t→∞ w(t) exists, it is easy to show that (24) holds. Now, we assume that the limit lim t→∞ w(t) does not exist. Then there is an increasing sequence {s n } such that w(s n ) →w as n → ∞ and w (s n ) = 0. From (3), we have Given that {w(s n − τ )}, {g(s n )}, and {g(s n − τ )} are bounded, they all have convergent subsequences. For convenience and without loss of generality, we let lim n→∞ w(s n − τ ) = w 1 , lim n→∞ w(s n ) = g 1 , and lim n→∞ g(s n − τ ) = g 2 . Then g 1 , g 2 > g * and by taking limit on both sides of (25), we have which leads to Q(w, g * ) < 0, a contradiction to the non-negativity of Q(w, g * ) for all nonnegative w. The proof is complete.
We provide an example below to confirm the results from Theorems 4.1 and 4.2. The parameters are the same as in Example 3.1 but the release functions are no longer constant functions.

Concluding remarks
We formulated a new model (4) for the interactive dynamics between wild and sterile mosquito populations. This model includes a time delay due to the larval stage of the mosquito life cycle. Adopting modelling methods in [28], which are based on the short sexual lifespan of male mosquitoes, we ignored the dynamics of the sterile mosquito population. In place of a dynamical equation for sterile mosquitoes, we introduced a function g(t) which is the amount of sterile mosquitoes released into the wild population. We then first assumed that the releases of sterile mosquitoes are kept as a constant so that g(t) = c, and gave complete mathematical analysis for the model dynamics in Section 3. We derived a threshold value g * for the releases and showed that if the release amount exceeds the threshold with c > g * , the origin E 0 is globally uniformly asymptotically stable; that is, all wild mosquitoes are wiped out regardless of their initial sizes. If the release amount is below the threshold with 0 < c < g * , the origin E 0 is locally uniformly asymptotically stable and there exist two positive equilibria E ∓ such that E − is unstable and E + is uniformly asymptotically stable. In this case, wild mosquito populations of high density will not be eliminated, whereas those of low density will be. Although it may rarely happen biologically, for mathematical completion, we also showed that at the critical situation where the release amount is exactly as g * , for all t ≥ 0, E 0 = 0 is uniformly asymptotically stable, and there exists a unique positive equilibrium which is globally uniformly asymptotically stable from the right-side.
When the releases of the sterile mosquitoes are given as non-constant functions of time, the model equation becomes a non-autonomous equation with time delay, and its mathematical analysis is much more challenging. We showed, in Section 4, that if inf t>0 g(t) exceeds the threshold g * , the origin E 0 is globally uniformly asymptotically stable, whereas if g(t) ∈ (0, g * ) and inf t>0 g(t) > 0, E 0 is only locally but not globally uniformly asymptotically stable. Note that the release function g(t) can be any given function of t. The model equation opens a door for modelling other ecological situations such as periodic and impulsive release functions. Further studies and new investigations are to appear in the near future.