Dynamic analysis of wild and sterile mosquito release model with Poincaré map

1.   Introduction

Pulsed semi-dynamic systems are widely used in the study of biological systems with thresholds, such as biological resource management ^[1,2,3,4,5] and control of epidemics^{[6,7,8,9,10,11]}, etc. These systems include interactions of continuous and discrete dynamics, where the portion that exhibits discontinuity is called a pulse set, thus producing pulse dynamics. In recent years, the qualitative theory of pulsed semi-dynamic systems has been extensively developed^{[1,12,13,14,15]}, using a variety of analytical methods, such as successor functions, Bendixson theorem, etc^{[16,17,18,19]}. However, due to the complexity of state-dependent feedback control, the global dynamics of the state-dependent feedback control model is far from being solved. In order to better study the qualitative theory and dynamic complexity of the state-dependent feedback control model, we use the Poincar$\acute{\mbox{e}}$ map method to analyze and discuss the model.

The most classic and commonly used method for pest control is chemical control, which involves spraying insecticides to control pests. Due to the widespread use of this method, many pests have developed resistance to insecticides^{[20,21,22,23]}. Especially in the summer, when mosquito population get rampant, repeated use of insecticides reduces their effectiveness and control. It takes a long time to use only biological control methods, and factors such as temperature and weather may affect biological control^{[24,25,26,27]}. When the number of mosquitoes suddenly increases, it is impossible to control the mosquito population in a timely and effective manner. Therefore, in order to solve these problems, we adopted integrated mosquito control^{[28,29,30,31]}, that is, the use of both chemical control methods and biological control methods for wild mosquitoes.

The biological control method used in this paper is the sterile insect technique (SIT). In this way, mosquitoes are disturbed by the natural reproduction process, so that the wild mosquitoes is reduced ^{[32,33,34,35]}. Insect sterility technology uses some physical methods to make male mosquitoes sterile and unable to produce offspring. And these sterile mosquitoes are released into the environment to mate with wild females but not to reproduce. This repeated release of sterile mosquitoes or the release of large numbers of sterile mosquitoes may control wild mosquito population^{[32,36,37,38]}.

In terms of release methods, many articles used proportional release or continuous release^{[39,40,41,42]}, and these two release methods have some advantages only when the wild mosquitoes population is small. In the actual environment, the mosquito population may be very large. In order to control wild mosquito population more effectively, we adopt a new release method in which the release of sterile mosquitoes is in proportional to the number of wild mosquitoes when the number of wild mosquitoes is small, and while the release of sterile mosquitoes saturates and approaches a constant when the number of wild mosquitoes increases.

The structure of this paper is as follows: In Section 2, we briefly introduce the model and make qualitative analysis. In Section 3, the definition domain and main properties of Poincar$\acute{\mbox{e}}$ map are introduced and proved. In Section 4, the existence and stability of periodic solutions are analyzed and proved. In Section 5, we simulate the model to verify our conclusion. Finally, the article is analyzed and summarized.

2.   Model formulation

2.1.   Wild and sterile mosquitoes release model

In article ^[43], the release model of wild and sterile mosquitoes using only biological control is proposed as follows:

where $w(t)$ and $g(t)$ indicate the population densities of wild and sterile mosquitoes at time $t$ respectively. $r$ represents the wild mosquito population birth rate and $b$ represents the sterile mosquitoes release rate. Wild and sterile mosquito populations follow logistic growth without interaction then $a_{i}$ and $\varphi_{i}$, i = 1, 2, denote the density-dependent and density-independent death rates respectively.

The model (2.1) using only biological control, is more susceptible to external influences, and cannot control sudden growth or a large number of wild mosquito population in a timely and effective manner. Therefore, we propose the following wild and sterile mosquitoes release model using integrated mosquito control:

