Uniqueness and stability of periodic solutions for an interactive wild and Wolbachia-infected male mosquito model

We investigate a mosquito population suppression model, which includes the release of Wolbachia-infected males causing incomplete cytoplasmic incompatibility (CI). The model consists of two sub-equations by considering the density-dependent birth rate of wild mosquitoes. By assuming the release waiting period T is larger than the sexual lifespan of Wolbachia-infected males, we derive four thresholds: the CI intensity threshold , the release amount thresholds and , and the waiting period threshold . From a biological view, we assume throughout the paper. When , we prove the origin is locally asymptotically stable iff , and the model admits a unique T-periodic solution iff , which is globally asymptotically stable. When , we show the origin is globally asymptotically stable iff , and the model has a unique T-periodic solution iff , which is globally asymptotically stable. Our theoretical results are confirmed by numerical simulations.


Introduction
Mosquitoes can spread diseases between humans, or from animals to humans. They ingest pathogens from an infected host and then transmit diseases during subsequent bites. Common mosquito-borne diseases include dengue, chikungunya, Zika, etc. [37]. The World Health Organization has reported that an estimated 100-400 million people are infected with dengue each year [38]. In recent years, dengue has spread rapidly in more than 129 countries [3]. Traditional ways to prevent the transmission of dengue are integrated mosquito management which includes spraying insecticide or removing breeding sites [1,36]. However, besides environmental pollution, spraying insecticide has only a short-term effect on eradicating mosquitoes due to the emergence of insecticide resistance [32,34]. At the same time, removing all breeding sites of larva mosquitoes is an impossible mission. Although the dengue vaccine Dengvaxia was approved for use in 2015, the safety of the vaccine cannot be guaranteed due to the phenomenon of antibody-dependent enhancement [9,11].
Nowadays, an innovative biological control strategy involves an intracellular bacterium called Wolbachia, which exists in about 20% arthropod hosts [19,35]. Wolbachia was first identified in 1924 by Hertig and Wolbach, and has gained widespread attention of researchers after its role in cytoplasmic incompatibility (CI) was revealed by Laven [21] in 1956. CI occurs when Wolbachia-infected males mate with Wolbachia-uninfected females, resulting in eggs produced from mated females will not hatch [15,29]. Since then, scientists have made great efforts to inject Wolbachia into Aedes mosquitoes to offer possible biological tools to combat dengue and other mosquito-borne diseases. The groundbreaking work was credited to Xi et al. [39], who successfully established the first stable Wolbachia strain, named as WB1, in Aedes aegypti by embryonic microinjection in 2005. Subsequent experiments showed that WB1 induced complete CI. Owing to this breakthrough, various Wolbachia strains have been stably established in Aedes albopictus [40], which is the sole vector of dengue in Guangzhou. In 2015, Xi's team built the largest mosquito factory of the world in Guangzhou, China, and a three-year field trial in Shazai island and Dadaosha island by releasing factory-reared male mosquitoes eradicated more than 95% of Aedes albopictus mosquito populations [48].
As the successful establishment of various Wolbachia strains in mosquitoes, mathematical modelling and analysing on the interactive dynamics of wild and released mosquitoes has attracted plenty of attention [5][6][7]22,23,30,33,42,44]. Various mathematical models have been established, including discrete models [4,8,45,47,50], which are suited to cage mosquito populations, ordinary differential equation models [13,20,49] and delay differential equation models [41,43], which consider the overlapping generations in mosquito populations and the maturation delay of mosquitoes from egg to reproductive adults, respectively. Also, stochastic models [16,17] have been developed to include the random switch of environments, and reaction diffusion models [18,28,31] have been investigated to include the effect of mosquito dispersal on the interactive dynamics of wild and released mosquitoes.
Among most of the above models, the authors assumed that CI is complete, i.e. no eggs from the incompatible mating will hatch. Complete CI is the most ideal situation when we aim to sterilize wild females through CI by releasing Wolbachia-infected males. However, experiments showed that CI could be incomplete [10,40] due to different hosts or different Wolbachia strains, which could also be influenced by environmental conditions. In this paper, we aim to develop a mosquito suppression model which includes both the release of male mosquitoes and the incomplete CI mechanism. To this end, we introduce the first parameter s h ∈ (0, 1) to estimate the CI intensity. That is, s h is the relative reduction in egg hatch caused by incompatible crosses, and the hatch rate of the compatible case is s h = 1. Let w(t) and g(t) be the numbers of wild and Wolbachia-infected mosquitoes at time t, respectively. Then the ratio g(t)/(w(t) + g(t)) represents the probability of wild mosquitoes mating with Wolbachia-infected mosquitoes under the assumption of random mating behaviour [8]. Let a be the total number of offspring per wild mosquito per unit time. With the interference of Wolbachia-infected males, the birth rate of wild mosquitoes is reduced from a to a · w(t) .
Meanwhile, mosquitoes go through four stages of complete metamorphosis: egg, pupa, larva, and adult, with the first three stages named as aquatic stages. Based on the experimental observations that the crowding effect basically takes place in aquatic stages [12,14], many authors included the density-dependent effect in the growth of the mosquito population [24,26]. Under this consideration, we also assume that the density dependence (denoted by ξ ) of mosquitoes development process takes place in aquatic stages. Thus, the birth rate of wild mosquitoes is Together with the assumption that the death rate of wild mosquitoes is densityindependent, we reach the following model, where μ is the death rate of wild mosquitoes, and we assume that a > μ.
We should mention here that the modelling idea in our model which treated the released male mosquitoes as a given function g(t) was first proposed by Yu [41]. In [41], Yu developed an effective strategy to suppress wild mosquitoes by releasing Wolbachia-infected male mosquitoes, which is governed by a given function g(t). This strategy reduces the difficulties in theoretical researches by reducing dimensions of systems. Since then, this novel modelling idea has attracted a great attention of many scholars. Yu introduced a parame-terT as the sexual activity lifespan of Wolbachia-infected males. The release function g(t) is then determined by the relation ofT with the waiting period T between two consecutive releases. We assume that the Wolbachia-infected mosquitoes are released impulsively and periodically at discrete time points T k = kT, k = 0, 1, 2, . . ., and the release amount of each batch is maintained at a constant level c. The simplest case of g(t) happens when T =T, that is, once the released males lose their mating ability, another batch of infected males is released supplementarily. Thus, we have the release function and Model (1) takes the form We provide a complete dynamics of Model (2) in Section 2, which is the basis for the discussion of T >T and T <T. We only focus on the case T >T in this paper and leave the case T <T in future work. Thus, the release function g(t) becomes and Model (1) can be rewritten as two sub-equations where A = (a − μ)/aξ .
In this study, we formulate a mosquito population suppression model interfered by Wolbachia-infected males causing incomplete CI. In Section 2, the local and global dynamics of wild mosquitoes are studied completely when T =T. In Section 3, by using Poincaré map, we obtain the sufficient and necessary conditions for the stability of the origin E 0 and the existence and uniqueness of a globally asymptotically stable T-periodic solution for the case T >T, respectively. Finally, some numerical examples are provided to confirm our main results in Section 4.

