Existence and stability of two periodic solutions for an interactive wild and sterile mosquitoes model

In this paper, we study the periodic and stable dynamics of an interactive wild and sterile mosquito population model with density-dependent survival probability. We find a release amount upper bound , depending on the release waiting period T, such that the model has exactly two periodic solutions, with one stable and another unstable, provided that the release amount does not exceed . A numerical example is also given to illustrate our results.


Introduction
Mosquitoes have been listed as one of the deadliest animals around the world, due to their ability of transmitting mosquito-borne diseases through female mosquitoes' bitings. Take dengue as an example, the number of cases reported by World Health Organization increased dramatically over the last two decades, from 505,430 cases in 2000, to over 2.4 million in 2010, and 5.2 million in 2019 [30]. Since there have neither safe vaccines nor effective drugs to combat dengue, biological methods including Sterile Insect Technique (SIT) and Incompatible Insect Technique (IIT) have played vital roles in preventing the transmission of mosquito-borne diseases [1][2][3][4][5]8,11,19,[27][28][29]. Both SIT and IIT involve releasing sterile male mosquitoes into the wild to mate with wild females, which result in no progeny for those females. It occurs that the number of wild mosquitoes will decrease gradually.
Based on a thought-provoking modelling idea proposed in [31] that the number of sterile mosquitoes at time t can be treated as a given function rather than an independent variable constrained by a differential equation, and by taking the fact that density-dependent death of wild mosquitoes mainly occurs in the aquatic stages (eggs, larvae, pupae) [7,24] into consideration, the authors in [21] studied the following ordinary differential equation model for exploring the interactive dynamics of wild and sterile mosquitoes. In model (1), w = w(t) and g = g(t) denote the numbers of wild and sterile mosquitoes at time t, respectively, a is the maximum number of survived offspring produced per individual, w w+g is the mating probability between wild mosquitoes, ξ denotes the carrying capacity parameter such that 1 − ξ w describes the density-dependent survival probability, and μ is the natural mortality of wild mosquitoes.
Although model (1) looks simple, we find that it can generate rich dynamics. Obviously, g(t) plays a key role in determining the asymptotic behaviours of model (1), it is very important for us to make clear the structure of g(t). To this end, letT and T be the sexual lifespan of released sterile mosquitoes and the waiting period between two consecutive releases, respectively. Then there exists three different release strategies: T =T, T >T, and T <T. Moreover, we assume that the release begins at time t = 0 such that g(t) = 0 for t < 0, and a constant amount c of sterile mosquitoes is impulsively and periodically released at discrete time points T i = iT, i = 0, 1, 2, . . .. When T =T, g(t) ≡ c for all t ≥ 0. And model (1) with g(t) ≡ c has been investigated in [21]. For the case when T >T, we have which has recently been discussed in [43], where A = a−μ aξ . We study the release strategy of T <T in this paper. Let [x] represent the largest integer portion of x. Then there exists an integer j = j(T) and a nonnegative number r = r(T) ∈ [0, T) such thatT = jT + r. When t ∈ [0, jT), we get When t ≥ jT, if r = 0, then g(t) = jc, under which model (1) can be reduced to the one in [21]. If r = 0, then See Figure 1 for illustration. Then model (1) reads as We assume 0 < r < T in the remaining of this paper. Since our main concern is the asymptotic dynamics of the model, we may assume that the initial time point is t 0 = jT. Then g(t) takes the form of (2).
Let the release amount upper bound G * be defined as .
Then model (1) can be rewritten as i = j + 1, j + 2, . . .. We, in this paper, mainly study the long-time dynamics of model (3) under the assumptions of T <T and 0 < c ≤ G * . The paper is arranged as follows. In Section 2, we define the Poincaré map and give three lemmas. In Section 3, the theoretical results show that if 0 < c ≤ G * , then the model has exactly two T-periodic solutions, with one stable and another unstable. Finally, a brief discussion is given in Section 4.

Preliminaries
It is apparent that the origin, denoted by w 0 , is the unique equilibrium of model (3). We first give the definition of a solution of model (3). A function w(t) is said to be a solution of model (3) for n = 0, 1, 2, . . ., see [34] for details. Following the lines in [34], for sequences defined in (4), we may similarly obtain the next lemma. By a similar argument to that in the proof of Lemma 2.9 in [32], we arrive at the following lemma, which provides a necessary and sufficient condition for the origin w 0 to be asymptotically stable.