where $Z_{1}(D), Z_{2}(D)$ represent the survival part of wild and sterile mosquitoes when $D$ dose of insecticide is administered respectively, and $0\leq Z_{1}(D)\leq 1, 0\leq Z_{2}(D)\leq 1$. While killing wild mosquitoes with insecticides, we inevitably kill a certain number of sterile mosquitoes. Therefore we use the response curves of the two populations for insecticides to be expressed in $Z_{1}(D) = e^{-k_{1}D}$ and $Z_{2}(D) = e^{-k_{2}D}$ respectively ^[44], where $0\leq k_{i}\leq 1$ represents the pharmacokinetics of insecticides. We release sterile mosquitoes in a ratio-dependent manner, when the density of wild mosquitoes reaches the threshold, sterile mosquitoes are released at a rate of $\frac{bw}{1+w}$. Define $w(0^{+})$ and $g(0^{+})$ as the initial density of wild mosquitoes and sterile mosquitoes, and assumed the initial density of wild mosquitoes is less than the threshold value of $ET$. We can set $ET < \frac{r-\varphi_{1}}{a_{1}}$.

2.2.   The qualitative analysis

In the absence of a pulse, system (2.2) becomes the following:

By caculation, we can easily get system (2.3) with only one equilibria point $E(w_{0}, 0)$, where $w_{0} = \frac{r-\varphi_{1}}{a_{1}}$. Define the isocline of model (2.3) as follows:

For model (2.3), the following theorem is satisfied:

Theorem 2.1. When $r-\varphi_{1}-2a_{1}w < 0$, axial equilibrium point $E(w_{0}, 0)$ is locally asymptotically stable (Fig. 1).

Proof. Let

By calculation, we can obtain

So the jacobian matrix at $E(w_{0}, 0)$ point is

From the above results, we obtain that when $r-\varphi_{1}-2a_{1}w < 0,$

then $E(w_{0}, 0)$ is the point of locally asymptotically stable. This completes the proof.

The following article is discussed under the condition of $r-\varphi_{1}-2a_{1}w < 0$ and $b > 0$.

3.   The definition and properties of the Poincar$\acute{\mbox{e}}$ map

First we give the definition of Poincar$\acute{\mbox{e}}$ map, which is used to study the dynamics of model (2.2). Define the following two lines

By bringing the $w = ET$ into a straight line $L_{1}$, get the intersection of $L_{1}$ and $L_{3}$, defined as $M_{1}(ET, g_{M})$, where

similarly, the intersection of $L_{1}$ and $L_{2}$ is defined as $N_{1}(e^{-k_{1}D}ET, g_{N})$, where

The set is defined in $R_{+}^{2}$ as follows

Now we define $M$ as the impulse set for model (2.2) $M = \{(w, g)\in\Omega|x = ET, 0 < g < g_{M}\}$ and M is a closed subset of $R_{+}^{2}$.

While the continuous function is as follows

Then define the phase set as

Without losing generality, we suppose that the initial point $(w_{0}^{+}, g_{0}^{+})$ belongs to the phase set.

Define the following two sections:

Choose $W_{2}$ as the Poincar$\acute{\mbox{e}}$ section. Then the trajectory starts at point $S_{k}^{+}(e^{-k_{1}D}ET, g_{k}^{+})\in L_{2}$ and intersects with $L_{3}$ at the unique point $S_{k+1}(ET, g_{k+1})$. We get $g_{k+1}$ determined by $g_{k}^{+}$ and $g_{k+1} = \sigma(g_{k}^{+})$ from the Cauchy-Lipschitz theorem. Point $S_{k+1}$ goes through one pulse and then it reaches $L_{2}$ at point $S_{k+1}^{+}(e^{-k_{1}D}ET, g_{k+1}^{+})$ with $g_{k+1}^{+} = e^{-k_{2}D}g_{k+1}+\frac{bET}{1+ET}$. So we define the Poincar$\acute{\mbox{e}}$ section as follows

In order to better study the dynamic behavior of system (2.2), we can propose a Poincar$\acute{\mbox{e}}$ map determined by the pulse point in the phase set according to the definition of the phase diagram. Therefore, we define

It satisfies a scalar differential equation in a phase space

For model (3.2), we only focus on the regions

where $h(w, g)$ is continuous and differentiable. And then we define $w_{0}^{+} = e^{-k_{1}D}ET, g_{0}^{+} = X, X\in N$. Define

and from model (3.2), we have

Thus, the Poincar$\acute{\mbox{e}}$ map $P(X)$ can be represented by the following expression

Next, we prove the main properties of the Poincar$\acute{\mbox{e}}$ map.