The dynamical behaviour for the case T =T
In this section, we focus on a special release strategy, that is, the waiting period T between two consecutive releases is equal to the sexual activity lifespanT. The release amount of each batch is controlled at stable level c, and the loss of released mosquitoes during the short sexual activity lifespan can be ignored. Thus, the total amount of released mosquitoes with sexual activity is maintained at the level c. The dynamics of wild mosquitoes ω(t) are mainly governed by (2). Thus, to better understand the dynamical behaviours, we need to count the positive equilibria of (2) and analyse their stability.
We claim that Model (2) has no positive equilibrium when the release amount c ≥ 1/ξ . In fact, from (2) we easy to find that the part in the first bracket is positive and the part in the second bracket is negative when c ≥ 1/ξ . Taking together, we derive that the right part of (2) is negative for all ω(t) > 0 implying that (2) has no positive equilibrium. In the following, we only need to consider the case that 0 < c < 1/ξ and analyse the existence of the positive equilibria over the rectangular area in s h -c plane, as shown in Figure 1. We where The equilibria of Model (5) satisfy −wf (w)/(w + c) = 0. It is very obvious that Model (5) has at least one trivial equilibrium E 0 = 0. We then search other equilibria from f (w) = 0. Since f (w) = 0 is a quadratic equation, its discriminant is 2 , which is a quadratic polynomial with respect to c. It is easy to find that there are two positive real roots c 1 and c 2 such that = 0, where c 1 and c 2 are given as The two roots divide (0, ∞) into three intervals: (0, c 1 ], (c 1 , c 2 ) and [c 2 , ∞). We compare c 2 with 1/ξ and find which implies that Model (5) has no positive equilibrium when the release amount c ∈ [c 2 , ∞). When c ∈ (c 1 , c 2 ), the discriminant < 0 holds, then f (w) is positive for all w > 0, thus Model (5) also has no positive equilibrium. We only need to find the positive equilibria for c ∈ (0, c 1 ]. When c = c 1 , we have = 0 and derive a unique equilibrium of (5), that is, From the expression of c 1 , it is a function of s h when s h varies over (0, 1). This curve divides the rectangle in Figure 1 into two areas. For the case that the release amount c = c 1 , the equilibrium E * is positive only when When c ∈ (0, c 1 ), Model (5) has two equilibria, that is, Now, we need to determine the signs of the two equilibria. Since f (w) is a convex function, it is easy to get that f (w) = 0 has a unique positive root when f (0) < 0, or f (0) = 0 and aξ c(2 − s h ) − a + μ > 0. We consider f (0) = 0 and derive another curve which is tangent to c 1 at s * h . The two curves divide the rectangle in Figure 1 into six areas. When c < c 0 , we have f (0) < 0 and Model (5) has a unique equilibrium E 2 . When c = c 0 , Model (5) has a unique positive equilibrium E 2 if s h > s * h and has no positive equilibrium if s h ≤ s * h . By further analysing, we find f (0) > 0 when c > c 0 , and Model (5) has two positive equilibria (1) If c ∈ (0, c 0 ), then Model (5) has two non-negative equilibria: the origin E 0 and the unique positive equilibrium E 2 , where E 0 is unstable and E 2 is asymptotically stable. (5) has only one non-negative equilibrium E 0 , and E 0 is globally asymptotically stable.

