Global Dynamics of a Predator–Prey Model with Fear Effect and Impulsive State Feedback Control

Abstract

In this paper, a predator–prey model with fear effect and impulsive state control is proposed and analyzed. By constructing an appropriate Poincaré map, the dynamic properties of the system, including the existence, nonexistence, and stability of periodic solutions are studied. More specifically, based on the biological meaning, the pulse and the phase set are firstly defined in different regions as well as the corresponding Poincaré map. Subsequently, the properties of the Poincaré map are analyzed, and the existence of a periodic solution for the system is investigated according to the properties of the Poincaré map. We found that the existence of the periodic solution for the system completely depends on the property of the Poincaré map. Finally, several examples containing numerical simulations verify the obtained theoretical result.

1. Introduction

As an important branch of ecology, population ecology takes populations as the research object to study the changes and regulation mechanisms affecting population size or quantity in time and space. The complexity of ecological relationships in nature has increasingly required the application of mathematical methods to ecology, which has gradually made mathematical ecology one of the most mature branches of biological mathematics so far. The development of mathematical ecology can be traced back to the study of predator and prey models by the ecologist Lotka [1] and the mathematician Voltra [2] in the first half of the 20th century. Generally, there are four basic types of interactions between the two groups: competing type, mutually beneficial type, parasitic type, predator and prey type. Among them, the predator–prey model has received a great deal of attention and research interest [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]. The predator–prey interaction is one of the basic relationships between biological populations, and it is a popular research topic in ecology and biomathematics. A simple model describing the relationship between predator and prey can be written as follows:

$$\left\{\begin{array}{c}\frac{dx}{dt}=x(r_1-a_{11}x+a_{12}y),\hfill \\ \frac{dy}{dt}=y(r_2-a_{22}y+a_{21}x).\hfill \end{array}\right.$$

(1)

However, a new view is that the mere presence of a predator may change the behavior and physiology of the prey such that it can have a more powerful impact on the prey population than in direct predation [20,21,22,23]. Animals respond to perceived predation risks and exhibit various anti-predator responses. For example, when the prey is assessing the risk of predation, they may choose to abandon their original high-risk habitats and move to low-risk habitats, which may lead to energy consumption, especially when living conditions in the low-risk habitats are poor [21,24,25,26,27]. These changes also include changes in foraging behavior, vigilance, and physiology [21,24,25,26,27,28].

Similarly, fearful prey will reduce their outings for food, and the reduction in food may reduce birth and survival rates. The high risk of acute predation can even cause the prey to temporarily leave the habitat or foraging site [21,22]. Because of the complexity of the interaction between prey and predator, it should not be simply described in terms of direct predation, and the influence of fear should be considered [21,24,25,26,27]. As for the extent to which fear affects prey populations, biologists have performed a lot of experiments to show the impact of fear on population size. Recently, Zanette et al. [29] conducted an experiment under the assumption that any impact on reproduction can only be attributed to fear. They broadcast the voice of predators to some songbird populations and found that the fear of predators alone reduced the number of offspring that song sparrows and sparrows could produce by 40%. Biologists attribute this phenomenon to the fear effect, that is, the risk of predation has an impact on the birth rate and survival rate of offspring. It is manifested in the low egg-laying rate of female song sparrows, and the low egg-hatching rate and low survival rate of nestlings. This fear effect has also been found in studies of other species such as other birds [30], elk [31], snow rabbits [32], dugongs [33], and aphids [34,35].

Recently, Wang et al. [36] constructed a predator–prey model with the fear effect, which is written as follows:

$$\left\{\begin{array}{c}\frac{dx}{dt}=r_{0}xf(k,y)-dx-a{x}^2-g\left(x\right)y,\hfill \\ \frac{dy}{dt}=y(-d_1+cg\left(x\right)),\hfill \end{array}\right.$$

(2)

where $r_0$ represents the self-growth rate of the prey, d and $d_1$ are the death rate of preys and predators, respectively. a indicates the mortality rate of the prey due to interspecific competition, $g\left(x\right)$ represents the predation function of the predator to the prey, c represents the utilization of energy by predators after ingesting prey, k represents the degree of fear that the bait resists in the face of predation by the predator. It is assumed that the fear effect has an impact on the birth rate of the prey, the term $r_{0}xf(k,y)$ involves a factor $f(k,y)$, representing the cost of anti-predator defense due to fear, while the function $f(k,y)$ satisfies the following equation:

$$\left\{\begin{array}{c}f(0,y)=1,f(k,0)=1,\underset{k\to \infty }{lim}f(k,y)=0,\underset{y\to \infty }{lim}f(k,y)=0,\hfill \\ \frac{\partial f(k,y)}{\partial k}<0,\frac{\partial f(k,y)}{\partial y}<0.\hfill \end{array}\right.$$

Obviously, the function $f(k,y)=\frac{1}{1+ky}$ meets the conditions above. The authors discussed the dynamics of the model with the functional responses $px$ and $\frac{x}{1+qx},$ respectively. Based on model (2), Zhu et al. [37] proposed a predator–prey model incorporating the fear effect, with the prey and the predator having other food resources:

$$\left\{\begin{array}{c}\frac{dx\left(t\right)}{dt}=\frac{r_{0}x\left(t\right)}{1+ky\left(t\right)}-dx\left(t\right)-a{x}^2\left(t\right)-bx\left(t\right)y\left(t\right),\hfill \\ \frac{dy\left(t\right)}{dt}=cbx\left(t\right)y\left(t\right)+my\left(t\right)-h{y}^2\left(t\right),\hfill \end{array}\right.$$

(3)

where m is the utilization of nutrients obtained from other prey for the predator’s reproduction, and h is the mortality rate due to environmental density constraints.

State feedback control is the control of the system according to the state of system variables. Because rapid changes in system variables can be described by pulses, it is also called impulsive state feedback control. In the past two decades, impulsive state feedback control has received extensive attention which has been widely used in integrated pest management [38,39,40,41,42,43,44], microbial culture [45,46,47,48,49,50,51], disease control, and so on [52,53,54,55,56,57,58,59,60]. Among them, Tang et al. [38] established an impulsive state feedback control system based on a predator–prey system to model the implementation process of integrated pest management:

$$\left\{\begin{array}{c}\left.\begin{array}{c}\frac{dx\left(t\right)}{dt}=x\left(t\right)(a-by\left(t\right)),\hfill \\ \frac{dy\left(t\right)}{dt}=y\left(t\right)(cx\left(t\right)-d),\hfill \end{array}\right\}x\ne ET,\hfill \\ \left.\begin{array}{c}\Delta x\left(t\right)=-px\left(t\right),\hfill \\ \Delta y\left(t\right)=\tau ,\hfill \end{array}\right\}x=ET,\hfill \end{array}\right.$$

(4)

where x represents the number of pest populations, y represents the number of natural enemy populations, and $ET$ is the economic threshold, i.e., the number of pests when control measures must be taken to prevent the economic injury level ($EIL$) from being reached and exceeded, and where the $EIL$ is the lowest pest population density that will cause economic damage. It is assumed that integrated pest management is carried out, and measures such as physical killing, chemical treatment, and biological control are adopted to control the number of the pests when the number of pests reaches the economic threshold ($ET$).

In this paper, we consider an example of integrated pest management. The bird cherry oat aphid is a small grain aphid common in wheat fields that attacks wheat, oats, rye, and barley. Its common natural enemy is the seven-spotted lady beetle (Coccinella septempuctata). Research has shown that predator odor has an important effect on the growth of the bird cherry oat aphid, and aphids can sense the information left in the surroundings by the seven-spotted lady beetle through hearing, sight, and smell [34]. Moreover, aphids have an innate fear of natural enemies. When the presence of natural enemies is detected, aphids will take further defensive measures [35]. Thus, based on model (3) and model (4), we establish the following system:

$$\left\{\begin{array}{c}\left.\begin{array}{c}\frac{dx\left(t\right)}{dt}=\frac{r_{0}x\left(t\right)}{1+ky\left(t\right)}-dx\left(t\right)-a{x}^2\left(t\right)-bx\left(t\right)y\left(t\right),\hfill \\ \frac{dy\left(t\right)}{dt}=cbx\left(t\right)y\left(t\right)+my\left(t\right)-h{y}^2\left(t\right),\hfill \end{array}\right\}x\ne ET,\hfill \\ \left.\begin{array}{c}\Delta x\left(t\right)=-px\left(t\right),\hfill \\ \Delta y\left(t\right)=-qy\left(t\right)+\tau .\hfill \end{array}\right\}x=ET,\hfill \end{array}\right.$$

(5)

where x represents the number of bird cherry oat aphids, and y represents the number of the seven-spotted lady beetles. $\frac{r_{0}x\left(t\right)}{1+ky\left(t\right)}$ contains the effect of fear of natural enemies. When pest numbers are below ET, we rely on natural enemies in nature to control pest populations. When the pest population reaches ET, integrated pest management will be implemented to suppress the pest population. These include spraying pesticides and releasing predators. That is, the implementation of the control strategy is decided according to the number of pests, which can obviously be represented by a state feedback control system.

Our purpose is to explore the dynamics of the system under impulse state feedback control. The structure of the paper is as follows. Some basic results for model (3) is summarized in Section 2. The dynamics of the system with impulsive state feedback control is discussed using the Poincaré map in Section 3. In Section 4, we give some examples with numerical simulations to verify the obtained theoretical result. In Section 5, we summarize the full paper.

2. Some Basic Results for Model (3)

Without the impulse effect, Zhu et al. studied the dynamic behavior of model (3) (see Theorems 2.1, 2.2, and 3.1 in [37]). Denote $R=\frac{r_0{h}^2}{(km+h)(dh+mb)};$ below, we have summarized some basic results in Lemma 1 and

Table 1. Linear analysis of equilibrium points.

Points	$\mathit{R}>1$	$\mathit{R}<1$ and ${\mathit{r}}_0>\mathit{d}$	${\mathit{r}}_0<\mathit{d}$
$E_0(0,0)$	unstable node	unstable node	saddle
$E_1(0,\frac{m}{h})$	saddle	stable node	stable node
$E_2(\frac{r_0-d}{a},0)$	saddle	saddle	does not exist
$E_3(x^{\ast },y^{\ast })$	stable node or focus	does not exist	does not exist

.

Lemma 1.

The system equilibrium point has the following characteristics:

(a) 

$E_0(0,0)$ is always unstable, and it is an unstable node or a saddle point;

(b) 

If $R>1,$$E_3(x^{\ast },y^{\ast })$ is a locally asymptotically stable positive equilibrium, where $y^{\ast }=\frac{m+cb{x}^{\ast }}{h}$ and $x^{\ast }$ is the unique positive solution of the equation $A_1{x}^2+A_{2}x+A_3=0$ where $A_1=c^{2}k{b}^3+achkb,$$A_2=cdhkb+2ckm{b}^2+ahkm+ch{b}^2+a{h}^2,$$A_3=dhkm+k{m}^2+dh-h^2{r}_0+hmb;$

(c) 

If $R<1$ and $r_0>d,$ the system equilibrium point $E_1(0,\frac{m}{h})$ is locally asymptotically stable;

(d) 

If $r_0<d,$$E_0(0,0)$ is a saddle point, point $E_1(0,\frac{m}{h})$ is locally asymptotically stable, and point $E_2(\frac{r_0-d}{a},0)$ does not exist.

3. Dynamic Analysis of the Model with Impulsive State Feedback Control

3.1. The Impulsive Set and the Phase Set

For system (3), there exist two X-isoclines—$L_1:x=\frac{r_0}{a(1+ky)}-\frac{d}{a}-\frac{by}{a}$ and the X-axis, and two Y-isoclines—$L_2:x=\frac{hy}{cb}-\frac{m}{cb}$ and the Y-axis.

Denote the straight lines with $L_3:x=(1-p)ET$ and $L_4:x=ET.$ The intersection of the straight line $L_4$ and the isoclinic line $L_1$ is denoted by $Q(ET,y_{ET}),$ where:

$$y_{ET}=\frac{-b-A+\sqrt{{(-b+A)}^2+4kb{r}_0}}{2kb},A=kd+kaET.$$

The intersection of the straight line $L_3$ and the isoclinic line $L_1$ is denoted by $P((1-p)ET,y_{p_{ET}}),$ where:

$$y_{p_{ET}}=\frac{-b-A_1+\sqrt{{(-b+A_1)}^2+4kb{r}_0}}{2kb},A_1=kd+ka(1-p)ET.$$

The intersection of the straight line $L_3$ and the isoclinic line $L_2$ is denoted by $R((1-p)ET,y_{r_{ET}}),$ where:

$$y_{r_{ET}}=\frac{cb(1-p)ET+m}{h}.$$

Next, when the system (3) has a stable focus (see Figure 1), we define the impulsive set $I_s$ and the corresponding phase set $P_s$ of the system in two cases.

Case I: $(1-p)ET<ET<x^{\ast },$ that is, both the impulsive set and the phase set are on the left of equilibrium $E^{\ast }.$ Then, the impulsive set is defined as $I_s=\{(x,y)\in R^2|x=ET,0\le y\le y_{ET}\},$ and the corresponding phase set is defined as $P_s=\{(x^{+},y^{+})\in R^2|x^{+}=(1-p)ET,y^{+}\in D_0\},$ where $D_0=[\tau ,(1-q)y_{ET}+\tau ]$ (see Figure 2a).

Case II: $(1-p)ET<x^{\ast }<ET<\frac{r_0-d}{a},$ that is, the impulsive set and the phase set are on both sides of equilibrium $E^{\ast }$. Obviously, there is a trajectory $\Gamma$ and a straight line $L_4$ tangent to point $Q(ET,y_{ET}),$ the intersection of trajectory $\Gamma$ and isoclines; $L_1$ is denoted by $P_2(x_{p2},y_{p2})$. According to the relative position of point $P_2(x_{p2},y_{p2})$ and straight line $L_3,$ we divide the discussion into two subcases.

Case IIa: If $(1-p)ET<x_{p2},$ the trajectory ${\Gamma }_1$ starting from point $P(ET,y_{p_{ET}})$ will intersect with the line $L_4:x=ET$ at the point $P_1(ET,y_{p_1}).$ At the same time, we found that the point $P((1-p)ET,y_{p_{ET}})$ is the tangent point between the trajectory ${\Gamma }_1$ and the line $L_3:x=(1-p)ET.$ Hence, we have $I_s=\{(x,y)\in R^2|x=ET,0\le y\le y_{p_1}\}$ and $P_s=\{(x^{+},y^{+})\in R^2|x^{+}=(1-p)ET,y^{+}\in D_1\},$ where $D_1=[\tau ,(1-q)y_{p_1}+\tau ]$ (see Figure 2b).

Case IIb: If $x_{p_2}<(1-p)ET,$ the trajectory $\Gamma$ and the straight line $L_3:x=(1-p)ET$ intersect at two points, namely $T_1((1-p)ET,y_{T_1})$ and $T_2((1-q)ET,y_{T_2}),$ respectively, where $y_{T_1}<y_{T_2}.$ Thus, $I_s=\{(x,y)\in R^2|x=ET,0\le y\le y_{ET}\},$ and we have $P_s=\{(x^{+},y^{+})\in R^2|x^{+}=(1-p)ET,y^{+}\in D_2\},$ where $D_2=\{[0,y_{T_1}]\cup [y_{T_2},+\infty )\}$ (see Figure 2c).

Figure 2. Impulsive sets and phase sets of different positions. (a) $(1-p)ET<ET<x^{\ast }$; (b) $(1-p)ET<x_{p_2}<x^{\ast }<ET$; (c) $x_{p_2}<(1-p)ET<x^{\ast }<ET$.

Figure 2. Impulsive sets and phase sets of different positions. (a) $(1-p)ET<ET<x^{\ast }$; (b) $(1-p)ET<x_{p_2}<x^{\ast }<ET$; (c) $x_{p_2}<(1-p)ET<x^{\ast }<ET$.

3.2. Dynamic Analysis of System (5) for Case I

3.2.1. Definition and Properties of Poincaré Mapping

To establish the Poincaré map of system (5), let us denote two sections as follows: $S_{ET}=\left\{\right(x,y\left)\right|x=ET,y\ge 0\},$$S_{p_{ET}}=\left\{\right(x,y\left)\right|x=(1-p)ET,y\ge 0\}.$ We choose section $S_{p_{ET}}$ as the Poincaré section and pick the point $P_k^{+}=((1-p)ET,y_k^{+})$ in the Poincaré section $S_{p_{ET}}.$ Then, any trajectory $\psi (t,t_0,(1-p)ET,y_k^{+})=(x(t,t_0,(1-p)ET,y_k^{+}),y(t,t_0,(1-p)ET,y_k^{+}))$ starting from point $P_k^{+}$ will intersect with the section $S_{ET}$ at point $P_{k+1}=(ET,y_{k+1})$. Assuming that the elapsed time is $t_1,$ a finite time, then we have $y_{k+1}=y(t_1,t_0,(1-p)ET,y_k^{+})=\mathrm{\Phi }\left(y_k^{+}\right).$

For the sake of brevity, let us denote $y((1-p)ET,y_k^{+})=y(t_1,t_0,(1-p)ET,y_k^{+})$ throughout this paper. When the trajectory $\psi (t,t_0,(1-p)ET,y_k^{+})$ of system (5) intersects with the line $x=ET,$$\psi (t,t_0,(1-p)ET,y_k^{+})$ will undergo a pulse process and jump to point $P_{k+1}^{+}=((1-p)ET,y_{k+1}^{+})$ with $y_{k+1}^{+}=(1-q)y_{k+1}+\tau$ on $S_{p_{ET}}.$ We need to establish the mapping relationship between $P_k^{+}$ and $P_{k+1}^{+}.$ This is accomplished in two steps. For the first step, let us denote the following:

$$\left\{\begin{array}{c}\frac{r_{0}x\left(t\right)}{1+ky\left(t\right)}-dx\left(t\right)-a{x}^2\left(t\right)-bx\left(t\right)y\left(t\right)=P(x\left(t\right),y\left(t\right)),\hfill \\ cbx\left(t\right)y\left(t\right)+my\left(t\right)-h{y}^2\left(t\right)=Q(x\left(t\right),y\left(t\right)).\hfill \end{array}\right.$$

(6)

Consider the following scalar differential equation in a phase space:

$$\left\{\begin{array}{c}\frac{dy}{dx}=\frac{cbx\left(t\right)y\left(t\right)+my\left(t\right)-h{y}^2\left(t\right)}{\frac{r_{0}x\left(t\right)}{1+ky\left(t\right)}-dx\left(t\right)-a{x}^2\left(t\right)-bx\left(t\right)y\left(t\right)}=g(x,y),\hfill \\ y\left((1-q)ET\right)=y_0^{+}.\hfill \end{array}\right.$$

(7)

For a given region $\mathrm{\Omega }=\left\{(x,y)\right|x>0,y>0,x<\frac{r_0}{a(1+ky)}-\frac{d}{a}-\frac{by}{a}\},$ obviously, $g(x,y)$ is continuously differentiable on $\mathrm{\Omega }.$ Let $x_0^{+}=(1-p)ET,y_0^{+}=S,S\in I_s,$ and $S<y_{p_{ET}},$ then $(x_0^{+},y_0^{+})\in \mathrm{\Omega }.$ Therefore, we have the following:

$$\begin{array}{cc}\hfill y\left(x\right)=& y(x;x_0^{+},y_0^{+})\hfill \\ \hfill =& y(x;(1-p)ET,S)\hfill \\ \hfill =& y(x;S)\hfill \\ \hfill =& S+{\int }_{(1-p)ET}^{x}g(s,y(s,S))ds,(1-p)ET\le x\le ET.\hfill \end{array}$$

(8)

For the second step, consdering the impulsive effect, we have $y_{k+1}^{+}=(1-q)y_{k+1}+\tau .$ Therefore, the Poincaré map $P_M$ on $\mathrm{\Omega }$ can be expressed as follows:

$$P_M\left(S\right)=(1-q)y(ET,S)+\tau .$$

The Poincaré map $P_M$ of model (5) satisfies the following properties.

Theorem 1.

The Poincaré map $P_M$ has the following properties:

(I) 

The domain and range of the $P_M$ are $[0,+\infty )$ and $[\tau ,(1-q)y((1-p)ET,y_{p_{ET}})+\tau ]$, respectively. They increase on $[0,y_{p_{ET}}]$ and decrease on $(y_{p_{ET}},+\infty );$

(II) 

$P_M$ is continuously differentiable;

(III) 

$P_M$ is concave on $[y_{r_{ET}},y_{p_{ET}}];$

(IV) 

$P_M$ has a unique positive fixed point $y_f.$ If $\tau >0$, and $P_M\left(y_{p_{ET}}\right)<y_{p_{ET}},$ then $y_f\in (0,y_{p_{ET}});$ if $\tau >0$, and $P_M\left(y_{p_{ET}}\right)>y_{p_{ET}},$ then $y_f\in (y_{p_{ET}},+\infty );$

(V) 

For all $y\in [0,+\infty ),$$P_M\left(y\right)$ has a maximum value of $P_{M_{max}}=P_M\left(y_{p_{ET}}\right).$ $P_M$ has a minimum value of $P_{M_{min}}=P_M\left(0\right)=\tau .$

Proof. 

(I)

According to the vector field of system (3), it is easy to see that the definition domain of $P_M$ is $[0,+\infty ).$ For any $y_i^{+},y_j^{+}\in [0,y_{p_{ET}}]$ with $y_i^{+}<y_j^{+},$ by the uniqueness of the solution of model (3), it can be easily obtained that $y((1-p)ET,y_i^{+})<y((1-p)ET,y_i^{+});$ therefore, we have $P_M\left(y_i^{+}\right)<P_M\left(y_j^{+}\right).$ For any $y_i^{+},y_j^{+}\in (y_{p_{ET}},+\infty )$ with $y_i^{+}<y_j^{+},$ according to the the vector field of system (3), the orbit $\varphi (t,t_0,(1-p)ET,y_i^{+})$ starting from $y_i^{+}$ will first intersect with line $L_3$ and then line $L_4,$ and the intersection with $L_3$ is denoted by $y_{i^{\prime }}^{+}.$ The case of the orbits $\varphi (t,t_0,(1-p)ET,y_j^{+})$ and $\varphi (t,t_0,(1-p)ET,y_i^{+})$ is similar, and the intersection with $L_3$ is denoted by $y_{j^{\prime }}^{+}.$ By the uniqueness of the solution for model (3), we have $y_{i^{\prime }}^{+}>y_{j^{\prime }}^{+}.$ Then, similarly to the previous case, it can be discussed that $P_M\left(y_i^{+}\right)=P_M\left(y_{i^{\prime }}^{+}\right)<P_M(y_j^{+})=P_M\left(y_{j^{\prime }}^{+}\right).$ Therefore, $P_M$ is increasing on $[0,y_{p_{ET}}]$ and decreasing on $(y_{p_{ET}},+\infty ).$ Meanwhile, the range of $P_M$ is $[\tau ,(1-q)y((1-p)ET,y_{p_{ET}})+\tau ].$

(II)

Since both functions $P(x,y)$ and $Q(x,y)$ are continuous and differentiable, then based on the theorem of Cauchy and Lipschitz along with the parameters we have, $P_M$ is continuous and differentiable.

(III)

From (6), we can deduce that:

$$\begin{array}{cc}\hfill \frac{\partial g}{\partial y}=& \frac{cbx+m-2hy}{\frac{r_{0}x}{1+ky}-dx-a{x}^2-bxy}+\frac{(bx+\frac{r_{0}kx}{{(1+ky\left(t\right))}^2})(cbxy+my-h{y}^2)}{{(\frac{r_{0}x}{1+ky}-dx-a{x}^2-bxy)}^2}\hfill \end{array}$$

(9)

and

$$\begin{array}{cc}\hfill \frac{{\partial }^{2}g}{\partial y^2}=& \frac{-2h}{\frac{r_{0}x}{1+ky}-dx-a{x}^2-bxy}+\frac{(\frac{r_{0}kx}{{(1+ky)}^2}+bx)(cbx+m-2hy)}{{(\frac{r_{0}x}{1+ky}-dx-a{x}^2-bxy)}^2}\hfill \\ \hfill & +\frac{-2\frac{r_0{k}^{2}x}{{(1+ky)}^3}(cbxy+my-h{y}^2)+(bx+\frac{r_{0}kx}{{(1+ky)}^2})(cbx+m-2hy)}{{(\frac{r_{0}x}{1+ky}-dx-a{x}^2-bxy)}^2}\hfill \\ \hfill & +2\frac{{(\frac{r_{0}kx}{{(1+ky)}^2}+bx)}^2(cbxy+my-h{y}^2)}{{(\frac{r_{0}x}{1+ky}-dx-a{x}^2-bxy)}^3}\hfill \\ \hfill =& \frac{-2h}{\frac{r_{0}x}{1+ky}-dx-a{x}^2-bxy}+2\frac{bx(cbx+m-2hy)}{{(\frac{r_{0}x}{1+ky}-dx-a{x}^2-bxy)}^2}\hfill \\ \hfill & +2\frac{\frac{r_{0}kx}{{(1+ky)}^3}(\frac{cbx+m-hy}{1+ky}-hy)}{{(\frac{r_{0}x}{1+ky}-dx-a{x}^2-bxy)}^2}+2\frac{{(\frac{r_{0}kx}{{(1+ky)}^2}+bx)}^2(cbxy+my-h{y}^2)}{{(\frac{r_{0}x}{1+ky}-dx-a{x}^2-bxy)}^3}.\hfill \end{array}$$

(10)

Since $y_{r_{ET}}<y<y_{p_{ET}},$ we have $\frac{r_{0}x}{1+ky}-dx-a{x}^2-bxy>0$ and $cbx+m-hy<0,$ then we obtain $cbx+m-2hy<0.$ Therefore, we have $\frac{\partial g}{\partial y}<0,$$\frac{{\partial }^{2}g}{\partial y^2}<0$ for $y_{r_{ET}}<y<y_{p_{ET}}.$ On the other hand, we have the following equation:

$$\begin{array}{cc}\hfill \frac{\partial y(x,S)}{\partial S}=& e^{{\int }_{(1-p)ET}^x\frac{\partial }{\partial y}\left(\frac{Q(z,y(z,S\left)\right)}{P(z,y(z,S\left)\right)}\right)dz}<0.\hfill \\ \hfill \frac{{\partial }^{2}y(x,S)}{\partial S^2}=& \frac{\partial y(x,S)}{\partial S}{\int }_{(1-p)ET}^x\frac{{\partial }^2}{\partial y^2}\left(\frac{Q(z,y(z,S\left)\right)}{P(z,y(z,S\left)\right)}\right)\frac{\partial y(z,S)}{\partial S}dz<0.\hfill \end{array}$$

(11)

Hence, $P_M$ is monotonically increasing and concave for $y_{r_{ET}}<y<y_{p_{ET}}.$

(IV) Since $P_M$ is decreasing on $(y_{p_{ET}},+\infty ),$ then there exists a $\overline{y}\in (y_{p_{ET}},+\infty )$ such that $P_M\left(\overline{y}\right)<\overline{y}.$ In addition, $P_M\left(0\right)>0$ for $\tau >0;$ therefore, there exists $y_f\in [0,\overline{y})$ such that $P_M\left(y_f\right)=y_f,$ i.e., there exists a fixed point on $[0,+\infty )$ for $P_M.$

(a)

If $\tau >0,P_M\left(y_{p_{ET}}\right)<y_{p_{ET}},$ then the fixed point $y_f\in (0,y_{p_{ET}}).$ Since $P_M$ is decreasing on $[y_{p_{ET}},+\infty ),$ we have $P_M\left(y_k^{+}\right)<P_M\left(y_{p_{ET}}\right)<y_{p_{ET}}$ for $y_k^{+}\in [y_{p_{ET}},+\infty ),$ which indicates that no fixed point exists for $P_M$ on $[y_{p_{ET}},+\infty ).$ Next, we prove that the fixed point is unique.

On the one hand, if $P_M\left(y_{p_{ET}}\right)<y_{p_{ET}}$ and $P_M\left(y_{r_{ET}}\right)>y_{r_{ET}},$ since $P_M$ is concave on $[y_{r_{ET}},y_{p_{ET}}],$ a unique fixed point thus exists for $P_M$ on $[y_{r_{ET}},y_{p_{ET}}].$

On the other hand, if $P_M\left(y_{p_{ET}}\right)<y_{p_{ET}}$ and $P_M\left(y_{r_{ET}}\right)<y_{r_{ET}},$ there is at least one fixed point on $[0,y_{r_{ET}}].$ If $P_M$ has two fixed points with $\overline{y_1},\overline{y_2}\in (y_{n1}^{+},y_{n2}^{+}),$ then $P_M\left(\overline{y_1}\right)=\overline{y_1}$ and $P_M\left(\overline{y_2}\right)=\overline{y_2}.$ For any point $y_0^{+}\in (y_{n1}^{+},\overline{y_1}),$ from the successor function, we have $y_0^{+}<P_M\left(y_0^{+}\right).$ Since $P_M\left(\overline{y_1}\right)=\overline{y_1},$ there is $P_M\left(y_0^{+}\right)<y_0^{+}$ with $y_0^{+}\in (\overline{y_1},\overline{y_2}).$ Similarly, according to $P_M\left(\overline{y_2}\right)=\overline{y_2}$ and $P_M\left(y_{r_{ET}}\right)<y_{r_{ET}},$ then for any point $y_0^{+}\in (\overline{y_2},y_{n2}^{+}),$ we have $P_M\left(y_0^{+}\right)<y_0^{+},$ which contradicts $P_M\left(y_0^{+}\right)<y_0^{+}$ at any point of $y_0^{+}\in (\overline{y_1},y_{n2}^{+}).$ Thus, $P_M$ has a unique fixed point $y_f$ on $(0,y_{r_{ET}}).$

(b)

If $\tau >0,P_M\left(y_{p_{ET}}\right)>y_{p_{ET}},$ similarly, using the same method in case (a), we can prove that system (5) has only one fixed point in the interval $[y_{p_{ET}},+\infty ).$

(V)

According to the monotonicity of the function $P_M$ and the vector field of system (3), we can easily get that for all $y\in [0,+\infty ),$$P_M\left(y\right)$ has a maximum value of $P_{M_{max}}=P_M\left(y_{p_{ET}}\right).$ $P_M$ has a minimum value of $P_{M_{min}}=P_M\left(0\right)=\tau .$ □

3.2.2. Existence and Stability of the Boundary Order-1 Periodic Solution for $\tau =0$

If $y=0$ and $\tau =0,$ system (5) has a boundary order-1 periodic solution. Then, we discuss the the following subsystem:

$$\left\{\begin{array}{cc}\begin{array}{c}\frac{dx\left(t\right)}{dt}=(r_0-d)x\left(t\right)-a{x}^2\left(t\right),\hfill \end{array}\hfill & x\ne ET,\hfill \\ \begin{array}{c}x\left(t^{+}\right)=(1-p)x\left(t\right),\hfill \end{array}\hfill & x=ET.\hfill \end{array}\right.$$

(12)

Take the initial condition $x\left(0^{+}\right)=(1-p)ET$ into (12), then we get the following equation:

$$\begin{array}{c}\hfill x\left(t\right)=\frac{r_0-d}{a}\frac{1}{1+(\frac{r_0-d}{a(1-p)ET}-1)e^{-(r_0-d)t}}.\end{array}$$

(13)

Suppose that the trajectory reaches line $L_4$ at time $T,$ then we obtain the following:

$$\begin{array}{c}\hfill ET=\frac{r_0-d}{a}\frac{1}{1+(\frac{r_0-d}{a(1-p)ET}-1)e^{-(r_0-d)T}},\end{array}$$

(14)

then, we get the expression of T as follows:

$$\begin{array}{c}\hfill T=\frac{1}{r_0-d}ln\frac{r_0-d-a(1-p)ET}{(1-p)(r_0-d-aET)}.\end{array}$$

(15)

Therefore, the periodic solution of $(x^T\left(t\right),y^T\left(t\right))$ for system (5) is as follows:

$$\left\{\begin{array}{c}{x}^T\left(t\right)=\frac{r_0-d}{a}\frac{1}{1+(\frac{r_0-d}{a(1-p)ET}-1)e^{-(r_0-d)t}},\hfill \\ y^T\left(t\right)=0.\hfill \end{array}\right.$$

(16)

With regard to the stability of the boundary order-1 periodic solution, we have the following theorem.

Theorem 2.

If $\tau =0,$ ${\mu }_2<1,$ then the boundary order-1 periodic solution $(x^T\left(t\right),0)$ of system (5) is stable. Furthermore, if $\tau =0,$${\mu }_2>1,$ the boundary order-1 periodic solution $(x^T\left(t\right),0)$ is unstable.

Proof. 

Let us denote the following:

$$\begin{array}{cc}\hfill P(x,y)=& \frac{r_{0}x\left(t\right)}{1+ky\left(t\right)}-dx\left(t\right)-a{x}^2\left(t\right)-bx\left(t\right)y\left(t\right),\hfill \\ \hfill Q(x,y)=& cbx\left(t\right)y\left(t\right)+my\left(t\right)-h{y}^2\left(t\right),\hfill \\ \hfill \alpha (x,y)=& -px,\hfill \\ \hfill \beta (x,y)=& -qy+\tau ,\hfill \\ \hfill \varphi (x,y)=& x-ET,\hfill \end{array}$$

(17)

with

$$(x^T\left(T\right),y^T\left(T\right))=(ET,0),(x^T\left(T^{+}\right),y^T\left(T^{+}\right))=((1-p)ET,0).$$

By calculating the partial derivatives with respect to $x,y$, we get the following equation:

$$\begin{array}{cc}\hfill \frac{\partial P}{\partial x}=& \frac{r_0}{1+ky\left(t\right)}-d-2ax\left(t\right)-by\left(t\right),\hfill \\ \hfill \frac{\partial Q}{\partial y}=& cbx\left(t\right)+m-2hy\left(t\right),\hfill \\ \hfill \frac{\partial \alpha }{\partial x}=& -p,\hfill \\ \hfill \frac{\partial \beta }{\partial y}=& -q,\hfill \\ \hfill \frac{\partial \varphi }{\partial x}=& 1,\hfill \\ \hfill \frac{\partial \alpha }{\partial y}=& \frac{\partial \beta }{\partial x}=\frac{\partial \varphi }{\partial y}=0.\hfill \end{array}$$

(18)

Then,

$$\begin{array}{cc}\hfill {▵}_1=& \frac{P_{+}(\left(\frac{\partial \beta }{\partial y}\right)\left(\frac{\partial \varphi }{\partial x}\right)-\left(\frac{\partial \beta }{\partial x}\right)\left(\frac{\partial \varphi }{\partial y}\right)+\frac{\partial \varphi }{\partial x})+Q_{+}(\left(\frac{\partial \alpha }{\partial x}\right)\left(\frac{\partial \varphi }{\partial y}\right)-\left(\frac{\partial \alpha }{\partial y}\right)\left(\frac{\partial \varphi }{\partial x}\right)+\frac{\partial \varphi }{\partial y})}{P\left(\frac{\partial \varphi }{\partial x}\right)+Q\left(\frac{\partial \varphi }{\partial y}\right)}\hfill \\ \hfill =& \frac{P_{+}\left((x^T\left(T^{+}\right),y^T\left(T^{+}\right))\right)(1+\frac{\partial \beta }{\partial y})}{P(x^T\left(T\right),y^T\left(T\right))}\hfill \\ \hfill =& \frac{(1-p)(1-q)(r_0-d-a(1-p)ET)}{r_0-d-aET}.\hfill \end{array}$$

(19)

Let $C=\frac{r_0-d-a(1-p)ET}{a(1-p)ET},$ and we get the equation below:

$$\begin{array}{cc}\hfill {\int }_0^T{x}^T\left(t\right)dt=& \frac{1}{a}[ln(e^{(r_0-d)t}+C)]{|}_0^T.\hfill \end{array}$$

(20)

Thus,

$$\begin{array}{c}\hfill exp\left({\int }_0^T[\frac{\partial P}{\partial x}(x^T\left(t\right),y^T\left(t\right))+\frac{\partial Q}{\partial y}(x^T\left(t\right),y^T\left(t\right))]dt\right)\\ \hfill =exp((r_0+m-d)T+\frac{cb-2a}{a}[ln(e^{(r_0-d)t}+C)]{|}_0^T).\end{array}$$

(21)

Let $\theta =\frac{r_0-d-a(1-p)ET}{(1-p)(r_0-d-aET)},$ thus,

$$\begin{array}{cc}\hfill exp\left({\int }_0^T[\frac{\partial P}{\partial x}(x^T\left(t\right),y^T\left(t\right))+\frac{\partial Q}{\partial y}(x^T\left(t\right),y^T\left(t\right))]dt\right)=& {\theta }^{\frac{r_0+m-d}{r_0-d}}{\left(\frac{\theta +C}{1+C}\right)}^{\frac{cb-2a}{a}}\hfill \\ \hfill =& {\theta }^{\frac{r_0+m-d}{r_0-d}}{\left(\frac{r_0-d-a(1-p)ET}{(r_0-d-aET)}\right)}^{\frac{cb-2a}{a}}\hfill \\ \hfill =& {\theta }^{\frac{m}{r_0-d}+\frac{cb}{a}-1}{(1-p)}^{\frac{cb}{a}-2}.\hfill \end{array}$$

(22)

Therefore, the Floquet multiplier ${\mu }_2$ can be calculated as follows:

$$\begin{array}{cc}\hfill {\mu }_2=& {▵}_{1}exp\left({\int }_0^T[\frac{\partial P}{\partial x}(x^T\left(t\right),y^T\left(t\right))+\frac{\partial Q}{\partial y}(x^T\left(t\right),y^T\left(t\right))]dt\right)\hfill \\ \hfill =& (1-q){(1-p)}^{\frac{cb}{a}}{\theta }^{\frac{m}{r_0-d}+\frac{cb}{a}},\hfill \end{array}$$

(23)

where obviously, if $|{\mu }_2|<1,$ then the boundary order-1 periodic solution $(x^T\left(t\right),0)$ is stable, and if $|{\mu }_2|>1,$ the boundary order-1 periodic solution $(x^T\left(t\right),0)$ is unstable. □

Remark 1.

If $q=0,$ from (23), we get the following:

$${\mu }_2={(1-p)}^{\frac{cb}{a}}{\theta }^{\frac{m}{r_0-d}+\frac{cb}{a}}={(1-p)}^{-\frac{m}{r_0-d}}{\left(\frac{r_0-d-a(1-p)ET}{r_0-d-aET}\right)}^{\frac{m}{r_0-d}+\frac{cb}{a}}.$$

Since $R>1,$ then ${(1-p)}^{-\frac{m}{r_0-d}}>1;$ therefore, $|{\mu }_2|>1$ always holds, hence the boundary order-1 periodic solution $(x^T\left(t\right),0)$ of system (5) is always unstable.

3.2.3. Existence and Stability of the Order-1 Periodic Cycle for $\tau >0$

Theorem 3.

When $\tau >0$ and $|{\mu }_2|<1,$ the order-1 periodic solution $\left(\xi \right(t),\eta (t\left)\right)$ is orbitally asymptotically stable.

Proof. 

Let T denote the period of the order-1 periodic solution $\left(\xi \right(t),\eta (t\left)\right),$ and the order-1 periodic solution $\left(\xi \right(t),\eta (t\left)\right)$ passes point $M(ET,y_1)$ and point $M^{+}((1-p)ET,$$(1-q)y_1+\tau ).$

Using the calculation method of Theorem 2, we get ${\mu }_2,$

$$\begin{array}{cc}\hfill {\mu }_2=& {▵}_{1}exp\left({\int }_0^T[\frac{\partial P}{\partial x}(\xi \left(t\right),\eta \left(t\right))+\frac{\partial Q}{\partial y}(\xi \left(t\right),\eta \left(t\right))]dt\right)\hfill \\ \hfill =& (1-q)(1-p)\frac{\frac{r_0}{1+k((1-q)y_1+\tau )}-d-a(1-p)ET-b((1-q)y_1+\tau )}{\frac{r_0}{1+k{y}_1}-d-aET-b{y}_1}exp\left({\int }_0^{T}G\left(t\right)dt\right).\hfill \end{array}$$

(24)

Therefore, if $|{\mu }_2|<1,$ the order-1 periodic solution $\left(\xi \right(t),\eta (t\left)\right)$ is orbitally asymptotically stable, and $G\left(t\right)=\frac{r_0}{1+k\eta \left(t\right)}+(cb-2a)\xi \left(t\right)-(2h+b)\eta \left(t\right)+m-d.$ □

Theorem 4.

If $P_M\left(y_{p_{ET}}\right)<y_{p_{ET}},$ then the Poincaré map $P_M$ has only one stable fixed point $y_f$ on $[0,y_{p_{ET}}],$ and so system (5) has a globally stable order-1 periodic solution.

Proof. 

In Theorem 1, $P_M$ is represented as having a single fixed point. Therefore, system (5) has an order-1 periodic solution. In order to narrate the characteristics of the order-1 periodic solution for system (5) more accurately, firstly, the local stability of the fixed point of system (5) is analyzed. From the outcome of Theorem 3, we know that if $\tau >0$ and $|{\mu }_2|<1,$ the order-1 periodic solution for system (5) is locally asymptotically stable. Subsequently, the global asymptotic stability of the order-1 periodic solution for system (5) is discussed in the next two intervals: (a) $[0,y_f),$(b) $[y_f,+\infty ).$

In case (a), any trajectory running from $((1-q)ET,y_0^{+}),$$y_0^{+}\in [0,y_f)$ will get a series of pulse sequences undergoing multiple pulse processes; the pulse sequence is as follows: $y_0^{+}<y_1^{+}<y_2^{+}<y_3^{+}<\cdots <y_n^{+}<\cdots <y_f.$ In summary, the pulse sequence $P_M^n\left(y_0^{+}\right)$ increases monotonically with the increase of $n,$ and there is ${lim}_{n\to \infty }{P}_M^n\left(y_0^{+}\right)=y_f.$ In case (b), for any trajectory from $((1-q)ET,y_0^{+}),$$y_0^{+}\in [y_f,+\infty )$, after multiple pulses, there are two situations: (i) for all $n,$$P_M^n\left(y_0^{+}\right)>y_f$ holds; (ii) for all $n,$$P_M^n\left(y_0^{+}\right)>y_f$ does not hold.

In case (i), the pulse sequence $P_M^n\left(y_0^{+}\right)$ is obtained as follows: $y_f<\cdots <y_n^{+}<\cdots <y_3^{+}<y_2^{+}<y_1^{+}<y_0^{+}.$ Therefore, the pulse sequence $P_M^n\left(y_0^{+}\right)$ decreases monotonically with the increase of $n.$ Clearly, there is ${lim}_{n\to \infty }{P}_M^n\left(y_0^{+}\right)=y_f.$

In case (ii), it is surely feasible to find the smallest positive integer $n_1$ of $P_M^{n_1}\left(y_0^{+}\right)<y_f.$ That is, the trajectory with $y_0^{+}\in [y_f,+\infty )$ experiences a finite number of pulses, and then intersects with the straight line $x=(1-p)ET$ at point $((1-p)ET,P_M^{n_2}\left(y_0^{+}\right)),$ in which $P_M^{n_2}\left(y_0^{+}\right)\in [0,y_f).$

Therefore, with $n_2$ increasing, $P_M^{n_2}\left(y_0^{+}\right)$ is monotonically increasing. Eventually, we can get ${lim}_{n_2\to \infty }{P}_M^{n_2}\left(y_0^{+}\right)=y_f$ after going through countless pulses. From the above, it can be seen that the order-1 periodic solution for system (5) is globally asymptotically stable. □

Theorem 5.

If $P_M\left(y_{p_{ET}}\right)>y_{p_{ET}},$ then system (5) has an order-1 periodic solution; if $P_M\left(y_{p_{ET}}\right)>y_{p_{ET}},$ and any point $y_0^{+}\in [y_{p_{ET}},y_f)$ satisfies $P_M^2\left(y_0^{+}\right)>y_0^{+},$ then the order-1 periodic solution for system (5) is globally stable.

Proof. 

From Theorem 1, we know that if $P_M\left(y_{p_{ET}}\right)>y_{p_{ET}},$ then $P_M$ has a unique fixed point $y_f.$ For any $y_0^{+}\in [0,y_f),$ there is $P_M\left(y_0^{+}\right)>y_0^{+}$; for any $y_0^{+}\in (y_f,+\infty ),$$P_M\left(y_0^{+}\right)<y_0^{+}$ holds. According to the convexity of $P_M$, it can be known that $P_M$ takes the maximum value at point $y_{p_{ET}}$, where $y_{p_{ET}}\in (0,y_f).$

Firstly, we prove that if all $y_0^{+}\in [y_{p_{ET}},y_f)$ have $P_M^2\left(y_0^{+}\right)>y_0^{+}$, then the order-1 periodic solution of the system is globally asymptotically stable. Therefore, the global asymptotic stability of the order-1 periodic solution of the system is analyzed in the following three intervals: (a) $y_0^{+}\in [y_{p_{ET}},y_f);$ (b) $y_0^{+}\in (0,y_{p_{ET}});$ (c) $y_0^{+}\in (y_f,+\infty ).$

For case (a), any trajectory starting from $((1-q)ET,y_0^{+}),$$y_0^{+}\in [y_{p_{ET}},y_f)$ has $P_M\left(y_{p_{ET}}\right)$$>P_M\left(y_0^{+}\right)>y_f$ because of the monotonicity of $P_M.$ From Theorem 1, $P_M$ is monotonically increasing on $[y_{p_{ET}},y_f),$ so $P_M^2$ is also monotonically increasing on $[y_{p_{ET}},y_f)$, and $y_0^{+}<P_M^2\left(y_0^{+}\right)<y_f$ holds. By induction, we can get $P_M^{2(j-1)}\left(y_0^{+}\right)<P_M^{2j}\left(y_0^{+}\right)<y_f$ (j is a positive integer.) From the pulse sequence, we know that $P_M^{2j}\left(y_0^{+}\right)$ is monotonically increasing, and there are ${lim}_{j⟶+\infty }{P}_M^{2j}\left(y_0^{+}\right)=y_f$ ($j\ge 1$ and j is a positive integer). In the same way, the monotonicity of $P_M^{2j-1}\left(y_0^{+}\right)$ can be proved.

For case (b), it is known from the pulse sequence that there is a positive integer m that makes $P_M^m\left(y_0^{+}\right)\in [y_{p_{ET}},y_f)$ or $P_M^m\left(y_0^{+}\right)\in (y_f,+\infty )$ true. If $P_M^m\left(y_0^{+}\right)\in [y_{p_{ET}},y_f),$ referring to the method in case (a), $P_M^{m+2j}\left(y_0^{+}\right)$ is monotonically increasing, and ${lim}_{j⟶+\infty }{P}_M^{m+2j}\left(y_0^{+}\right)=y_f$ holds.

If $P_M^m\left(y_0^{+}\right)\in (y_f,+\infty ),$ it is known from the pulse relationship that there is $P_M^{m+1}\left(y_0^{+}\right)\in (y_{p_{ET}},y_f),$ and ${lim}_{j⟶+\infty }{P}_M^{m+1+2j}\left(y_0^{+}\right)=y_f.$

For case (c), from the discussion method of (b) in Theorem 4, it can be seen that any trajectory starting from point $((1-q)ET,y_0^{+}),$$y_0^{+}\in (y_f,+\infty )$ will approach a fixed point $y_f.$ In conclusion, if $P_M\left(y_{p_{ET}}\right)>y_{p_{ET}},$ and any point $y_0^{+}\in [y_{p_{ET}},y_f)$ satisfies $P_M^2\left(y_0^{+}\right)>y_0^{+},$ $P_M$ has an order-1 periodic solution which is globally asymptotically stable.

Finally, we prove that if there is a globally asymptotically stable periodic solution, then for any point $y_0^{+}\in [y_{p_{ET}},y_f)$, all holds $P_M^2\left(y_0^{+}\right)>y_0^{+}.$ Conversely, if there is at least one point ${\widehat{y}}^{+}\in [y_{p_{ET}},y_f)$ in $[y_{p_{ET}},y_f),$$P_M^2\left({\widehat{y}}^{+}\right)<{\widehat{y}}^{+}$ is true. In a small neighborhood of the fixed point $y_f,$ there is a point ${\widehat{y}}^{++}$ in which $P_M^2\left({\widehat{y}}^{++}\right)>{\widehat{y}}^{++}$ holds. Therefore, we know from the successor function that there is another point $y_f^{\ast }\in ({\widehat{y}}^{+},{\widehat{y}}^{++})$ that makes $P_M^2\left(y_f^{\ast }\right)=y_f^{\ast }$ true. This contradicts the global asymptotic stability of the order-1 periodic solution. Therefore, if there is a globally asymptotically stable periodic solution, it holds $P_M^2\left(y_0^{+}\right)>y_0^{+}$ for any point $y_0^{+}\in [y_{p_{ET}},y_f).$ □

3.3. Dynamic Analysis of System (5) for Case II

In Case II, we define impulse sets and phase sets in two sub-cases named Case IIa and Case IIb, respectively. Through the analysis of the system trajectories, we find that the definition and properties of the Poincaré map in Case IIa are the same as those in Case I. Then, we can obtain the same properties of the system as in Case I. Here, we do not demonstrate in detail, but we only give relevant conclusions for Case IIa and focus on Case IIb.

Theorem 6.

If $\tau =0,$${\mu }_2<1,$ then the boundary order-1 periodic solution $(x^T\left(t\right),0)$ for system (5) is stable. Furthermore, if $\tau =0,$ ${\mu }_2>1,$ the boundary order-1 periodic solution $(x^T\left(t\right),0)$ is unstable.

Theorem 7.

When $\tau >0$ and $|{\mu }_2|<1,$ the order-1 periodic solution $\left(\xi \right(t),\eta (t\left)\right)$ is orbitally asymptotically stable.

Theorem 8.

If $P_M\left(y_{p_{ET}}\right)<y_{p_{ET}},$ then the Poincaré map $P_M$ only has the globally asymptotically stable fixed point $y_f$ on the interval $[0,y_{p_{ET}}],$ that is, system (5) has an order-1 periodic solution.

Theorem 9.

If $P_M\left(y_{p_{ET}}\right)>y_{p_{ET}},$ system (5) has an order-1 periodic solution; if $P_M\left(y_{p_{ET}}\right)>y_{p_{ET}}$ and any point $y_0^{+}\in [y_{p_{ET}},y_f)$ satisfies $P_M^2\left(y_0^{+}\right)>y_0^{+},$ then the order-1 periodic solution for system (5) is globally stable.

Next, we analyze the characteristics of Poincaré mapping in case IIb. Different from Case I and Case IIa, the phase set of Case IIb is divided into two disjoint parts (see Figure 2c), such that the Poincaré map of Case IIb has different properties from Case I and Case IIa although their definitions are the same.

3.3.1. Properties of Poincaré Mapping for Case IIb

Theorem 10.

The Poincaré map $P_M$ has the following properties:

1 

The domain and range of $P_M$ are $D_2=\{[0,y_{T_1}]\cup [y_{T_2},+\infty )\}$ and $[\tau ,(1-q)y((1-p)ET,y_{T_1})+\tau ]$, respectively. It is increasing on $[0,y_{T_1}]$ and decreasing on $[y_{T_2},+\infty );$

2 

$P_M$ is continuously differentiable on $D_2;$

3 

$P_M$ has a unique fixed point $y_f$ or has no fixed point. If $\tau >0$ and $P_M\left(y_{T_1}\right)<y_{T_1},$ then $y_f\in (0,y_{T_1});$ if $\tau >0$ and $y_{T_1}<P_M\left(y_{T_1}\right)<y_{T_2},$ then the system has no fixed point; if $\tau >0$ and $P_M\left(y_{T_1}\right)>y_{T_2},$ then $y_f\in (y_{T_2},+\infty );$

4 

For all $y\in [0,+\infty ),$$P_M\left(y\right)$ has maximum value of $P_{M_{max}}=P_M\left(y_{T_1}\right)=P_M\left(y_{T_2}\right).$ $P_M$ has a minimum value of $P_{M_{min}}=P_M\left(0\right)=\tau .$

Proof. 

The proof is similar to the process in Theorem 1, and here we omit it. □

3.3.2. Existence and Stability of the Order-1 Limit Cycle for Case IIb

From Theorem 10, the Poincaré map has one fixed point or zero fixed points in case IIb. When the Poincaré map has a fixed point, then system (5) has an order-1 periodic solution. In the following, we will discuss the stability of the order-1 periodic solution for the system.

Theorem 11.

For case IIb, if $P_M\left(y_{T_1}\right)\le y_{T_1}$ or $P_M\left(y_{T_1}\right)\ge y_{T_2},$ then there exists a stable order-1 periodic solution for system (5). If $y_{T_1}<P_M\left(y_{T_1}\right)<y_{T_2},$ there is no order-1 periodic solution for system (5).

Proof. 

If $P_M\left(y_{T_1}\right)=y_{T_1}$ or $P_M\left(y_{T_1}\right)=y_{T_2},$ then $P_M$ has a fixed point on $(0,+\infty )$ according to Theorem 10.

If $P_M\left(y_{T_1}\right)<y_{T_1},$ then any trajectory initiating from $((1-p)ET,y_0^{+})$ experiences impulsive effects n times (infinitely), and we obtain $y_n^{+}\in [0,y_{T_1}].$ By Theorem 10, the impulsive point series $y_n^{+}$ is monotonically increasing on $[0,y_{T_1}]$ and monotonically decreasing on $[y_{T_2},+\infty ).$ Hence, we have ${lim}_{n\to \infty }{y}_n^{+}=y_f$ and $y_f\in [0,y_{T_1}].$ This proves that system (5) has a stable order-1 periodic solution.

If $P_M\left(y_{T_1}\right)>y_{T_2},$ then consider two cases: $y_{T_2}<\tau$ and $y_{T_2}>\tau .$ When $y_{T_2}<\tau ,$ then for any trajectory starting from $((1-q)ET,y_0^{+}),$ $y_0^{+}\in D_2$ infinitely passes many pulses, and from the properties of the Poincaré map, there is a unique fixed point on $[y_{T_2},+\infty ).$

When $y_{T_2}>\tau ,$ there exists a $y_c$ such that $P_M\left(y_c\right)=y_{T_2}.$ Moreover, any trajectory initiating from $((1-q)ET,y_0^{+}),$ $y_0^{+}\in [0,y_{T_2}]$ (or $y_0^{+}\in [y_c^1,y_c]$ and $0<y_c^1<y_c)$ will be free from impulsive effects after one single impulsive effect. Thus, system (5) has a stable order-1 periodic solution on $[y_{T_2},+\infty ).$

If $y_{T_1}<P_M\left(y_{T_1}\right)<y_{T_2},$ then any trajectory beginning with $((1-q)ET,y_0^{+}),$$y_0^{+}\in D_2$ will experience finite pulses and will then leave impulsive effects; finally, they will tend to the open set $(y_{T_1},y_{T_2}).$ That is, $P_M$ has no order-1 periodic solution. □

By using the same method as Theorem 5, we can prove the following theorem.

Theorem 12.

When $P_M\left(y_{T_1}\right)>y_{T_2}$ and any point $y_0^{+}\in [y_{T_2},y_f)$ satisfies $P_M^2\left(y_0^{+}\right)>y_0^{+},$ then the order-1 periodic solution of system (5) is globally stable.

4. Numerical Simulation

In this section, we employ some numerical simulations to verify the main results. The basic parameters are chosen as

$$r_0=1.5,k=0.03,a=0.02,d=0.03,b=0.04,c=0.5,m=0.1,h=0.01.$$

Then, we get the following system:

$$\left\{\begin{array}{c}\frac{dx\left(t\right)}{dt}=\frac{1.5x\left(t\right)}{1+0.03y\left(t\right)}-0.03x\left(t\right)-0.02x^2\left(t\right)-0.04x\left(t\right)y\left(t\right),\hfill \\ \frac{dy\left(t\right)}{dt}=0.02x\left(t\right)y\left(t\right)+0.1y\left(t\right)-0.01y^2\left(t\right).\hfill \end{array}\right.$$

(25)

Through calculation, we get $R=2.6834>1,$ then there exist four equilibria, i.e., $E_0=(0,0),$$E_1=(0,10),$$E_2(\frac{r_0-d}{a},0)=(73.5,0)$, and $E_3(x^{\ast },y^{\ast })=(4.996,20.139).$ By Theorems 2.1, 2.2, and 3.1 in Zhu et al. [37], we conclude that $E_0$ is an unstable node, $E_1$ and $E_2$ is an unstable saddle point, $E_3(x^{\ast },y^{\ast })=(4.996,20.139)$ is a stable focus.

Next, we will discuss the dynamics of the system under the effect of the state-dependent impulse control through the following cases.

For the Case I, we will firstly verify the existence and stability of the boundary order-1 periodic solution $(x^T\left(t\right),0).$ Let $p=0.5,ET=4,$ then we get the following system:

$$\left\{\begin{array}{c}\left.\begin{array}{c}\frac{dx\left(t\right)}{dt}=\frac{1.5x\left(t\right)}{1+0.03y\left(t\right)}-0.03x\left(t\right)-0.02x^2\left(t\right)-0.04x\left(t\right)y\left(t\right),\hfill \\ \frac{dy\left(t\right)}{dt}=0.02x\left(t\right)y\left(t\right)+0.1y\left(t\right)-0.01y^2\left(t\right).\hfill \end{array}\right\}x\ne 4,\hfill \\ \left.\begin{array}{c}\Delta x\left(t\right)=-0.5x\left(t\right),\hfill \\ \Delta y\left(t\right)=-qy\left(t\right)+\tau .\hfill \end{array}\right\}x=4.\hfill \end{array}\right.$$

(26)

On the one hand, if $\tau =0,$ let $q=0.2,$ and simple calculations show ${\mu }_2=0.86442<1,$ then by Theorem 2, we get that system (5) has a stable boundary order-1 periodic solution $(x^T\left(t\right),0)$ (see Figure 3). Let $q=0.01,$ by Theorem 2, there exists an unstable boundary order-1 periodic solution $(x^T\left(t\right),0)$ for system (5) (see Figure 4, where ${\mu }_2=1.06972>1$).

On the other hand, if $\tau >0,$ by calculation, we have $y_{p_{ET}}=21.177.$ Let us choose $\tau =3;$ simple calculations show that $P_M\left(y_{p_{ET}}\right)=(1-q)y(ET,y_{p_{ET}})+\tau =17.2104<y_{p_{ET}}$ holds. Then, by Theorem 4, the Poincaré map $P_M$ only has a globally asymptotically stable fixed point $y_f$ on $[0,y_{p_{ET}}]$ (see Figure 5a); as a result, system (5) has an order-1 periodic solution (see Figure 6a–c). If we let $\tau =10,$ then $P_M\left(y_{p_{ET}}\right)=(1-q)y(ET,y_{p_{ET}})+\tau =24.2104>y_{p_{ET}}$ holds, and by Theorem 5, the Poincaré map $P_M$ has a fixed point $y_f$ on $(y_{p_{ET}},+\infty )$(see Figure 5b). Therefore, system (5) has an order-1 periodic solution (see Figure 6d–f).

For case IIa, let $ET=6.4$ and $p=0.85,$ then we get the following system:

$$\left\{\begin{array}{c}\left.\begin{array}{c}\frac{dx\left(t\right)}{dt}=\frac{1.5x\left(t\right)}{1+0.03y\left(t\right)}-0.03x\left(t\right)-0.02x^2\left(t\right)-0.04x\left(t\right)y\left(t\right),\hfill \\ \frac{dy\left(t\right)}{dt}=0.02x\left(t\right)y\left(t\right)+0.1y\left(t\right)-0.01y^2\left(t\right).\hfill \end{array}\right\}x\ne 6.4,\hfill \\ \left.\begin{array}{c}\Delta x\left(t\right)=-0.85x\left(t\right),\hfill \\ \Delta y\left(t\right)=-qy\left(t\right)+\tau .\hfill \end{array}\right\}x=6.4.\hfill \end{array}\right.$$

(27)

Obviously, it is easy to calculate that $Q(ET,y_{ET})=(6.4,19.682),$$P_2(x_{p_2},y_{p_2})=(1.632,21.312).$ For the case of $\tau =0,$ let $q=0.3,$ and we compute that ${\mu }_2=0.86557<1;$ then, by Theorem 6, system (5) has a stable boundary periodic solution (see Figure 7). If we let $q=0.1,$ it is calculated that ${\mu }_2=1.11289>1,$ and by Theorem 6, the boundary periodic solution for system (5) is unstable (see Figure 8).

For the case of $\tau >0,$ for example, let $\tau =2,$ then we get $P_M\left(y_{p_{ET}}\right)=16.39405<y_{p_{ET}}=B=21.312.$ By Theorem 8, the Poincaré map $P_M$ has a globally asymptotically stable fixed point (see Figure 9a), and system (5) has a stable order-1 periodic solution (see Figure 10a–c)). If we let $\tau =10,$ then we get $P_M\left(y_{p_{ET}}\right)=24.38607>y_{p_{ET}}.$ It can be obtained from Theorem 9 that the Poincaré map $P_M$ has a globally asymptotically stable fixed point (see Figure 9b), and system (5) has an order-1 periodic solution (see Figure 10d–f).

In Case IIb, let $ET=6.4,$$p=0.5$, and $q=0.2;$ we then have the following system:

$$\left\{\begin{array}{c}\left.\begin{array}{c}\frac{dx\left(t\right)}{dt}=\frac{1.5x\left(t\right)}{1+0.03y\left(t\right)}-0.03x\left(t\right)-0.02x^2\left(t\right)-0.04x\left(t\right)y\left(t\right),\hfill \\ \frac{dy\left(t\right)}{dt}=0.02x\left(t\right)y\left(t\right)+0.1y\left(t\right)-0.01y^2\left(t\right).\hfill \end{array}\right\}x\ne 6.4,\hfill \\ \left.\begin{array}{c}\Delta x\left(t\right)=-0.5x\left(t\right),\hfill \\ \Delta y\left(t\right)=-0.2y\left(t\right)+\tau .\hfill \end{array}\right\}x=6.4.\hfill \end{array}\right.$$

(28)

By direct calculation, we obtain $T_1((1-p)ET,y_{T_1})=(3.2,17.132),$$T_2((1-p)ET,y_{T_2})=(3.2,30.510).$

In this case, any trajectory beginning with $(x_0^{+},y_0^{+}),y_0^{+}\in I_s$ either experiences a finite number of pulses or an infinite number of pulses. Furthermore, we draw images of the Poincaré mapping under case IIb (see Figure 11). When $\tau =1$, we get $P_M\left(y_{T_1}\right)=16.41391<y_{T_1}=17.132.$ From Figure 11a, we get that the fixed point $y_f$ of system (5) is on $[0,17.132]$, which confirms Theorem 10, and system (5) has a globally stable order-1 periodic solution (see Figure 12a–c).

From Figure 11b, when $\tau =5>0$ and $17.132<P_M\left(y_{T_1}\right)=P_M\left(17.132\right)=20.41717<30.510,$$P_M$ has no fixed point. Similarly, from Figure 12d–f, we know that there may be no periodic solution provided that $y_{T_1}<P_M\left(y_{T_1}\right)<y_{T_2}.$ From Figure 11c, when $\tau =25>0$ and $P_M\left(y_{T_1}\right)=P_M\left(17.132\right)=40.41717>30.510,$ the fixed point of $P_M$ is on $(y_{T_2},+\infty ).$ Figure 12g–i verify that the order-1 periodic solution of system (4) is stable when $P_M\left(y_{T_1}\right)=40.41717>y_{T_2}=30.510.$

Finally, we discuss the effect of the fear coefficient k on the dynamics of the system by numerical simulations. Since $R>1,$ we can then conclude that $k<\frac{1}{m}(\frac{r_0{h}^2}{dh+mb}-1).$ Let us take the following as parameters: $r_0=10,d=0.3,a=0.07,b=0.1,c=0.6,m=0.1,$$h=0.03$, $p=0.3,q=0.6,ET=4,\tau =10;$ then we get $k<4.4368.$ Numerical simulations show that for $0<k<1.83,$ system (5) has a stable order-1 periodic solution and for $k>1.83,$ the order-1 periodic solution disappears, and the solutions of system (5) tend to the positive equilibrium (see Figure 13).

Moreover, Figure 14 shows that when system (5) has a stable order-1 periodic solution, as the value of k increases, the oscillation of the prey population intensifies.

5. Conclusions

State-dependent impulse control has been applied to many aspects, and subsequently, it has been used in population ecology to study changes in biological populations in the last 20 years. Few papers have considered the problem of state feedback control with the fear effect on the prey or predator with additional food resources [61]. In reference [61], the author considered the state feedback control of a predator–prey model with the fear effect on prey; the model is as follows:

$$\left\{\begin{array}{c}\left.\begin{array}{c}\frac{dx\left(t\right)}{dt}=\frac{bx\left(t\right)}{1+ky\left(t\right)}-dx\left(t\right)-c{x}^2\left(t\right)-\frac{px\left(t\right)y\left(t\right)}{1+h_{1}x\left(t\right)},\hfill \\ \frac{dy\left(t\right)}{dt}=\frac{epx\left(t\right)y\left(t\right)}{1+(h_1+p{h}_2)x\left(t\right)}-my\left(t\right),\hfill \end{array}\right\}x\ne l,\hfill \\ \left.\begin{array}{c}\Delta x\left(t\right)=-q_{1}x\left(t\right),\hfill \\ \Delta y\left(t\right)=-q_{2}y\left(t\right)+\tau .\hfill \end{array}\right\}x=l,\hfill \end{array}\right.$$

(29)

where x and y are the prey and the predator, respectively, b is the birth rate of the prey, k represents the fear effect of the predator on the prey, d represents the natural mortality rate of the prey, c is the mortality due to intraspecific competition of the prey, p is the predator for the predation rate of the prey, $h_1$ represents the time spent by the predator to capture the prey, $h_2$ represents the time spent on digesting the captured prey, e is the utilization rate of the predator by the predator, and m is the death rate of the predator. $q_1$ and $q_2$ are the catchability coefficient of prey and predator, respectively. l represents the optimal capture level of the prey. The author assumed that when the number of prey reaches the optimal capture level l, the prey and predator are harvested. In order to avoid long-term capture leading to the extinction of the predator population, an appropriate amount of $\tau$ predators is placed. For system (29), by establishing the corresponding Poincaré map, the authors proved the existence of the order-k periodic solutions for the system.

In system (29), the authors took into account the fear effect of predators and did not consider the additional food resources of predators, which created a model that ignores the effect of impulses simpler in structure than our model. More specifically, our model generates a new equilibrium point $E_2(0,\frac{m}{h}),$ that although the prey eventually becomes extinct, because the predator has additional food resources, the predator may eventually exist. Therefore, the dynamics of our model are more complex and difficult to study for the problem of impulse feedback control, with both the fear effect on the prey and the additional food resources for the predator.

In this paper, we developed the application of state-dependent impulses by applying it to a predator–prey model, with the fear effect on the prey population and the additional food source for the predator population. In the model we built, we considered controlling the number of prey and predator populations based on the number of prey populations, hoping to achieve a certain balance of predators and prey. We constructed the Poincaré map of the state-dependent impulse system, and by analyzing the properties of Poincaré maps of pulse sets and phase sets at different positions, the dynamic properties of the system, including the existence, nonexistence, and stability of periodic solutions were investigated.

We found that the existence of the periodic solution for the system completely depends on the property of the Poincaré map. To verify the theoretical results we obtained, we gave some examples and numerical simulations for various cases. At the same time, numerical simulations also showed that the influence of the fear factor on the system is significant, that the increase of the fear factor makes the system periodic solution disappear, and that the predator and prey populations gradually stabilize to a level near the positive equilibrium.

Compared with the model without the state feedback control, the introduction of state feedback control makes the system dynamics more complicated. In [37], when $R<1,$ the system only had one stable equilibrium $E_1;$ biologically, the pest population x would eventually tend to extinction, while the natural enemy population y would eventually tend to a constant $\frac{m}{h}.$ Although the pests were eventually eradicated, this was a departure from the original purpose of integrated pest management. Note that when $R<1,$ the equilibrium $E_2$ is an unstable saddle point. Biologically, because the natural enemy population y has sufficient food sources, including pests and additional food resources, the population y will not become extinct in the end. However, for system (5), we controlled the pest population through an integrated pest management strategy, that is, by spraying pesticides and releasing natural enemies to control the pests. When $\tau =0,$ pests are controlled only by spraying pesticides. According to Theorems 1 and 2, when $\tau =0$ and ${\mu }_2<1,$ there is a stable boundary order-1 periodic solution for system (5). Biologically, the natural enemy population y may become extinct due to the overspraying of pesticides. When $R>1,$ system (4) has a stable equilibrium $E_3$ (stable node or focus); biologically, the pest population x and the natural enemy population y eventually tend to constants $x^{\ast }$ and $y^{\ast }$, respectively. In order to achieve the maximum benefit of economic harvest without destroying the ecological balance, we adopted an integrated pest management strategy, that is, when the pest population x reached $ET,$ the pests were controlled by spraying pesticides and releasing natural enemies. According to Theorems 4–9, there exists a globally asymptotically stable order-1 periodic solution for system (5). Biologically, the comprehensive implementation of pest control strategies makes the pest population x and natural enemy population y exhibit periodic oscillations.