Theorem 3.1. Poincar$\acute{\mathit{\mbox{e}}}$ map $P$ has the following properties (Fig. 2)

(I) The domain of $P$ is $[0, +\infty)$. And the Poincar$\acute{\mathit{\mbox{e}}}$ map $P$ is increasing on $[0, g_{N}]$, decreasing on $[g_{N}, +\infty)$.

(II) Poincar$\acute{\mathit{\mbox{e}}}$ map $P$ is continuously differentiable.

(III)Poincar$\acute{\mathit{\mbox{e}}}$ map $P$ has a unique fixed point $\overline g$.

(IV) When $g_{k}^{+}\rightarrow +\infty$, the Poincar$\acute{\mathit{\mbox{e}}}$ map $P$ is bounded and there exists horizontal asymptote $g = \frac{bET}{1+ET}$.

(V)Poincar$\acute{\mathit{\mbox{e}}}$ map $P$ takes the maximum value at $g = g_{N}$ and the maximum value is $e^{-k_{2}D}g(e^{-k_{1}D}ET, g_{N})+\frac{bET}{1+ET}$; the minimum value is taken at $g = 0$, and the minimum value is $\frac{bET}{1+ET}$.

Proof. (Ⅰ) Through qualitative analysis of model (2.3), we get that all trajectories tend to point $E_{0}$. After having a pulse, all the trajectories from the phase set eventually reach the pulse set. Therefore, the domain of $P$ is $[0, +\infty)$.

For any point $P_{1}^{+}(e^{-k_{1}D}ET, g_{p_{1}}^{+})$, $P_{2}^{+}(e^{-k_{1}D}ET, g_{p_{2}}^{+})$ with $g_{p_{1}}^{+}, g_{p_{2}}^{+} \in [0, g_{N}]$, we assume that $g_{p_{1}}^{+} < g_{p_{2}}^{+}$ and $g_{p+1} = \sigma(g_{p})$. From the Cauchy-Lipschitz Theorem, $g_{p_{1}+1} < g_{p_{2}+1}$ can be obtained. Therefore after one pulse, we get

Thus $P$ is increasing on $[0, g_{N}]$.

Similarly, we arbitrarily select two points $Q_{1}^{+}(e^{-k_{1}D}ET, g_{q_{1}}^{+})$, $Q_{2}^{+}(e^{-k_{1}D}ET, g_{q_{2}}^{+})$ with $g_{q_{1}}^{+}, g_{q_{2}}^{+} \in [g_{N}, +\infty)$ and assume that $g_{q_{1}}^{+} < g_{q_{2}}^{+}$. The trajectory starts from point $Q_{1}^{+}$, $Q_{2}^{+}$ across the $L_{2}$ and then hits the $L_{3}$ at points $(ET, g_{q_{1}+1})$ and $(ET, g_{q_{2}+1})$, where $g_{q_{1}+1} > g_{q_{2}+1}$. After the pulse, we get

Thus $P$ is decreasing on the $[g_{N}, +\infty)$.

(Ⅱ) From model (2.3), we get that $F(w, g), G(w, g)$ are continuous and differentiable in the first quadrant. The continuity and differentiability theorems of solutions of ordinary differential equations with respect to their initial conditions are used to determine the differentiability of $P$. That is, $P$ is continuously differentiable in $\Omega$ according to the Cauchy-Lipschitz theorem with parameters.

(Ⅲ) Since $P$ decreases on $[g_{N}, +\infty)$, therefore their exists a $g^{'}\in[g_{N}, +\infty)$ such that $P(g^{'}) < g^{'}$. We can easily obtain $P(0) = \frac{bET}{1+ET} > 0$, then there exists $\overline g\in(0, g^{'})$ such that $P(\overline g) = \overline g$. Therefore $P$ has fixed point on $[0, +\infty)$.