Theorem 2.2:
For the case s h ∈ (s * h , 1), the following conclusions hold.
(1) If c ∈ (0, c 0 ], then Model (5) has two non-negative equilibria: the origin E 0 and the unique positive equilibrium E 2 , where E 0 is unstable and E 2 is asymptotically stable. (5) has three non-negative equilibria: the origin E 0 , two positive equilibria E 1 and E 2 , where both E 0 and E 2 are asymptotically stable, and E 1 is unstable.
(3) If c = c 1 , then Model (5) has two non-negative equilibria: the origin E 0 and the unique positive equilibrium E * , where E 0 is asymptotically stable and E * is semi-stable, namely, the right-side of E * is stable and the left-side is not. (4) If c ∈ (c 1 , 1/ξ ), then Model (5) has only one equilibrium E 0 , and E 0 is globally asymptotically stable.
Proof: The existence of equilibria in each domain can be derived from the above analysis, and we only need to prove the stability. We consider the most complicated case, that is, s h ∈ (s * h , 1) and c ∈ (c m , c 1 ). In this case, Model (5) has three non-negative equilibria E 0 , E 1 and E 2 , hence (5) can be rewritten as The interval (0, ∞) is divided into three intervals: which implies that E 0 and E 2 are asymptotically stable, and E 1 is unstable. The proof of Theorem 2.2(2) is completed. For other cases, the proof is similar and omitted here.
Furthermore, combining Theorems 2.1 and 2.2, we can obtain the CI intensity threshold s * h and the release amount threshold (see Figure 1), In fact, what we most like to see is that E 0 is globally asymptotically stable, which means that no matter what the initial number of wild mosquitoes is, the population will be eventually eradicated in real life. The above theorems state that wild mosquitoes will tend to be extinct when the release amount of Wolbachia-infected males is larger than g * , otherwise, the mosquito population can only be suppressed locally.
On this basis, we continue to study the case T >T. When T >T and c > g * , the global dynamics behaviour of Model (1) will be quite different from the case T =T. For biological relevance, we only consider s h ∈ (s * h , 1) in the following.

The dynamical behaviour for the case T >T
Since the release period T is greater than the sexual activity lifespanT of Wolbachiainfected males, we easy to find that g(t) is a periodic piecewise function and jumps between two constants c and 0. The release strategy makes the time evolution of w(t) switch between (3) and (4). Then Model (1) has a complex dynamical behaviour.

