Periodic Orbits of a Mosquito Suppression Model Based on Sterile Mosquitoes

: In this work, we investigate the existence and stability of periodic orbits of a mosquito population suppression model based on sterile mosquitoes. The model switches between two sub-equations as the actual number of sterile mosquitoes in the wild is assumed to take two constant values alternately. Employing the Poincaré map method, we show that the model has at most two T -periodic solutions when the release amount is not sufﬁcient to eradicate the wild mosquitoes, and then obtain some sufﬁcient conditions for the model to admit a unique or exactly two T -periodic solutions. In particular, we observe that the model displays bistability when it admits exactly two T -periodic solutions: the origin and the larger periodic solution are asymptotically stable, and the smaller periodic solution is unstable. Finally, we give two numerical examples to support our lemmas and theorems.


Introduction
Mosquito-borne diseases, such as dengue, malaria, Zika, yellow fever, West Nile fever, seriously imperil the health of people in tropical and sub-tropical areas. Recently, the rising global temperatures and uncontrolled urbanization greatly expand the threat range of such diseases [1][2][3]. The most direct approach to prevent and control these mosquito-borne diseases is to control the vector mosquitoes. Sterile Insect Technique (SIT) and Incompatible Insect Technique (IIT) are two promising and environment-friendly methods for controlling mosquitoes [4][5][6][7][8][9], which involve rearing massive number of sterile or Wolbachia-infected mosquitoes (these two kinds of mosquitoes are called sterile mosquitoes hereafter) in the laboratories or factories, then releasing them into the field to sterilize wild mosquitoes, and thus, the mosquito density can be asymptotically controlled.
Mathematical modeling has played an important role in the control and prevention of mosquito-borne diseases. One is to study the Wolbachia spread dynamics in mosquito populations aiming to study the possible invasion of Wolbachia in mosquito populations [10][11][12][13][14][15][16][17][18][19], and the other is to study the interactive dynamics of wild and sterile mosquitoes aiming to eliminate or control the wild mosquitoes below the risk of mosquito-borne disease outbreaks. Most early studies on the interactive dynamics mainly concentrate on models where the number of sterile mosquitoes was assumed to be an independent variable satisfying an independent dynamical equation, see, for example, [20][21][22][23][24], whereas recent investigations have proposed a motivation-attention modeling idea where the number of sterile mosquitoes is treated as a control function instead of an independent variable, since the fundamental and only role of sterile mosquitoes is to mate with wild mosquitoes to induce the infertility of wild females, see, for example, [25][26][27][28][29][30][31][32][33][34][35][36][37][38].
In [21], Li formulated the following model to explore the interactive dynamics between wild and sterile mosquitoes. In this model, w = w(t), g = g(t) are the numbers of wild and sterile mosquitoes at time t, respectively, a is the total number of offspring produced per wild mosquito, per unit of time, w w+g is the fraction of mates with wild mosquitoes, ξ represents the carrying capacity of wild mosquitoes such that 1 − ξw characterizes the density-dependent survival probability, µ is the density-independent death rate coefficient of wild as well as sterile mosquitoes, and B(·) is the rate of the releases of sterile mosquitoes. Under three different release strategies: constant releases (B(w) = b), proportional releases (B(w) = bw), and fractional releases (B(w) = bw 1+w ), the author obtained the corresponding conditions for the existence and stabilities of positive equilibria.
Note that the sterile male mosquitoes are released to sterilize wild mosquitoes, along with the two facts: the lifespan of mosquitoes is short, and the sexual lifespan of those released sterile mosquitoes is much shorter than their longevity [39,40], it seems reasonable to ignore the death of those released sterile mosquitoes when they are sexually active and energetic. Moreover, the number of sterile mosquitoes that have been released into the field at time t is known in advance [9]. Hence, by regarding the number g(t) of released sterile mosquitoes at time t as an arbitrarily given nonnegative continuous function, the authors in [28] omitted the second equation in (1), and studied the global dynamics of wild mosquitoes by the following model This treatment greatly simplifies the model and makes the analysis mathematically more tractable.
In designing the release strategies of sterile mosquitoes, there are three important parameters: c, the release number in each batch; T, the waiting period between two adjacent releases;T, the sexual lifespan of sterile mosquitoes. From the perspective of the relations between T andT, there are three scenarios: T >T, T =T and T <T. In [38], we studied the first case, and split Equation (2) into the following two sub-equations and where i = 0, 1, . . . We defined the first release amount threshold Under the assumption of c > g * , we further obtained the second release amount threshold c * and the waiting period threshold T * as follows We investigated Equations (3) and (4) for exploring the dynamical behaviors of wild mosquitoes under the interferences of released sterile mosquitoes.
However, rearing sterile mosquitoes needs huge economic input and puts the density of wild mosquitoes below a given value which is more realistic than eradicating them. Taking these facts into consideration, in this paper, we mainly discuss the dynamics of the model (3) and (4) under the release strategy of c ≤ g * , which concentrates on reducing the number of wild mosquitoes with a lower cost than that in [38]. The remaining of the paper is structured as follows. In Section 2, we present some preparations, which mainly contain the introduction of the Poincaré map for determining the number of periodic solutions of the model (3) and (4), and some lemmas that are of vital importance to the proof of our main results. We give our main results in Sections 3 and 4. In Section 5, we provide two numerical examples to verify our lemmas and theorems. Finally, some discussions on the current as well as future works are offered in Section 6.