Next we prove the uniqueness of the fixed point. When $P(g_{N}) < g_{N}$, the intersection of the phase set and the $w$ axis is defined as $n_{0}(e^{-k_{1}D}ET, 0)$, then point $n_{0}$ reaches point $(e^{-k_{1}D}ET, g_{n_{0}+1}^{+})$ after one pulse. And the intersection point $N_{1}(e^{-k_{1}D}ET, g_{N})$ of $L_{1}$ and $L_{2}$ reaches the point $(e^{-k_{1}D}ET, g_{N+1}^{+})$ after one pulse. We assume that $P$ has two fixed points $G_{1}(e^{-k_{1}D}ET, \overline{g}_{1}), G_{2}(e^{-k_{1}D}ET, \overline{g}_{2})$ and $\overline{g}_{1}, \overline{g}_{2}\in(g_{n_{0}+1}^{+}, g_{N+1}^{+})$ such that $P(\overline{g}_{1}) = \overline{g}_{1}, P(\overline{g}_{2}) = \overline{g}_{2}$. Since $P(0) = \frac{bET}{1+ET} > 0$, then by the definition of the Poincar$\acute{\mbox{e}}$ map, we get $P(g) > g$ for $g\in(g_{n_{0}+1}^{+}, \overline{g}_{1})$. From $P(\overline{g}_{1}) = \overline{g}_{1}$, we get for any $g_{1}^{'}\in(\overline{g}_{1}, \overline{g}_{2})$, $P(g_{1}^{'}) < g_{1}^{'}$ holds. And $P(g_{2}^{'}) < g_{2}^{'}$ for $g_{2}^{'}\in(\overline{g}_{2}, g_{N+1}^{+})$ due to $P(g_{N}) < g_{N}$ and $P(\overline{g}_{2}) = \overline{g}_{2}$. Therefore we obtain $g\in(\overline{g}_{1}, g_{N+1}^{+})$ such that $P(g) < g$, which is a contradiction. So $P$ has a unique fixed point $\overline g$ and $\overline g\in(0, g_{N})$.

Similarly, when $P(g_{N}) > g_{N}$, we get that $P$ has unique fixed point $\overline g$ and $\overline g\in(g_{N}, +\infty)$.

(Ⅳ) Defining the closure of $\overline{\Omega}_{1}$ as

$\overline{\Omega}_{1}$ is the invariant set on system (2.2).

Let

if

where $\cdot$ is the product of two scalar vectors, the vector field will enter in $\overline{\Omega}_{1}$ at the end, that $\overline{\Omega}_{1}$ is the invariant set. By calculating

For all $(w, g)\in\overline{\Omega}_{1}$, $\frac{dw(t)}{dt} > 0, \frac{dg(t)}{dt} < 0$ holds. So we get $\sigma(+\infty) = 0$ with $(e^{k_{1}D}ET, 0)\in N$ and $P(+\infty) = \frac{bET}{1+ET}$. We assume that there exists a positive $g_{*}$ that makes $\sigma(+\infty) = g_{*}$ and $P_{*} = (ET, g_{*})\in M$. Take another point $P_{1} = (ET, g_{1})$ and $0 < g_{1} < g_{*}$. By the invariance of set $\overline{\Omega}_{1}$ and uniqueness of solution, we get that the orbit initiating $P_{1}$ will arrive at point $P_{0}^{+} = (e^{-k_{1}D}ET, g_{0}^{+})\in N$ and $g_{0}^{+} > +\infty$, which is a contradiction. Therefore, we obtain $\sigma(+\infty) = 0$ and $P(+\infty) = \frac{bET}{1+ET}$, so $P$ has a horizontal asymptote $P = \frac{bET}{1+ET}$.

(Ⅴ) From the proof (Ⅰ), we get that $P$ is increasing on $[0, g_{N}]$ and decreasing on $[g_{N}, +\infty)$. And for $\forall \overline{n}\in(0, +\infty)$, $P(\overline{n})\leq P(g_{N})$ satisifty. Therefore $P$ takes the maximal value at $g_{N}$, which is also the maximum value. And the maximum $P(g_{N}) = e^{-k_{2}D}g(e^{-k_{1}D}ET, g_{N})+\frac{bET}{1+ET}$.

By proof (Ⅱ), we can easily get that $P$ takes the minimum value at $g = 0$, and the minimum value is $\frac{bET}{1+ET}$. This completes the proof.

4.   Periodic solution of system (2.2)

4.1.   Boundary periodic solution

For system (2.2), if sterile mosquitoes are no longer released, the system has boundary periodic solution. At this point, we get the following system:

The initial condition is $w(0^{+}) = e^{-k_{1}D}ET$, and solving the above equation we get

Suppose $w(t)$ reaches $L_{3}$ at time T, let $w(t) = ET$, we obtain:

Further solve the above equation on $T$, where $T$ represents the boundary periodic solution, we have

Therefore, the boundary periodic solution of model (2.2) with period $T$ is

Theorem 4.1. If the condition $\left|\chi_{1}\right| < 1$ holds, then the boundary periodic solution of system (2.2) is orbitally asymptotically stable, where

where $\Delta_{1} = \frac{e^{-k_{2}D}e^{-k_{1}D}(r-\varphi_{1}-a_{1}e^{-k_{1}D})}{r-\varphi_{1}-a_{1}ET},$ $I_{1} = \frac{r-\varphi_{1}-\varphi_{2}}{r-\varphi_{1}}ln\frac{a_{1}ET-e^{k_{1}D}r+e^{k_{1}D}\varphi_{1}}{a_{1}ET-r+\varphi_{1}},$ $I_{2} = (\frac{2a_{1}+a_{2}}{a_{1}})ln\left[1+\frac{r-\varphi_{1}-e^{k_{1}D}r+e^{k_{1}D}\varphi_{1}}{(a_{1}ET-r+\varphi_{1})(r-\varphi_{1})e^{-k_{1}D}}\right].$

Proof. Let $F(w, g) = w(\frac{rw}{w+g}-\varphi_{1}-a_{1}(w+g))$, $G(w, g) = -g(\varphi_{2}+a_{2}(w+g))$, $\alpha(w, g) = w(e^{-k_{1}D}-1)$, $\beta(w, g) = g(e^{-k_{1}D}-1)+\frac{bET}{1+ET}$, $\Phi(w, g) = w-ET, (w^{T}(T), g^{T}(T)) = (ET, 0), (w^{T}(T^{+}), g^{T}(T^{+})) = (e^{-k_{1}D}ET, 0).$

Then

and

where

Furthermore,

From article [45,Lemma 1], we know that if the condition $\left|\chi_{1}\right| < 1$ holds, the order-1 periodic solution of system (2.2) is orbitally asymptotically stable. This completes the proof.

4.2.   The existence and stability of periodic solution

From the Theorem 3.1, we get that system (2.2) must exist a fixed point of the Poincar$\acute{\mbox{e}}$ map, which represents system (2.2) exists an order-1 periodic solution. So we first talk about the stability of order-1 periodic solution $(\xi(t), \eta(t))$. Without loss of generality, we assume that the period of the order-1 periodic solution is $T$.

Theorem 4.2. The order-1 period solution $(\xi(t), \eta(t))$ is orbitally asymptotically stable if

where $Q(t) = \frac{rwg}{(w+g)^{2}}+\frac{rw}{w+g}-2a_{1}w-a_{1}g-\varphi_{1}-\varphi_{2}-a_{2}w-2a_{2}g.$

Proof. We define the starting and the ending point of the order-1 periodic solutions as $(\xi(T), \eta(T)) = (ET, 0)$ and $(\xi(T^{+}), \eta(T^{+})) = (e^{-k_{1}D}ET, e^{-k_{2}D}\eta_{0}+\frac{bET}{1+ET})$ respectively.

Thus,

and

The multiplier $\chi_{2}$ is obtained

From the (4.4), then $|\chi_{2}| < 1$ holds. By article [45,Lemma 1], we get that the order-1 periodic solution is orbitally asymptotically stable. This completes the proof.

Theorem 4.3. If $P(g_{N}) < g_{N}$, then system (2.2) has an globally asymptotically stable order-1 periodic solution.

Proof. If $P(g_{N}) < g_{N}$, by Theorem 3.1 available Poincar$\acute{\mbox{e}}$ map $P$ has a unique fixed point on $\overline g$, and $0 < \overline g\leq g_{n}$. This means that system (2.2) has a unique order-1 periodic solution that is orbitally asymptotically stable.