The Poincaré map of model (1)
For the case T >T, the dynamical behaviour of (1) mainly depends on the constantly switching between (3) and (4). According to (3) and (4), Model (1) has only one equilibrium E 0 . Let w(t) = w(t; 0, u) be a solution of (1) with the initial condition w(0) = u > 0. Clearly, w(t) is continuous and piecewise differentiable. Initial values satisfy w(0) = u on [0,T) for (3) and w(T) = w(T − ) on [T, T) for (4). In reality, w(T; 0, u) and w(T; 0, u) are continuously differentiable with respect to u and w(T), respectively, where w(T; 0, u) can be regarded as a solution of the initial value problem of (4) with w(T) = w(T − ). The solution can be defined in the same way in other intervals. Furthermore, throughout this paper. For Equations (3) and (4), we introduce two new thresholds which are the release amount threshold c * and the waiting period threshold T * , that is, and Since c * is strictly increasing with respect to s h and c * = g * when s h = s * h , we have c * > g * for s h ∈ (s * h , 1). Next, we study the initial value problem of Equation (3) with w(0) = u on [0,T) and Equation (4) where Equation (11) is equivalent to where Integrating (12) from 0 toT, we obtain We rewrite (4) as Integrating it fromT to T, we haveh where m = e −Aaξ(T−T) = e −(a−μ)(T−T) . Here h(u) is implicitly determined by (13) and (14). Following (13) and (14), we have sinceh(u) → 0 and h(u) → 0 as u → 0. Thus, we have According to the Poincaré map of Equations (3) and (4) By using a similar discussion as in [44], we obtain the following lemma.

Release amount c ∈ (g * , c * )
In this case, the release amount c of mosquitoes is controlled between g * and c * . To study the dynamical behaviour, we establish two lemmas before giving the main result.