Preliminaries
In this section, we give some lemmas, which provide a springboard to the proof of our main result. Lemma 1 ([38]). Assume that c > g * and T >T. Then, the following conclusions are true.
(1) If c ≥ c * and T ≤ T * . Then, the origin, denoted by E 0 , is globally asymptotically stable; (2) If T > T * or g * < c < c * and T = T * . Then, the model (3) and (4) has a unique T-periodic solution, which is globally asymptotically stable; (3) If g * < c < c * . Then, E 0 is locally asymptotically stable if and only if T < T * .
We next give the definition of a solution of the model (3) and (4). A function w(t) is said to be a solution of the model (3) and (4), if it satisfies both (3) and (4) Then, w(t) is a continuous and piecewise differentiable function defined on [0, +∞). Denote the solution of (3) and (4) with an initial value w(0) = u > 0 by w(t) = w(t; 0, u). Then, w(t) is a T-periodic solution of the model (3) and (4) provided that w(T) = w(T; 0, u) = u. For convenience, we set h(u) := w(T; 0, u),h(u) := w(T; 0, u) throughout this paper. Obviously, both h(u) andh(u) are continuously differentiable functions in u > 0.
As in [33], we define the following two sequences h n (u) andh n (u) by   Similar to Lemma 2.9 in [30], we give the following necessary and sufficient condition for E 0 to be asymptotically stable. The proof is also similar to that in [30] and we omit it here.