For any trajectory that starts at point $(e^{-k_{1}D}ET, g_{0}^{+})$, if $g_{0}^{+}\in[0, g_{N}]$, then $g_{0}^{+} < P(g_{0}^{+}) < \overline g$ from the Theorem 3.1. After $n$ times pulses, we get monotonically bounded sequence $P^{n}(g_{0}^{+})$, then $\lim\limits_{n\rightarrow+\infty}P^{n}(g_{0}^{+}) = \overline g$.

For any $g_{0}^{+} > g_{N}$, we have two situations as follows: (a) for all n, if $P^{n}(g_{0}^{+}) > \overline{g}$, then with the increases of $n$, $P^{n}(g_{0}^{+})$ monotonically decreases and $\lim\limits_{n\rightarrow+\infty}P^{n}(g_{0}^{+}) = \overline g$ from the Theorem 3.1. (b) If $P^{n}(g_{0}^{+}) > \overline{g}$ is not necessarily holds with all $n$, we make $n_{1}$ is the smallest positive integer and $P^{n_{1}}(g_{0}^{+}) < \overline{g}$. From the case (a), we get there must have a positive integer $n_{2}(n_{2} > n_{1})$, which makes $P^{n_{2}}(g_{0}^{+})$ monotonically increases with the increase of $n_{2}$, therefore $\lim\limits_{n_{2}\rightarrow+\infty}P^{n_{2}}(g_{0}^{+}) = \overline{g}$. Thus, the unique order-1 periodic solution is globally asymptotically stable. This completes the proof.

Theorem 4.4. When $P(g_{N}) > g_{N}, P^{2}(g_{N})\geq g_{N}$, the system (2.2) has a stable order-1 periodic solution or a stable order-2 periodic solution.

Proof. According to the Theorem 3.1, we get that $P$ increases on $[0, g_{N}]$. So when $P(g_{N}) > g_{N}$, $P$ has no fixed point on $[0, g_{N}]$. Then there is a positive integer j which satisfies $g_{j-1}^{+} < g_{N}, g_{j}^{+}\geq g_{N}$, by the definition of Poincar$\acute{\mbox{e}}$ map we obtain $g_{j}^{+} = P(g_{j-1}^{+})\leq P(g_{N})$ and $g_{j}^{+}\in [g_{N}, P(g_{N})]$.

For $g_{0}^{+}\in(g_{N}, +\infty)$, $P$ is monotonically decreasing on $[g_{N}, +\infty)$, so $g_{1}^{+} = P(g_{0}^{+})\leq P(g_{N})$. Then to any integer $j$ such that $g_{j}^{+}\in [g_{N}, P(g_{N})]$. Therefore we get $g_{j}^{+}\in [g_{N}, P(g_{N})]$ always hold.

From $P$ is monotonically decreasing on $[g_{N}, P(g_{N})]$, $P^{2}$ is increasing on $[g_{N}, P(g_{N})]$, thus

Based on the above conclusions, for any $g_{0}^{+}\in [g_{N}, P(g_{N})]$, we assume that the $g_{1}^{+} = P(g_{0}^{+})\neq g_{0}^{+}, g_{2}^{+} = P^{2}(g_{0}^{+})\neq g_{0}^{+}, g_{n}^{+} = P^{n}(g_{0}^{+})$. This means the solution of model (2.2) starts from $(e^{-k_{1}D}ET, g_{0}^{+})$, which is not the order-1 and the order-2 periodic solution. In other words, if $g_{1}^{+} = P(g_{0}^{+}) = g_{0}^{+}, g_{2}^{+} = P^{2}(g_{0}^{+}) = g_{0}^{+}$, then $P$ takes the fixed point at $g_{0}^{+}$, system (2.2) has both order-1 periodic solution and order-2 periodic solution. So we should discuss the parameters of $g_{N}, P(g_{N}), g_{0}^{+}, g_{1}^{+} and\ g_{2}^{+}$. Thus, we discuss the following four cases:

(Ⅰ) $P(g_{N})\geq g_{1}^{+} > g_{0}^{+} > g_{2}^{+}\geq g_{N}$. In this case (Fig. 3a), $g_{3}^{+} = P(g_{2}^{+}) > P(g_{0}^{+}) = g_{1}^{+}, g_{2}^{+} = P(g_{1}^{+}) < P(g_{3}^{+}) = g_{4}^{+}$. So $g_{3}^{+} > g_{1}^{+} > g_{0}^{+} > g_{2}^{+} > g_{4}^{+}$. By summing up the above relations, we obtain