Lemma 3.2:
If T > T * and c > g * hold, then (1) has a unique T-periodic solution.
Proof: First, we show the existence of T-periodic solution. According to (16), we have Thus, there exists a sufficiently small δ 0 > 0 such that According to (1), we know that the solution w(t) = w(t; 0, A) is strictly decreasing on [0,T), that is, we haveh(A) < A, which means that h(A) < A though w(t) = w(t; 0, A) is strictly increasing on [T, T). Hence, there must exist a u 1 ∈ [δ 0 , A) such that h(u 1 ) = u 1 , h (u 1 ) ≤ 1, and h(u) > u, for u ∈ (0, u 1 ).
Next, the uniqueness of the T-periodic solution is showed by contradiction. Assume that there exists another T-periodic solution u 2 in (u 1 Without loss of the generality, we assume that u 1 and u 2 are the minimal and maximal Tperiodic solutions on [δ 0 , A), respectively. There must exist av ∈ [u 1 , u 2 ] such that h(v) = v due to the continuous differentiability of h(u). From the fact that 0 < u 1 ≤v ≤ u 2 < A, there are two cases:v ∈ (u 1 , u 2 ), orv = u 1 (or u 2 ). For the casev ∈ (u 1 , u 2 ), let Differentiating P(u) with respect to u gives · P(u).
Substituting (21) into (13), we obtain Taking the derivative with respect to u in (22), we have At the same time, taking the derivative of Equation (14), we havē By (14) again, we geth .
(26) According to (19), (20) and (26), we have which is equivalent to which is a contradiction. For the casev = u 1 or u 2 , without loss of generality, we assume thatv = u 1 . Then, we have and h(u) < u, for u ∈ (u 2 , A).
Let k be a small constant such that k > 1, mk < 1, and h(u) − ku has three roots According to (26), we have Equation (28) is equivalent to is in contradiction to the concavity of Q k (u). The proof is completed.

Lemma 3.3: If T = T * and g * < c < c * hold, then (1) has a unique T-periodic solution.
Proof: To prove the existence of T-periodic solution, similar to the proof of Lemma 3.2, we only need to show (18) is true. Applying (13) and (14), we have .
Since m = e −aξT/α when T = T * , we obtain , that is, Set
Thus, H(u) is strictly increasing for u ∈ (0,ū), that is, H(u) > H(h(u)). Therefore, there exists a u ∈ (0,ū) such that h(u) > u. Then Model (1) has at least one T-periodic solution for u ∈ (0, A) with the fact that h(u) < u for u ≥ A. Now, we show the uniqueness of the T-periodic solution. From the existence of Tperiodic solution, we choose a v 1 ∈ (ū, A) such that

If (1) has another T-periodic solution, then there exists a
According to (30)
Obviously, w(nT; 0, u) is the maximum of w(t) for t ∈ [(n − 1)T +T, nT +T], n = 1, 2, . . ., hence we only need to prove that By (32) and Lemma 3.1, we know that the sequence {h n (u)} is strictly decreasing for u ∈ (0, δ). Therefore, the limit lim n→∞ h n (u) = l exists. Choosing l = 0, otherwise there will exist a T-periodic solution w(t; 0, l). In conclusion, E 0 is locally asymptotically stable. Now, we give the proof of (2). Assume that Model (1) has a unique globally asymptotically stable T-periodic solution. If T < T * , then there is a contradiction with the above conclusion. Hence, we have T ≥ T * .
Inversely, we assume that T ≥ T * and show that Model (1) has a unique globally asymptotically stable T-periodic solution. Employing Lemmas 3.2 and 3.3, we know that Model (1) has a unique T-periodic solution. Define the T-periodic solution as w(t) = w(t; 0,ũ), forũ ∈ (0, A). Obviously, We first show the stability of the T-periodic solution. For any ε > 0, t 0 ≥ 0, select δ = ε, it suffices to show that |u −ũ| < δ, that is, Assume that (34) is not true. Then there must exist at > t 0 such that Without loss of generality, we assume that there must be a nonnegative integer p such thatt = pT, ort = pT +T, ort ∈ (pT, pT + T) ∪ (pT +T, (p + 1)T). Assume that t 0 ≤ (p − 1)T. Otherwise, the solutions w(t) and w(t) can reverse to the point t = (p − 1)T. For the caset = pT, we get w(pT) =w(pT) + ε, and w(t) <w(t) + ε, for t ∈ [t 0 , pT), which means that From Lemma 3.1, we know that sequence {w(nT)} is strictly increasing and lim n→∞ w(nT) = l where l ∈ (0, A). Hence, w(t; 0, l) is another T-periodic solution of Model (1), which leads to a contradiction. For the caset = pT +T, we obtain w(pT +T) =w(pT +T) + ε, and w(t) <w(t) + ε, for t ∈ [t 0 , pT +T).
By (14), we have which implies that the sequence {w(nT)} is strictly increasing by Lemma 3.1. Thus, there is another T-periodic solution, a contradiction again.

Release amount c ≥ c *
Before giving the main result, we first give a lemma for the case c ≥ c * .

Lemma 3.4:
If T ≤ T * and c ≥ c * hold, then h(u) < u for any u > 0.

Proof:
The proof here is divided into two cases: T = T * and T < T * . For the case T = T * , by Lemma 3.2 and c ≥ c * , we haveū ≤ 0, so H (u) < 0 for u ∈ (0, A). Therefore, H(u) is strictly decreasing. Furthermore, from (29), we obtain Thus, we have h(u) < u for u ∈ (0, A).
In conclusion, with the fact h(u) < u for u ≥ A, h(u) < u for u > 0 is true under the given conditions. The proof, therefore, is completed. Proof: (1) Assume that E 0 is globally asymptotically stable. Therefore, to prove T ≤ T * , it suffices to show that h (0) ≤ 1. Similar to the proof of Theorem 3.1(1), the conclusion can be proved by contradiction. Assume that T ≤ T * . The stability can be shown in the same way as Theorem 3.1(1). Using Lemma 3.4, we have h(u) < u for any u > 0, which means that the sequence {h n (u)} is strictly decreasing. Hence, we have lim n→∞ h n (u) = 0, same as the proof of Theorem 3.1(1), which implies that E 0 is globally attractive. Thus, the origin E 0 is globally asymptotically stable.
(2) Assume that the T-periodic solution of Model (1) is unique globally asymptotically stable. If T ≤ T * , according to the above result, we know that E 0 is globally asymptotically stable, which contradicts the assumption. Assume that T > T * . By Lemma 3.2, we have a unique T-periodic solution. Similar to the proof of Theorem 3.1(2), the global asymptotic stability of T-periodic solution can be proved. The proof is completed. Proof: Assume that Model (1) has a unique globally asymptotically stable T-periodic solution denoted by w * (t) = w(t; 0, u * ) with u * ∈ (0, A). We prove that E 0 is unstable by contradiction. If E 0 is stable, by the definition of stability, we know that the solution w(t; t 0 , u 0 ) does not approach w * (t), which is a contradiction. In turn, assume that E 0 is unstable. Similar to the proof of Theorems 3.1(2) and 3.2(2), we can show that Model (1) has a unique globally asymptotically stable T-periodic solution. Detailed steps are omitted here.

Numerical simulations
According to the data in [2,25,27,46] and the method of choosing data in [49], we get the values of parameters as  The most difficult parameter to determine is the density-dependent death rate ξ . Note that ξ is primarily determined by the environmental conditions of the area studied. Normally, when the environment is suitable for mosquitoes breeding, ξ will be very small. Not to change the dynamics significantly, we can take ξ = 0.001 as a typical example. We provide two numerical simulations to confirm the analytical results of Theorems 3.1 and 3.2. We choose s h = 0.9 > s * h and g * < c = 790 < c * . The release period is set atT < T = 14.04 < T * . As shown in Figure 2(a), when the initial amount is small, the wild mosquitoes can be suppressed finally, implying that the origin E 0 is locally asymptotically stable, as stated by Theorem 3.1(1). If the release period T is set at T = T * or T = 15 > T * , the amount of wild mosquitoes oscillates periodically around some level, as shown in Figures 2(b) and 2(c), which verifies the result given in Theorem 3.1(2). And we can find the amplitude increases with the release period. Example 4.2: Same as the above example, we select the same parameter values a, μ, ξ and T, and get four thresholds s * h , g * , c * , and T * . Choosing s h = 0.9 > s * h and c = c * , if the release periodT < T = 14.01 < T * orT < T = T * , then the origin E 0 is globally asymptotically stable as shown in Figures 3(a) and 3(b). If T = 15 > T * , then the amount of wild mosquitoes oscillates periodically as shown in Figure 3(c). Next, we choose c = 900 > c * and selectT < T = 14.01 < T * ,T < T = T * or T = 15 > T * , respectively. The same conclusions are reached (see Figure 3(d,e,f)).

Conclusion
Based on the modelling methods of Yu [41] and Li [24], and the incomplete CI mechanism, we considered a population suppression model consisting of two sub-equations switching each other in this study by releasing Wolbachia-infected male mosquitoes. In modelling, we introduced the CI intensity s h ∈ (0, 1) and the sexual lifespanT of Wolbachia-infected mosquitoes, and ignored the independent dynamical equation for Wolbachia-infected mosquitoes.
We assumed that the release amount g(t) of Wolbachia-infected mosquitoes was periodic and impulsive at discrete time points T k = kT, k = 0, 1, 2, . . ., and focused on the case that the waiting period T was greater than the sexual lifespanT. In this paper, we defined four thresholds: the CI intensity threshold s * h , the release amount thresholds g * and c * , and the release waiting period threshold T * . We obtained a relatively complete analysis for the dynamics of Model (1) when CI intensity s h is greater than its threshold s * h . Several important lemmas were established and used to prove the main results. For many dynamical systems, it is challenging to obtain the conditions for the existence and uniqueness of periodic solutions. We established necessary and sufficient conditions for both the stability of the origin E 0 and the existence and uniqueness of the globally asymptotically stable T-periodic solution in Theorems 3.1 and 3.2, respectively.
In Theorem 3.1, we showed that when the conditions s h > s * h and g * < c < c * hold, Model (1) has either a locally asymptotically stable origin E 0 or a unique globally asymptotically stable T-periodic solution depending on T < T * or T ≥ T * . That is, if we don't release enough Wolbachia-infected mosquitoes making the release amount be less than the threshold c * , then whether the mosquito population is successfully suppressed depends on its initial value even though we release at a sufficient frequency T < T * (see Figure 2(a)). At the same time, when the mosquitoes are released infrequently, namely, T ≥ T * , then the suppression of mosquito population is unsuccessful (see Figure 2(b,c)). In Theorem 3.2, we obtained that when the conditions s h > s * h and c ≥ c * hold, the globally asymptotically stable origin E 0 or the unique globally asymptotically stable T-periodic solution of Model (1) are obtained if T ≤ T * or T > T * . In other words, if we release enough Wolbachia-infected mosquitoes (c ≥ c * ) and release them frequently (T ≤ T * ), then the mosquito population will be suppressed (see Figure 3(a,b,e,d)). On the contrary, if the frequency of release is not frequent (T > T * ), then the mosquito population is persistent (see Figure 3(c,f)). Combining Theorems 3.1 and 3.2, under the conditions s h > s * h and c > g * , the necessary and sufficient condition for the existence of the unique globally asymptotically stable T-periodic solution is that the origin E 0 is unstable.
In summary, we studied the dynamical behaviours of the wild mosquito population when s h > s * h and c > g * , and gave complete conclusions and detailed steps. These studies have reference values for guiding mosquito control. The cases of s h ≤ s * h or c ≤ g * have not been studied, and these investigations are left to future work.

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

Funding
This work was supported by the National Natural Science Foundation of China [grant number 12171113].