Lemma 3.
The origin E 0 is asymptotically stable if and only if there is δ 0 > 0 such that h(u) < u, for u ∈ (0, δ 0 ).
Lemma 2 tells us that the sign of h (u) − 1 plays a key role in determining the number of T-periodic solutions of the model (3) and (4). To obtain this sign, we first solve Equation (3) for t ∈ [0,T) to get the expression forh(u), and solve Equation (4) for t ∈ [T, T) to get the relation betweenh(u) and h(u). Then, we take the derivative of both sides of the above two expressions to get the expressions forh (u) and h (u). We first consider Equation (4) as it is independent of c.
When t ∈ [T, T), Equation (4) is where A = a−µ aξ . Further computations offer Integrating (7) fromT to T, we obtain where m = e −aAξ(T−T) . When t ∈ [0,T), we solve Equation (3) to obtain the expression forh(u). Since a solution may depends on the relations between c and g * , we consider the following three possible cases. Case (1): 0 < c < g * . In this case, we can rewrite Equation (3) in the form As 0 < c < g * , we know that (w − A 2 ) 2 − µ aξ (g * − c) = 0 has two real roots: Then, we write (9) as By simple computations, we have where Integrating (11) from 0 toT, along with the facts Equation (13) is equivalent to where We then transform (14) to which gives, by integrating from 0 toT, where Owing to the above preparations, we now give the following lemma, which manifests that the model (3) and (4) has no T-periodic solutions when u ≥ A.
Proof. Due to the continuous dependence of the solutions on the initial data, the solution w(t; 0, u) is strictly increasing with respect to u, and it is strictly decreasing with respect to Figure 1A. Thus, when u > u A , we have h(u) = w(T; 0, u) < u, which shows that the lemma holds in this case. When u ∈ [A, u A ), we have thus, the lemma also holds in this case. This completes the proof.  (3) and (4) such that the solution w(t) is strictly decreasing for t ∈ [0, T] when u > u A ; (B): When 0 < c < g * , the movement trend of the solution of the model (3) and (4) The following lemma tells us that the model (3) and (4) (6) and (13) that h(E * ) > E * , which proves the lemma in this case.
Proof. We divide the proof into two cases of 0 < c < g * and c = g * . Define H(u) = h(u) − u. For the case when 0 < c < g * , assume by contradiction that H(u) = 0 has at least two positive real roots for u ∈ (0, E 1 (c)), denote the largest two adjacent roots by u 1 and u 2 , Figure 2 for illustration. For convenience, we let E i (c) = E i with i = 1, 2 hereafter. From (12), we havē Calculating the derivative on both sides of (19), we obtain which, combining with (19), gives From (8), we get which indicates Moreover, by taking the derivative on both sides of (21), we arrive at which yields, along with (20), Substituting (22) into (23), we obtain From the two facts that h(u 2 ) = u 2 and h (u 2 ) ≥ 1, we get which, along with E 1 + E 2 = A, gives where the three coefficients B, C, D are defined as Define Then, (26) becomes χ(u 2 ) ≥ 0. For the relation between h (u 1 ) and 1, we need to consider two cases: (I) h (u 1 ) ≤ 1 (see panels (A), (B) and (C) in Figure 2 for illustration), and (II) h (u 1 ) > 1 (see panel (D) in Figure 2 for illustration). Next, we prove that both cases (I) and (II) are impossible.
For the case when c = g * , using similar arguments that led to (25), we have the remaining proof is similar to that of case 0 < c < g * and is omitted. We complete the proof.
The following lemma tells us that the model (3) and (4) has a unique T-periodic solution when u ∈ (E 2 , A).
For case (I), we have where χ(u) is defined in (28). To complete the proof, we consider two cases of χ(ū) > 0 and χ(ū) = 0 separately. Since we already know that B > 0 and χ(u) is a quadratic polynomial, the first case does not hold. Now, we focus on the latter one. There exist two sub-cases for χ(ū) = 0, that is, (a)ū ∈ (ū 1 ,ū 2 ); and (b)ū =ū 1 orū =ū 2 . By similar arguments used in the situation χ(ū) > 0, we can show that sub-case (a) does not hold. For sub-case (b), we first use the perturbation method to discuss the case of u =ū 1 .  To this end, let k > 1 be small enough such that h(u) − ku = 0 has exactly three positive real roots, denoted by u (i) See Figure 4B for illustration.
which yields, by (33), with B(k), C(k), D(k) being defined as follows: Hence, we have B(k) > 0, from the fact m < 1 < k. Set Then, we have, from (34), which contradicts B(k) > 0 and χ k (u) is a quadratic polynomial. This proves thatū =ū 1 is impossible. Moreover, ifū =ū 2 , then H (ū) = 0, and hence χ(ū) = 0. However, we find that which derives a contradiction, sinceū 1 <ū < A and B > 0 and χ(u) is a quadratic polynomial. Thus, sub-case (b) is also impossible, and hence case (I) is impossible. For case (II), we obtain H (ū i ) = 0, i = 1, 2, and H (ū) < 0, which leads to χ(A) > 0, a contradiction to χ(A) < 0, and hence case (II) is also impossible, which manifests that the lemma holds for the case when 0 < c < g * . For the case c = g * , by using similar arguments to those in the proof of the case 0 < c < g * , we can prove the lemma also holds in this case. The proof is thus completed.