Lemma 2.2: The origin w 0 is asymptotically stable if and only if there is
Since our focus is to find the fixed points of h, we need to know the different relationships between h(u) and u for all u > 0. Moreover, from Lemmas 2.1 and 2.2, we find that the sign of h (u) − 1 is crucial to the comparison between h(u) and u. To obtain the expression for h (u), we will first solve the first equation of model (3) with initial value w(jT) = u for getting the expression ofh (u), then solve the second equation of model (3) with initial value w(jT + r) =h(u) to get the relationship betweenh (u) and h (u). To go ahead without the burden of frequently switching between the first equation and the second equation of model (3), we rewrite model (3) as follows we know that f (w) = 0 has two positive real roots. Denote the two roots by E 1 (k, c) and E 2 (k, c), then Moreover, we have By letting δ 0 = E 1 (j, c) in Lemma 2.2, and using (3) and (6), we have the following lemma.

Lemma 2.3:
The origin w 0 is always locally asymptotically stable.
Biologically, Lemma 2.3 means that the suppression on wild mosquitoes is always successful provided that the initial density of wild mosquitoes is contained within the region of attraction of w 0 .
Recall that the sign of h (u) − 1 is our deepest concern, and to get this sign, it suffices to obtain the expression for h (u). In the following, we are devoted to seeking the expression. To this end, we need to consider the following two cases.
Making the identity transformation for (8), we get Integrating (9) from jT to jT + r with k = j + 1, and recall that we havē For notation simplicity, let Taking the partial derivative with respect to u of both sides of (12), we get Furthermore, substituting (12) and k = j + 1 into (11), we have To get the expression forh (u), we take the derivative with respect to u of both sides of (14), which yields Substituting (13) with k = j + 1 and (14) into (15), we obtain For the solution going fromh(u) to h(u), we integrate (9) from jT + r to (j + 1)T with k = j, which reaches By (12) and k = j, (17) becomes Similarly, to get h (u), we take the derivative with respect to u of both sides of (18), which gives Substituting (12) with k = j and (18) into (19), we get Case 2: c = G * . In this case, Equation (5) is We first observe that the relationship betweenh (u) and h (u) is the same as case 1, i.e. with c = G * , (20) still holds in this case. Thus, for getting the expression for h (u), we only need to know the expression forh (u). To this end, we turn to analyse the next equation By the methods of separation of variables and partial fraction decomposition of the rational function, we obtain where Further computations from (21) offer Integrating (22) from jT to jT + r with k = j + 1, and by (10), we reach Set Then, we have, by calculating the partial derivative of L 2 (u, k) with respect to u, and (23) becomes Differentiating on both sides of (25) with respect to u yields Substituting (24) and (25) into (26) with k = j + 1, we obtain The above preliminaries are sufficient to pave the way for presenting our main results.

Exact two periodic solutions
In this section, we give a sufficient condition for model (3) to admit exactly two periodic solutions, with one stable and the other unstable, which is harboured in the following theorem.
Theorem: Assume that 0 < c ≤ G * . Then model (3) has exactly two T-periodic solutions, with the smaller one unstable and the larger one asymptotically stable.
The theoretical proof of the Theorem will be given after the following three indispensable lemmas are introduced.
Next, we prove the uniqueness of the corresponding T-periodic solutions initiated from (E 1 (j, c), E 1 (j + 1, c)) and (E 2 (j + 1, c), E 2 (j, c)) in the following two lemmas.

Proof: The existence of u
can be deduced from Lemma 3.1. Denote w 1 (t) = w(t; jT, u 1 ). Then w 1 (t) is a Tperiodic solution of model (3). Next, we prove the uniqueness of T-periodic solutions by contradiction.
Assume that there exists u 2 ∈ (u 1 , see Figure 2(B) and (C) for illustration. Then, the growth rate of the Poincaré map h at u i with i = 1, 2 satisfies one of the following two conclusions: which correspond to Figure 2(B) and (C), respectively. To obtain a contradiction, we need to take the derivative of h(u). Set 2 (k, c)) .
Then we have, from (16) and (20), Then at points u i , we have h (u i ) = N(h(u i )) N(u i ) . Since we have N(u i ) < 0, i = 1, 2. From (27), cases (i) and (ii) imply and respectively. However, from (28), we have Since Hence, from (31) and the factsh(u i ) < u i , we have N(h(u i )) < N(u i ), i = 1, 2. A contradiction to (29) and (30). Obviously, (31) can also exclude the possibility of three or more T-periodic solutions to model (3). The proof is complete.