(Ⅱ) $P(g_{N})\geq g_{1}^{+} > g_{2}^{+} > g_{0}^{+}\geq g_{N}$. In this case (Fig. 3b), we have $g_{1}^{+} = P(g_{0}^{+}) > P(g_{2}^{+}) = g_{3}^{+} > g_{2}^{+} = P(g_{1}^{+})$ and $P(g_{2}^{+}) = g_{3}^{+} > g_{4}^{+} = P(g_{3}^{+}) > P(g_{1}^{+}) = g_{2}^{+}$, therefore $g_{1}^{+} > g_{3}^{+} > g_{4}^{+} > g_{2}^{+} > g_{0}^{+}$. By induction, one obtains

(Ⅲ) $P(g_{N})\geq g_{0}^{+} > g_{2}^{+} > g_{1}^{+}\geq g_{N}$. Like case (2), we obtain

(Ⅳ) $P(g_{N})\geq g_{2}^{+} > g_{0}^{+} > g_{1}^{+}\geq g_{N}$. By using the similar process as case (1), we get

For case (Ⅱ) and (Ⅲ), we know that $P^{2n}(g_{0}^{+}) = g_{2n}$ is monotonically increasing and $P^{2n+1}(g_{0}^{+}) = g_{2n+1}^{+}$ is monotonically decreasing. So there exists a unique $\overline{g}$ that makes

Or there exists two fixed values $\overline{g}_{1}$, $\overline{g}_{2}\in[g_{N}, P(g_{N})]$ and $\overline{g}_{1}\neq\overline{g}_{2}$ such that

For case (Ⅰ) and (Ⅳ), only the second case is satisfied.

Thus, system (2.2) has an order-1 or order-2 periodic solution for case (Ⅱ) and case (Ⅲ), has an order-2 periodic solution for case (Ⅰ) and case (Ⅳ). This completes the proof.

Theorem 4.5. When $P(g_{N}) > g_{N}, P^{2}(g_{N}) < g_{min}^{+}$, where $g_{min}^{+} = min\{g^{+}, P(g^{+}) = g_{N}\}$. Then model (2.2) has a nontrivial order-3 periodic solution. And system (2.2) exist order-k $(k\geq3)$ periodic solution.

Proof. When $P(g_{N}) > g_{N}$, the Poincar$\acute{\mbox{e}}$ map have a unique fixed point $\overline{g}$ and $\overline{g}\in (g_{N}, P(g_{N}))$. In order to study the existence of order-3 periodic solution, we need to get a fixed point $\widetilde{g}\in [0, +\infty)$ such that $P^{3}(\widetilde{g}) = \widetilde{g}$ and $P(\widetilde{g})\neq \widetilde{g}$. By the Theorem 3.1 $P^{3}(g)$ is continuous on $[0, +\infty)$, therefore $P^{3}(0) = P^{2}(\frac{bET}{1+ET}) > 0$ and $P^{3}(g_{min}^{+}) = P^{2}(g_{N}) < g_{min}^{+}$.

Further, it follows from the intermediate value theorem and the continuity of $P^{3}$, system (2.2) must have a positive point $\widetilde{g}\in (0, g_{min}^{+})$ where $P^{3}(\widetilde{g}) = \widetilde{g}$ and $g_{min}^{+} < g_{N}$. We can easily get that $\widetilde{g}\neq g_{N}$. Therefore, system (2.2) has an order-3 periodic solution. From Sarkovskii^[46] theorem and [45,Definition 2] that the order-k periodic solution exists on system (2.2). This complete the proof.

5.   Numerical simulations

In the section 4, we obtain the boundary periodic solution of system (2.2) and the expressions of $T$ and $D$. Therefore, we can analyze the key factors $T$ and $ET$ that affect the dose D of insecticides. Let $r = 2.5, \varphi_{1} = 0.2, a_{1} = 0.3, \varphi_{2} = 0.01, a_{2} = 0.05$, we get when the threshold $ET$ increases, the insecticides dose $D$ decreases, and the downward trend becomes more and more obvious (Fig. 4a). Further, when the period T increases, the insecticides dose D also increases (Fig. 4b). Therefore, the effects of $ET$ and $T$ should be taken into account in the actual dose of insecticide sprayed.