At Most Two Periodic Solutions
In this section, we give the main theorem, whose proof involves Lemmas 1-7, and tackles three different cases under the release strategy of T >T. Proof. We divide the proof into the following three cases: For case (i), Lemmas 4-7 tell us that the model (3) and (4) has no T-periodic solutions when u ∈ [E 1 , E 2 ] ∪ [A, +∞), exactly one T-periodic solution when u ∈ (E 2 , A), at most one T-periodic solution when u ∈ (0, E 1 ); these facts manifest that the model (3) and (4) has at most two T-periodic solutions in this case.
For case (ii), from Lemma 1(2), we know that the model (3) and (4) has a unique T-periodic solution for the case when T ≥ T * . Hence, we only need to consider the case of T < T * .
Assume by contradiction that the model (3) and (4) has at least three T-periodic solutions. Let w 1 (t) = w(t; 0, u 0 1 ) and w 2 (t) = w(t; 0, u 0 2 ) be the two periodic solutions with the smallest and largest initial value, respectively. From Lemma 4, we obtain that u 0 i ∈ (0, A), i = 1, 2. From Lemma 1(3), we know that E 0 is locally asymptotically stable, which gives If the model (3) and (4) has exactly three T-periodic solutions, we denote the third Tperiodic solution by w 3 (t) = w(t; 0, u 0 3 ). For the relations between h (u 0 3 ) and 1, as shown in Figure 5, we need to consider two possible cases: h (u 0 3 ) < 1 and h (u 0 3 ) ≥ 1. Here, we only consider the case of h (u 0 3 ) ≥ 1, as the analysis for the case of h (u 0 3 ) < 1 is similar and is omitted. There exists u 0 Then, we have By (35), (16) becomes which gives, by taking the derivative on both sides, substituting (36) and (37) into (38), we obtain Taking the derivative on both sides of (8), we get Substituting (8) and (39) into (40), we reach Then, we have and Define Then, (41) and (42) imply Note that B 0 > 0 and Ω(u) is a quadratic polynomial and u 0 1 ≤ u 0 4 ≤ u 0 3 ≤ u 0 2 , which violate (43). The proof is completed in this case.
For case (iii), Lemma 1(1) and Lemma 1(2) manifest that the model (3) and (4) has at most one T-periodic solution. We complete the proof.

A Unique and Exact Two Periodic Solutions
Lemmas 4-7 imply that the model (3) and (4) has at least one and at most two Tperiodic solutions when 0 < c ≤ g * . In this section, we give the corresponding sufficient and necessary conditions for the model (3) and (4) to admit a unique or exactly two Tperiodic solutions, which shows that the wild mosquitoes can be eradicated only if their initial population size is less than a given value.
The proof of the global asymptotic stability ofŵ(t) is similar to that of Theorem 3.4 in [38] and is omitted.
Finally, we prove (49) is true. It follows from (47) that h(u) < u for u ∈ (û, +∞). By Lemma 2, we know that sequence {h n (u)} is strictly decreasing. Then, lim n→∞ h n (u) =û, which means that (49) is true if u >û. For the case when u ∈ ( u,û), the proof is similar to that in the case of u >û and is omitted. This completes the proof.
The following theorem provides a sufficient condition which ensures that the model (3) and (4) admits exactly two T-periodic solutions.
Theorem 3. Assume that c ∈ (0, g * ] and T < T * . Then, the model (3) and (4) admits exactly two T-periodic solutions. Furthermore, the larger one is asymptotically stable, the smaller one is unstable, and the origin is also asymptotically stable.
Proof. From Theorem 2, we only need to prove that the origin is asymptotically stable. When u → 0, we haveh(u) → 0 and h(u) → 0, these two facts, combining with (8), (12) and (15) Hence, we obtain which yields, along with Lemma 3, the origin E 0 is asymptotically stable. We complete the proof.