Lemma 3.3:
Denote w 2 (t) = w(t; jT, v 1 ). Then w 2 (t) > w 1 (t), and w 2 (t) is a T-periodic solution of model (3). Hence, we only need to prove the uniqueness of T-periodic solutions to model (3) for u ∈ (E 2 (j + 1, c), E 2 (j, c)). Similarly, assume that model (3) has another T-periodic solution in this case, and we denote it by w(t; jT, v 2 ). Then as illustrated in Figure 3(A) and (B), there are two cases: and and hold, which respectively correspond to cases (i) and (ii) in (32). However, taking the derivative of N(u) yields and hence N (u) < 0 for u ∈ (E 2 (j + 1, c), E 2 (j, c)), which leads to a contradiction to N(h(v 2 )) = N(v 2 ) in (33) and a contradiction to N(h(v 1 )) = N(v 1 ) in (34). Further, N (u) < 0 is sufficient to exclude the possibility of three or more T-periodic solutions to model (3). We complete the proof.
With the above three lemmas, we are now ready to give the detailed proof of the theorem.
We then prove that (II) is also impossible. If not, then we have Furthermore, from (12) and (13), with k = j, we can derive that L 1 u, j > 0, and which manifests that function L 1 (u, j) is strictly decreasing with respect to u. Thus, from (18) and (37), we have that is, we have w((i + 1)T) < w(iT), or, equivalently, h i+1 (u 0 ) < h i (u 0 ), which gives a similar contradiction to that in the proof of case (I). Thus, (II) is also impossible. Finally, we prove (III) is impossible as well. We first prove that the case oft ∈ (iT, iT + r) is impossible. Integrating (9) from iT tot with k = j + 1, we obtain similarly, integrating (9) from (i − 1)T tot − T with k = j + 1, we get from (12) and (13), with k = j + 1, we have which implies that function L 1 (u, j + 1) is strictly increasing with respect to u. Furthermore, since (38) and (39) tell us that L 1 (w(iT), j + 1) < L 1 (w((i − 1)T), j + 1), therefore, follows from the increasing monotonicity of L 1 (u, j + 1) with respect to u, we have w(iT) < w((i − 1)T), a similar contradiction to that in the proof of case (I) occurs again. Thus, the case of t ∈ (iT, iT + r) is impossible.
We then prove the case oft ∈ (iT + r, (i + 1)T) is also impossible. Integrating (9) from t to (i + 1)T with k = j, we have Similarly, integrating (9) fromt − T to iT with k = j, we reach Recall that function L 1 (u, j) is strictly decreasing with respect to u, thus, (40), (41) and (42) give L 1 (w((i + 1)T), j) > L 1 (w(iT), j), hence, we obtain w((i + 1)T) < w(iT). A similar contradiction as that in the proof of case (I) can be achieved. Thus, the case oft ∈ (iT + r, (i + 1)T) is also impossible. Hence, (36) is true, which implies that w 2 (t) is stable. For the attractivity of w 2 (t), we need to prove For the case when u 0 ∈ (u 1 , v 1 ), from (35) this shows that (43) is true for the case when u 0 ∈ (u 1 , v 1 ). For the case when u 0 ∈ (v 1 , +∞), from (35), we have h(u 0 ) < u 0 , the remaining proof is similar to that of u 0 ∈ (u 1 , v 1 ) and is omitted. We complete the proof.
In the following, we will carry out a numerical example to illustrate our theoretical results.

Concluding remarks
Based on the modelling idea in [31] that the number of sterile mosquitoes released is treated as a given nonnegative function rather than a variable satisfying an independent dynamical equation, the interactive dynamics of wild and sterile mosquitoes with or without time delay has been recently studied in [18,21,[32][33][34]41]. A common assumption is that the sterile mosquitoes are impulsively and periodically released at discrete time points kT, k = 0, 1, 2, . . .. Such impulsive and periodic release strategy is consistent with the actual field release situation [42]. Model (1) is derived by combining the ideas in [6] and [31]. We introduce two release amount thresholds g * and c * , and the waiting period threshold T * . In [43], for the case of c ≥ c * , the authors obtained sufficient conditions for the global asymptotic stability of the origin and the existence of a globally asymptotically stable T-periodic solution, respectively. It is well known for obtaining sufficient conditions on the existence of exact two periodic solutions is mathematically challenging. In this paper, for the case of c ∈ (0, G * ], we show that model (1) has exactly two T-periodic solutions and analyse their stability, the bigger one is stable and the smaller one is unstable. The result, combining with that of [43], indicates the better understanding of how to determine a better release strategy for biological control of wild mosquitoes from practical release perspective. Further investigations, for the case of G * < c < G * * , are understudy, this will be a topic of our future study. We hope that our theoretical results can provide real guidance to help the practical workers make the better release strategy to gain a better effect for the wild mosquito population suppression.

Disclosure statement
No potential conflict of interest was reported by the author(s).