We use the boundary periodic solution stability condition (4.3) to determine whether the wild mosquito population can be controlled by chemical control alone, that is, $\chi_{1} < 1$ means that the wild mosquito population can be controlled below $ET$ under the conditions of only chemical control. Therefore, we use numerical simulation to study the effect of threshold $ET$ on stability condition $\chi_{1}$ (Fig. 5), where the parameter values are $r = 2.5, \varphi_{1} = 0.2, a_{1} = 0.3, \varphi_{2} = 0.01, a_{2} = 0.05, k_{1} = 1.6, D = 1.6$. From the Fig. 5, we get $\chi_{1} < 1$ when $ET$ is relatively small, and $\chi_{1} > 1$ when $ET$ increases to a certain value. This indicates that when the value of $a$ is fixed, a smaller $ET$ value is more conducive to controlling the wild mosquito population. And the early use of chemical control makes the control of wild mosquitoes. But when the $ET$ get large, repeated chemical controls can make wild mosquito populations resistant and cannot be effectively controlled. So the use of integrated mosquito control is a reasonable and effective method.

Let parameters $r = 2.5, \varphi_{1} = 0.2, a_{1} = 0.3, \varphi_{2} = 0.01, a_{2} = 0.05, ET = 2.5$, we get that system (2.2) has order-1 periodic solution and it is globally asymptotically stable (Fig. 6a). Figure (6b) and (6c) are time series diagram of $w, g$. The results show that the wild mosquito population can be controlled below the threshold $ET$ when using integrated mosquito control.

6.   Conclusion

In this paper, we present the wild and sterile mosquito system with state-dependent feedback control and study its global dynamics. According to the introduction and main results of the article, system (2.2) using state-dependent feedback control is not only practical, but also rich in dynamic behavior.

In the article we define the Poincar$\acute{\mbox{e}}$ map $P$ and study its main properties such as monotonicity, continuity, extremum, and fixed point. We obtain an expression for the boundary period solution with the period $T$ and the pesticide dose $D$, which proves that it is stable under certain conditions. Based on the properties of Poincar$\acute{\mbox{e}}$ map, we prove the existence and stability of the order-1 periodic solution of system(2.2), and solve the dynamic complexity of the system, such as the proof of the existence of the order-3 periodic solution.

We use numerical simulations to verify the impact of key parameters and validate our conclusions. By analyzing the key influencing factors of pesticide $D$, we get that when the threshold $ET$ increases, $D$ decreases continuously, and when the period $T$ increases, $D$ increases continuously. Therefore, the dose of pesticide should consider the effects of $ET$ and $T$ at the same time. By studying the influence of threshold $ET$ on the stability condition $\chi_{1}$ of the boundary period, we get that the population of wild mosquitoes cannot be effectively controlled when $ET$ is large. This verifies the necessity of adopting an integrated mosquito control. Through numerical simulation, we verify the existence and stability of the order-1 periodic solution. That is, through integrated mosquito control, wild mosquitoes can be controlled below the threshold.

Compared with the previous state-dependent feedback control model, we summarize some of the highlights of this paper: (1) We study the global dynamics of the model through the properties of the Poincar$\acute{\mbox{e}}$ map, and the existence of periodic solution is proved by studying the fixed point of Poincar$\acute{\mbox{e}}$ map. (2) The biological control method used in this paper is sterile insect technique. And wild mosquito populations can be more effectively controlled by using a new proportional release method when releasing sterile mosquitoes. (3) We study the effects of insecticide dose on integrated mosquito control and the influence of key parameters is analyzed by numerical simulation. In future research work, we will add optimization problems to reduce control costs.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (11371230), the SDUST Research Fund (2014TDJH102), Shandong Provincial Natural Science Foundation, China (ZR2015AQ001), Joint Innovative Center for Safe And Effective Mining Technology and Equipment of Coal Resources, SDUST Innovation Fund for Graduate Students (SDKDYC170225), SDUST Innovation Fund for Graduate Students(SDKDYC190351).

Conflict of interest

The authors declare that there is no conflict of interests regarding the publication of this paper.