Numerical Simulations
In this section, we give some numerical simulations to support our lemmas and theorems.
Finally, from a different point of view, we give the following numerical example to verify Theorem 3.

Example 2.
Let the parameters be specified in (51). Then, g * = 9506.25 and T * ≈ 14.36. If we select T = 14.2 < T * and c = 8000 < g * or c = g * , then the conditions in Theorem 3 are satisfied, and hence the model (3) and (4) admits exactly two T-periodic solutions. Moreover, the origin E 0 and the larger T-periodic solution are both asymptotically stable, while the other T-periodic solution is unstable, which is shown in Figure 7.

Conclusions
Blood-feeding female mosquitoes are responsible for many life-threatening mosquitoborne diseases, and the most effective way to prevent and control mosquito-borne diseases is to control the vector mosquitoes, and many scholars have been focusing on this project [13,16,32,[41][42][43][44]. The applications of insecticides may rapidly reduce mosquito population size to mitigate the epidemics, however, it also causes environmental pollution and provides only a short-term effect as the insect can develop resistance. SIT and IIT are efficient and environment-friendly methods to control mosquitoes, and they share two common steps: (i) rear a huge number of male mosquitoes in laboratories or factories, and (ii) release them into the field to sterilize wild mosquitoes and thus suppress the wild mosquito population. Considering the cost of mosquito feeding and artificial release, the release project, including release amount and period, requires careful analysis.
Periodic and impulsive phenomena are universal in artificial intervention on biological populations [30,31,33]. The suppression based on SIT or IIT can reduce the wild mosquito population to a low level quickly, but it is difficult to completely eliminate wild mosquitoes and the wild mosquito amount changes periodically in a relatively small range [9]. A periodic solution means that the wild mosquito population is persistent and can be controllable under a low level, which can guarantee that the mosquito population size is below the epidemic risk threshold. In this paper, we assumed that the sterile mosquitoes are released at fixed time points T k = kT(k = 0, 1, 2, . . .) with size c. In addition, we set that the sexual lifespanT of sterile mosquitoes is less than the release period, and the release amount is not larger than g * , which is defined in (5), i.e., T >T and 0 < c ≤ g * . We found that the number of T-periodic solutions depends on the stability of the origin E 0 : if E 0 is unstable, then the model (3) and (4) has a unique globally asymptotically stable T-periodic solution; if E 0 is asymptotically stable, then the model (3) and (4) admits exactly two T-periodic solutions, and the one with a larger initial value is asymptotically stable, while the other is unstable. This manifests the difference between the results in this article and in [38]: the current work implies that the model (3) and (4) has at least one periodic solution, which biologically means that complete suppression depends on both release period and the initial density of wild mosquitoes. Specifically, the wild mosquitoes can be eradicated only when the release period is smaller such that T < T * and the initial size of wild mosquitoes is less than a given threshold value. While [38] shows that E 0 could be globally asymptotically stable, which biologically means that the wild mosquitoes can be eradicated under certain conditions regardless of their initial density. Nevertheless, the release strategy in this paper complements the situation not mentioned in [38]. Combining the conclusions in this paper and in [38], we find that the model (3) and (4) has at most two T-periodic solutions when T >T.
It is well known that mosquitoes undergo four distinct stages (eggs, larvae, pupae, and adults) during their lifetime, and the first three stages occur in water while the final one in air. In our ODE model (3) and (4), we ignored the developmental process and let w(t) denote the number of mosquitoes in the adult stage. However, developmental lags may affect the dynamical behavior of mosquito populations [45]. Thus, it is more reasonable and realistic to take the lag factor into consideration [24,31], which will be harbored in our future work.