Dynamic Behaviours and Bifurcation Analysis of a Filippov Ecosystem with Fear Effect

Abstract

Evidences of the fear effect of the prey are well documented, which can greatly affect the dynamics of the predator-prey system. In this study, considering that the fear effect of the prey is triggered on as the density of the predator reaches and exceeds a threshold value, we develop a Filippov system of predator-prey model with the fear effect. In addition, we also include a modify factor of the growth rate of the prey when they adopt the antipredator behaviours due to the fear effect. We initially analyze the dynamics of the two subsystems, including the existence and stability of the equilibria. Utilizing the theory of the Filippov system, we discuss the sliding dynamics, i.e., the existence of sliding region and sliding equilibria. By choosing the threshold as the bifurcation parameter, we investigate the bifurcations near the regular equilibria. The solution curve has three cases: crossing the threshold curve, sliding on the threshold curve, and approaching the pseudoequilibrium. Finally, we numerically verified the existence of the global sliding bifurcation near the regular equilibrium and also the touching bifurcation.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Traditionally, the researchers usually connect the predator and their prey only through direct predations [1, 2]. In the recent decades, lots of experiments show that the relationship between the predator and the prey is far more complex than hunting. For one aspect, due to the fear effect of the predation [3–5], the prey will choose to adopt the antipredator behaviours, such as moving to new habitats, searching for safer foods, and other physiological changes, aiming at avoiding the predation. This will definitely help to decrease the death rate of the prey. However, the growth rate of the preys will also be decreased because of the extra costs of adopting antipredator behaviours. Therefore, there could be the trade-off of the fear effect on the growth of the prey. Also, it will be important to study the impacts of the fear effect of prey on the dynamic behaviours between the prey and their predators.

Recently, some scholars have also considered the influence of fear effect on the dynamic behaviours of predator-prey models. For example, Wang et al. [6] first proposed a predator-prey model with the fear effect. In their paper, they have demonstrated a correlation between the fear effect and the direct effect of predators on prey and gave three specific expressions for the fear effect. They also demonstrated that the fear effect does not affect the dynamic behaviours of the system with the Holling-I response function. However, when the functional response is Holling-II, higher levels of fear can stabilize the system by eliminating the presence of oscillatory behaviours. When the appropriate birth rate and fear effect are selected, the system will have a limit cycle. The authors in [7] established a Beddington–DeAngelis predator-prey model with fear effect, refuge, and harvest. They found that there is a critical value of the fear effect, and the two species can continue to exist for a long time through positive density levels. There are also other critical values so that both species can exist at a positive density, but their density levels swing periodically over time. Zhang et al. [8] proved that the system exhibits multiple spatiotemporal patterns due to spatial memory delay and nonlocal fear effect delay. Researchers explored the impact of the fear effect and the Leslie–Gower function on the dynamical behaviours of the predator-prey model and found that as the fear effect increases, the dynamical behaviours of the system switches multiple times until the prey eventually becomes extinct, while the predator survives due to the presence of alternative prey in [9]. For more relevant research on the fear effect, refer to the literature [10–14].

When the number of predators reaches a certain level, the prey will have a fear effect and show antipredator behaviours, resulting in a change in the functional response between the predator and the prey. Therefore, we need to establish a model with threshold conditions. But most predator-prey models are described by ordinary differential equations with continuous right-hand sides, but these models cannot reflect the influence of refuge, group defense, and other factors on population dynamics. Therefore, a discontinuous right-hand nonsmooth dynamic system widely used in mechanics has attracted the attention of ecologists [15–19]. This kind of system is called the Filippov system or switched system, which provides a basic framework for the establishment of many mathematical models with practical significance [20–22]. Therefore, based on the above research, this study proposes a Filippov predator-prey model, which depends on the number of predators using the threshold strategy. Our model extends the existing predator-prey model with the fear effect by introducing a threshold strategy to describe the impact of the fear effect when the number of predators exceeds the threshold. When the number of predators is below the threshold, the functional response is Holling-I. When the number of predators is above the threshold, the functional response changes to Holling-II, and the fear effect is produced.

This study is organized as follows: In Section 2, a Filippov predator-prey model with fear effect is established, and the basic theory and related definitions of the Filippov system will be provided. In Section 3, the dynamic behaviours of the two subsystems are analyzed. In Section 4, the sliding region, sliding dynamics, and the existence of various equilibria of the model are analyzed theoretically. In Section 5, the regular/virtual equilibrium bifurcation, boundary equilibrium bifurcation, and global sliding bifurcation are numerically verified. Finally, the conclusions are presented in Section 6.

2. Model Formulation and Preliminaries

2.1. Model Formulation

We use the classical Lotka–Volterra model to describe the interaction between the predator and prey [23], i.e.,$$\begin{array}{}\text{(1)}& \begin{array}{l}\frac{dx}{dt}=rx-dx-a{x}^2-pxy\text{,}\\ \frac{dy}{dt}=upxy-my\text{,}\end{array}\end{array}$$where $x$ and $y$ represent the population density of prey and predator, respectively. In model (1), the growth of prey follows the logistic mode, and predator-prey interactions follow the Holling-I functional response. $r$ and $d$ are the natural birth rate and natural mortality rate of the prey, respectively, $a$ indicates the decay rate of the prey due to intraspecific competition, and $p$ denotes the rate of search for the prey by the predator. The rate of conversion of prey biomass to predator biomass is denoted by $u$, and $m$ expresses the natural mortality of the predator. In addition, the assumption that $r>d$ will apply to the entire article.

We assume that the fear effect of the prey to their predator is triggered on by a threshold value of the density of predators $Y_T$. That is, there is no fear effect before the predator reaches the threshold value $Y_T$, and the interactions between the predator and prey is assumed to follow model (1). When $y>Y_T$, due to the fear effect, the prey may adopt the antipredator behaviours, such as moving to new habitats or searching for safer foods and other physiological changes, aiming at avoiding the predation; consequently, there should be a limitation of the predation of the predator to the prey. Hence, we use the Holling-II response function as the existence of fear effects of the prey. In addition, the antipredator behaviours will also lead to the extra costs of the prey, reflecting on the decreasing of its growth rate. Also, we use the factor to modify the growth rate of prey due to the antipredator behaviours induced by the fear effect. Then, the predator-prey model with the fear effect becomes [6]$$\begin{array}{}\text{(2)}& \begin{array}{l}\frac{dx}{dt}=\varphi k,yrx-dx-a{x}^2-\frac{pxy}{1+qx}\text{,}\\ \frac{dy}{dt}=\frac{upxy}{1+qx}-my\text{,}\end{array}\end{array}$$where $q$ is the half-saturation constant and $\varphi k,y$ is a decreasing function with respect to $k$ and $y$. $k$ reflects the degree of fear that drives antipredatory behaviours in prey. According to the biological significance of $k$, $y$, and $\varphi k,y$, the function $\varphi k,y$ should satisfy the following four assumptions:

(1) $\varphi 0,y=1$: if there is no fear, there will be no impact on prey. In other words, the natural birth rate of the prey population will not decrease

(2) $\varphi k,0=1$: if there are no predators, there is no fear effect on prey, and therefore the natural birth rate of prey does not decrease

(3) $\partial \varphi k,y/\partial k<0$: the natural birth rate of prey decreases with the increasing fear effect

(4) $\partial \varphi k,y/\partial y<0$: if the density of predators increases, the fear effect increases, and therefore the natural birth rate of prey decreases

A simple form of $\varphi k,y$ with the above assumptions satisfied is $\varphi k,y=1/1+ky$, proposed by [6]. So, system (2) can be rewritten as follows:$$\begin{array}{}\text{(3)}& \begin{array}{l}\frac{dx}{dt}=\frac{rx}{1+ky}-dx-a{x}^2-\frac{pxy}{1+qx}\text{,}\\ \frac{dy}{dt}=\frac{upxy}{1+qx}-my\text{.}\end{array}\end{array}$$

So, systems (2) and (3) can be integrated into the following nonsmooth dynamical system:$$\begin{array}{}\text{(4)}& \begin{array}{l}\frac{dx}{dt}=\frac{rx}{1+\epsilon ky}-dx-a{x}^2-\frac{pxy}{1+\epsilon qx}\text{,}\\ \frac{dy}{dt}=\frac{upxy}{1+\epsilon qx}-my\text{,}\end{array}\end{array}$$with$$\begin{array}{}\text{(5)}& \epsilon =\begin{array}{l}0,\mathrm{y}<Y_T,\\ 1,\mathrm{y}>Y_T.\end{array}\end{array}$$

2.2. Preliminaries

To further analyze the Filippov system, here are some symbols [24–26]. Let $R_{+}^2=Z={x,y}^{T}x\ge 0,y\ge 0,F_{S_i}:R_{+}^2⟶R_{+}^2,i=\mathrm{1,2}$ be the smooth vector fields and $\epsilon$ be a partition function that depends on $y-Y_T$. Let $HZ=y-Y_T$ be a smooth scalar function with a nonzero gradient whose threshold value depends on the number of individuals in the predator. Thus, (5) can be rewritten as given in the following equation:$$\begin{array}{}\text{(6)}& \epsilon =\begin{array}{l}0,HZ<0,\\ 1,HZ>0.\end{array}\end{array}$$

For convenience, we denote$$\begin{array}{c}\text{(7)}& F_{S_1}={rx-dx-a{x}^2-pxy,upxy-my}^T\text{,}{F}_{S_2}={\frac{rx}{1+ky}-dx-a{x}^2-\frac{pxy}{1+qx},\frac{upxy}{1+qx}-my}^T\text{.}\end{array}$$

Then, system (4) with (6) can be rewritten as the following Filippov system:$$\begin{array}{}\text{(8)}& \dot{Z}t=\begin{array}{l}{F}_{S_1}Z,Z\in S_1,\\ F_{S_2}Z,Z\in S_2,\end{array}\end{array}$$where $S_1=Z\in R_{+}^{2}HZ<0$ and $S_2=Z\in R_{+}^{2}HZ>0$. $\mathrm{\Sigma }$ be called as the switching manifold, which is the separation boundary of the region and the discontinuity boundary set $\mathrm{\Sigma }$ can be defined as$$\begin{array}{}\text{(9)}& \mathrm{\Sigma }=Z\in R_{+}^{2}HZ=0\text{.}\end{array}$$

Furthermore, $\mathrm{\Sigma }$ divides $R_{+}^2$ into two distinct regions $S_1$ and $S_2$. Then, the Filippov system (8) in region $S_i$ is called subsystem $S_i$, $i=\mathrm{1,2}$.

Let$$\begin{array}{}\text{(10)}& \sigma Z=H_{Z}Z,F_{S_1}Z\cdot H_{Z}Z,F_{S_2}Z=F_{S_1}HZ\cdot F_{S_2}HZ\text{,}\end{array}$$where $\cdot$ denotes the standard scalar product, and Lie derivative $F_{S_i}HZ=H_{Z}Z,F_{S_i}Zi=\mathrm{1,2}$ is used to denote the directional derivative of $H$ concerning the vector field $F_{S_i}i=\mathrm{1,2}$ at $Z$. Therefore, we partition the discontinuity boundary according to the sign case of $\sigma Z$.

Accordingly, the boundary will be classified as follows:

(i) ${\mathrm{\Sigma }}_e\in \mathrm{\Sigma }$ is called the escaping region if $F_{S_1}HZ<0\mathrm{a}\mathrm{n}\mathrm{d}{F}_{S_2}HZ>0$ on ${\mathrm{\Sigma }}_e$

(ii) ${\mathrm{\Sigma }}_s\in \mathrm{\Sigma }$ is called the sliding region if $F_{S_1}HZ>0\mathrm{a}\mathrm{n}\mathrm{d}{F}_{S_2}HZ<0$ on ${\mathrm{\Sigma }}_s$

(iii) ${\mathrm{\Sigma }}_c\in \mathrm{\Sigma }$ is called the crossing region if $F_{S_1}HZ\cdot F_{S_2}HZ>0$ on ${\mathrm{\Sigma }}_c$

For the convenience of the later study, it is useful to list the definitions of the various equilibria of the Filippov system [16, 17] for this paper.

Definition 1.

If a point $Z^{\ast }$ makes $F_{S_1}{Z}^{\ast }=0,H{Z}^{\ast }<0$ or $F_{S_2}{Z}^{\ast }=0,H{Z}^{\ast }>0$, the point $Z^{\ast }$ is called the regular equilibrium of the Filippov system (8). If a point $Z^{\ast }$ makes $F_{S_1}{Z}^{\ast }=0,H{Z}^{\ast }>0$ or $F_{S_2}{Z}^{\ast }=0,H{Z}^{\ast }<0$, the point $Z^{\ast }$ is called the virtual equilibrium of the Filippov system (8).

Definition 2.

If a point $Z^{\ast }$ is the equilibrium of the sliding region of the Filippov system (8), the point $Z^{\ast }$ is called pseudoequilibrium and $Z^{\ast }\in {\mathrm{\Sigma }}_s$, that is, $1-\lambda F_{S_1}{Z}^{\ast }+\lambda F_{S_2}{Z}^{\ast }=0,H{Z}^{\ast }=0$, and $0<\lambda <1$, where$$\begin{array}{}\text{(11)}& \lambda =\frac{H_{Z}Z,F_{S_1}Z}{H_{Z}Z,F_{S_1}Z-F_{S_2}Z}\text{.}\end{array}$$

Definition 3.

If $F_{S_1}{Z}^{\ast }=0,H{Z}^{\ast }=0$ or $F_{S_2}{Z}^{\ast }=0,H{Z}^{\ast }=0$, the point $Z^{\ast }$ is called a boundary equilibrium of the Filippov system (8).

Definition 4.

If $Z^{\ast }\in {\mathrm{\Sigma }}_s$, $H{Z}^{\ast }=0$, and $H_Z{\mathrm{Z}}^{\ast },F_{S_1}{\mathrm{Z}}^{\ast }=0$ or $H_Z{\mathrm{Z}}^{\ast },F_{S_2}{\mathrm{Z}}^{\ast }=0$, the point $Z^{\ast }$ is called a tangent point of the Filippov system (8).

3. Dynamics Analysis of Subsystems

The Filippov system (8) consists of three parts: two smooth subsystems and a discontinuity boundary, so it is necessary to study the dynamical behaviours of two subsystems before analyzing the complete Filippov system. In this section, we analyze the existence and stability of the equilibria of the two subsystems separately.

3.1. Dynamics of the Subsystem $S_1$

When $y<Y_T$, the dynamical behaviours of system (8) is determined by the subsystems defined by (1). The subsystem $S_1$ has the following three equilibria: $E_{01}=\mathrm{0,0}$, $E_{11}=r-d/a,0$, and $E_1^{\ast }=x_1^{\ast },y_1^{\ast }$ where $x_1^{\ast }=m/\mathrm{u}\mathrm{p},y_1^{\ast }=upr-d-ma/\mathrm{u}{\mathrm{p}}^2$. The equilibria $E_{01}$ and $E_{11}$ always exist. The point $E_1^{\ast }$ is in the first quadrant if $r>d+am/up$. The following theorem will give the possible dynamical behaviours of all equilibria and the conditions under which they occur.

Theorem 5.

(I) The trivial equilibrium $E_{01}=\mathrm{0,0}$ is a saddle

(II) The boundary equilibrium $E_{11}=r-d/a,0$ is locally asymptotically stable if$$\begin{array}{}\text{(12)}& r<d+\frac{am}{up}\text{.}\end{array}$$

(III) The interior equilibrium $E_1^{\ast }=x_1^{\ast },y_1^{\ast }$ is locally asymptotically stable when (12) is violated

Proof.

We analyze the stability of the equilibria by analyzing the corresponding Jacobian matrix for each equilibrium. Now, the Jacobian matrix at $E_{01}=\mathrm{0,0}$, $E_{11}=r-d/a,0$, and $E_1^{\ast }=x_1^{\ast },y_1^{\ast }$ are given by$$\begin{array}{}\text{(13)}& J_{E_{01}}=\begin{array}{ll}r-d& 0\\ 0& -m\end{array},J_{E_{11}}=\begin{array}{ll}d-r& \frac{pd-r}{a}\\ 0& \frac{\mathrm{u}\mathrm{p}r-d}{a}-m\end{array},\text{and}{J}_{E_1^{\ast }}=\begin{array}{ll}r-d-2a{x}_1^{\ast }-p{y}_1^{\ast }& -p{x}_1^{\ast }\\ up{y}_1^{\ast }& up{x}_1^{\ast }-m\end{array}\text{.}\end{array}$$

Since $tr{J}_{E_{01}}=r-d-m$ and $\det J_{E_{01}}=-mr-d<0$, the trivial equilibrium $E_{01}$ is a saddle.

The trace and determinant corresponding to the matrix $J_{E_{11}}$ are $tr{J}_{E_{11}}=up-ar-d/a-m$ and $\det J_{E_{11}}=d-rupr-d/a-m$. Both $tr{J}_{E_{11}}$ and $\det J_{E_{11}}$ can be positive or negative, therefore, from the Routh–Hurwitz criterion, it follows that $E_{11}$ is locally asymptotically stable when $tr{J}_{E_{11}}<0$ and $\det J_{E_{11}}>0$. Thus, the boundary equilibrium is locally asymptotically stable when $r<d+am/up$. Similarly, the condition of local asymptotic stability of internal equilibrium $E_1^{\ast }$ is $r>d+am/up$.

Theorem 6.

When $r\in d,d+am/up$, the boundary equilibrium $E_{11}$ is globally asymptotically stable. When $r\in d+am/up,+\infty$, the interior equilibrium $E_1^{\ast }$ is globally asymptotically stable.

Proof.

The two functions on the right-hand side of system (1) are denoted by $Px,y$ and $Qx,y$. We consider the Dulac function as $Bx,y=1/xy$. After calculations, we obtain$$\begin{array}{}\text{(14)}& D=\frac{\partial PB}{\partial x}+\frac{\partial QB}{\partial y}=-\frac{a}{y}<0\text{.}\end{array}$$

For $x,y\in \mathrm{int}{R}_{+}^2=0,\infty \times 0,\infty$. Thus, according to the Dulac–Bendixson theorem, system (1) has no periodic orbits in $0,\infty \times 0,\infty$. Besides, when $r>d+am/up$, the positive equilibrium in $0,\infty \times 0,\infty$ has only the interior equilibrium $E_1^{\ast }$. So all positive solutions will tend to $E_1^{\ast }$. Combining the $E_1^{\ast }$ proved in Theorem 5 is locally asymptotically stable, it can be concluded that $E_1^{\ast }$ is globally asymptotically stable if condition (12) does not hold.

When $r\in d,d+am/up$, there are only two equilibria $E_{01}$ and $E_{11}$ in $R_{+}^2$. In addition, there is no periodic orbit in $R_{+}^2$ which means that every positive solution converges to $E_{01}$ or $E_{11}$. It is easy to obtain that $E_{01}$ is repelling when $r>d$, so every positive solution converges to $E_{11}$. Combining the $E_{11}$ proved in Theorem 5 is locally asymptotically stable, it can be concluded that $E_{11}$ is globally asymptotically stable if $r\in d,d+am/up$.

3.2. Dynamics of the Subsystem $S_2$

For the subsystem $S_2$, it is easy to calculate the following three equilibria: $E_{02}=\mathrm{0,0},E_{12}=r-d/a,0$ and $E_2^{\ast }=x_2^{\ast },y_2^{\ast }$ where $x_2^{\ast }=m/up-mq,y_2^{\ast }=-C_2+\sqrt{C_2^2-4C_1{C}_3}/2C_1$ with$$\begin{array}{}\text{(15)}& C_1=k{up-mq}^2,C_2={up-mq}^2+ukam+dup-dm,and{C}_3=uam+d-rup-mq\text{.}\end{array}$$

Notice that the equilibria $E_{01}$ and $E_{11}$ are always present and the equilibrium $E_2^{\ast }$ is feasible if $r>d+am/up-mq$ with $up>mq$.

Theorem 7.

The trivial equilibrium $E_{02}=\mathrm{0,0}$ is a saddle.

Proof.

The Jacobian matrix at $E_{02}=\mathrm{0,0}$ is$$\begin{array}{}\text{(16)}& J_{E_{02}}=\begin{array}{ll}r-d& 0\\ 0& -m\end{array}\text{.}\end{array}$$

Obviously, we can get $tr{J}_{E_{02}}=r-d-m$ and $\det J_{E_{02}}=-mr-d<0$, so by the Routh–Hurwitz criterion it follows that the trivial equilibrium $E_{02}$ is a saddle.

The following lemmas [6] give the stability analysis regarding the boundary equilibrium $E_{02}$ and the internal equilibrium $E_2^{\ast }$.

Lemma 8.

The boundary equilibrium $E_{02}$ is locally asymptotically stable if$$\begin{array}{}\text{(17)}& r-dup-mq<am,\end{array}$$is satisfied and is unstable if$$\begin{array}{}\text{(18)}& r-dup-mq>am,\end{array}$$holds.

Lemma 9.

The boundary equilibrium $E_{02}$ is globally asymptotically stable if $up\le mq$ holds.

Lemma 10.

The internal equilibrium $E_2^{\ast }$ is locally asymptotically stable if$$\begin{array}{}\text{(19)}& \begin{array}{l}r>d+\frac{am}{up-mq}\text{,}\\ r\le d+\frac{aup+mq}{qup-mq}\text{,}\\ up>mq\text{,}\end{array}\end{array}$$or$$\begin{array}{}\text{(20)}& \begin{array}{l}r>d+\frac{aup+mq}{qup-mq}\text{,}\\ k>\frac{q{up-mq}^{2}r-dqup-mq-aup+mq}{u^{2}paqdup-mq+aup+mq}\text{,}\\ up>mq\text{.}\end{array}\end{array}$$

It is unstable if$$\begin{array}{}\text{(21)}& \begin{array}{l}r>d+\frac{aup+mq}{qup-mq}\text{,}\\ k<\frac{q{up-mq}^{2}r-dqup-mq-aup+mq}{u^{2}paqdup-mq+aup+mq}\text{,}\\ up>mq\text{.}\end{array}\end{array}$$

Lemma 11.

The internal equilibrium $E_2^{\ast }$ is globally asymptotically stable if$$\begin{array}{}\text{(22)}& \begin{array}{l}r>d+\frac{am}{up-mq}\text{,}\\ r\le d+\frac{aup+mq}{qup-mq}\text{,}\\ up>mq\text{,}\\ up-mq\le rq\text{.}\end{array}\end{array}$$

4. Sliding Region and Equilibria of the Filippov System (8)

In this section, we calculate the sliding region and various equilibria of the Filippov system (8) according to the definitions.

4.1. Sliding Segment and Region

Through simple calculation, the following Lie derivatives can be obtained:$$\begin{array}{}\text{(23)}& F_{S_1}HZ=Y_{T}upx-m\text{,}\text{(24)}& F_{S_1}HZ=Y_T\frac{upx}{1+qx}-m\text{.}\end{array}$$Let the right-hand sides of the above two functions be zero to get $x_{S_1}=m/up$ and $x_{S_2}=m/up-mq$. By definitions, we can obtain the following theorem for the sliding region, escaping region, and crossing region.

Theorem 12.

(I) If $x_{S_2}>0$, the sliding and crossing segment of the Filippov system (8) can be defined as$$\begin{array}{c}\text{(25)}& {\mathrm{\Sigma }}_s=x,y\in \mathrm{\Sigma }{x}_{S_1}<x<x_{S_2}\text{,}{\mathrm{\Sigma }}_c=x,y\in \mathrm{\Sigma }0<x<x_{S_1}\cup x,y\in \mathrm{\Sigma }{x}_{S_2}<x\text{.}\end{array}$$

(II) If $x_{S_2}<0$, the escaping and crossing segment of the Filippov system (8) can be defined as$$\begin{array}{c}\text{(26)}& {\mathrm{\Sigma }}_e=x,y\in \mathrm{\Sigma }0<x<x_{S_1}\text{,}{\mathrm{\Sigma }}_c=x,y\in \mathrm{\Sigma }{x}_{S_1}<x\text{.}\end{array}$$

Proof.

(I) If $x_{S_2}>0$, we can know that $0<x_{S_1}<x_{S_2}$. According to (23), (24), and the definition of sliding region, we obtain the sliding region$$\begin{array}{}\text{(27)}& {\mathrm{\Sigma }}_s=x,y\in \mathrm{\Sigma }{x}_{S_1}<x<x_{S_2}\text{.}\end{array}$$

  Moreover, there are two possible scenarios that can lead to the presence of transversal intersection and crossing region:$$\begin{array}{}\text{(28)}& F_{S_1}HZ>0\mathrm{a}\mathrm{n}\mathrm{d}{F}_{S_2}HZ>0,\end{array}$$

  or$$\begin{array}{}\text{(29)}& F_{S_1}HZ<0\mathrm{a}\mathrm{n}\mathrm{d}{F}_{S_2}HZ<0\text{.}\end{array}$$

  From (28), we calculate that $x<x_{S_1}$ and $x<x_{S_2}$, i.e., $x<\min x_{S_1},x_{S_2}$. So, we get the crossing region ${\mathrm{\Sigma }}_{c1}=x,y\in \mathrm{\Sigma }0<x<x_{S_1}$. Also, from (29), we can get $x>x_{S_1}$ and $x>x_{S_2}$, i.e., $x>\max x_{S_1},x_{S_2}$. Thus, the crossing region is ${\mathrm{\Sigma }}_{c2}=x,y\in \mathrm{\Sigma }{x}_{S_2}<x$. In summary, the crossing region is$$\begin{array}{}\text{(30)}& {\mathrm{\Sigma }}_c=x,y\in \mathrm{\Sigma }0<x<x_{S_1}\cup x,y\in \mathrm{\Sigma }{x}_{S_2}<x\text{.}\end{array}$$

(II) If $x_{S_2}<0$, we can know that $x_{S_2}<0<x_{S_1}$. By definition, we can get $x<x_{S_1}\mathrm{a}\mathrm{n}\mathrm{d}x>x_{S_2}$. Thus, the escaping region is$$\begin{array}{}\text{(31)}& {\mathrm{\Sigma }}_e=x,y\in \mathrm{\Sigma }0<x<x_{S_1}\text{.}\end{array}$$

Similarly, (28) and (29) can be obtained for the crossing region. Simplifying (28) can yield $x<x_{S_1}$ and $x<x_{S_2}$, i.e., $x<\min x_{S_1},x_{S_2}$. To ensure that the solution must satisfy biological significance, $x$ must be a non-negative constant, so this situation does not exist. For (29), the calculation gives $x>x_{S_1}$ and $x>x_{S_2}$, i.e., $x>\max x_{S_1},x_{S_2}$. So the crossing region is$$\begin{array}{}\text{(32)}& {\mathrm{\Sigma }}_c=x,y\in \mathrm{\Sigma }{x}_{S_1}<x\text{.}\end{array}$$

We next examine the existence of four types of equilibria of the Filippov system (8) according to the definitions in Section 2.2.

4.2. Pseudoequilibrium

Based on Definition 2, we first calculate$$\begin{array}{}\text{(33)}& \lambda =\frac{upx-mqx+1}{pqu{x}^2}\text{.}\end{array}$$So,$$\begin{array}{}\text{(34)}& 1-\lambda =\frac{mq-upx+m}{pqu{x}^2}\text{.}\end{array}$$

Substituting the expressions for $\lambda$ and $1-\lambda$ into $1-\lambda F_{S_1}{Z}^{\ast }+\lambda F_{S_2}{Z}^{\ast }=0$, we can obtain the following equation:$$\begin{array}{}\text{(35)}& \begin{array}{l}\frac{dx}{dt}=\frac{C_4{x}^3+C_5{x}^2+C_{6}x+C_7}{pqux1+k{Y}_T}\text{,}\\ \frac{dy}{dt}=0\text{,}\end{array}\end{array}$$where $C_4=-apqu1+k{Y}_T$, $C_5=-pqud1+k{Y}_T-r$, $C_6=-mpq{Y}_{T}1+k{Y}_T+kr{Y}_{T}up-mq$, and $C_7=kmr{Y}_T$.

We know from [27] that we only need to discuss the case of the molecular roots of $dx/dt$.

Let$$\begin{array}{}\text{(36)}& Fx=C_4{x}^3+C_5{x}^2+C_{6}x+C_7\text{.}\end{array}$$

It can be seen that the intersection of the function with the vertical axis is $A=0,C_7{C}_7>0$. Its derivative is$$\begin{array}{}\text{(37)}& F^{\prime }x=3C_4{x}^2+2C_{5}x+C_6\text{.}\end{array}$$

From this, it follows that the discriminant of the roots of the derivative equation and the axis of symmetry of the image are$$\begin{array}{}\text{(38)}& \Delta =4C_5^2-12C_4{C}_6,X=-\frac{C_5}{3C_4}>0\text{.}\end{array}$$

Since all parameters are non-negative numbers, $C_4$ is less than zero and $C_7$ is large than zero. From the uniform persistence of the system it follows that $r/1+k{Y}_T-d>0$, therefore $C_5$ is large than zero. However, the sign of $C_6$ can be positive or negative.

The next step is to determine the positive roots of the original function based on the derivative.

(A) $C_6$ is a non-negative number

  In this case, $\Delta >0$ and the intersection of $F^{\prime }x$ with the vertical axis is $B=0,C_6{C}_6⩾0$. It follows that $Fx$ monotone increasing and then decreasing on $0,+\infty$. $Fx$ intersects with the positive half axis of the longitudinal axis, so $Fx$ has one positive root. This solution is the horizontal coordinate of the pseudoequilibrium, denoted by $x_1^p$. So the coordinates of the pseudoequilibrium is $E_1^p=x_1^p,Y_T$. See Figure 1(a).

(B) $C_6$ is a negative number

  In this case, the sign of $\Delta$ can be positive or negative. Therefore, we will discuss these categories.

(I) $\Delta \le 0$

  In this case, since $C_4<0$ and $\Delta \le 0$, $F^{\prime }x\le 0$. It follows that $Fx$ is monotone decreasing on $0,+\infty$ and because $Fx$ intersects the positive semiaxis of the vertical axis, $Fx$ has one positive root. This solution is the horizontal coordinate of the pseudoequilibrium, denoted by $x_2^p$. So the coordinate of the pseudoequilibrium is $E_2^p=x_2^p,Y_T$. See Figure 1(b).

(II) $\Delta >0$

  Combining all the above conditions shows that $F^{\prime }x$ has two distinct positive roots. Suppose the two roots are $x_{-}$ and $x_{+}{x}_{-}<x_{+}$. Therefore, $Fx$ is monotone decreasing, then monotone increasing, and finally monotone decreasing on $0,+\infty$. There are five possible scenarios for the solution of $Fx$ as follows:

(i) When $F{x}_{-}>0$, $Fx$ has one positive root. Therefore, the solution is the horizontal coordinate of the pseudoequilibrium, denoted by $x_3^p$. So the coordinate of the pseudoequilibrium is $E_3^p=x_3^p,Y_T$. See Figure 1(c).

(ii) When $F{x}_{-}=0$, $Fx$ has two positive roots, one of which is $x_{-}$ and the other is set to $x_5^p$. Hence the coordinates of the pseudoequilibria are $E_4^p=x_{-},Y_T$ and $E_5^p=x_5^p,Y_T$. Beyond that, we know that $x_{-}<x_5$. See Figure 1(d).

(iii) When $F{x}_{-}<0$ and $F{x}_{+}>0$, $Fx$ has three positive roots. Let these three positive roots be $x_6^p$, $x_7^p$ and $x_8^p$. Therefore, the coordinates of the pseudoequilibria are $E_6^p=x_6^p,Y_T$, $E_7^p=x_7^p,Y_T$ and $E_8^p=x_8^p,Y_T$. Furthermore, we assume that $x_6<x_7<x_8$. See Figure 1(e).

(iv) If $F{x}_{-}<0$ and $F{x}_{+}=0$, then $Fx$ has two positive roots. One of the roots is $x_{+}$ and the other root is set to $x_9^p$. In addition, we know that $x_9<x_{+}$. Therefore, the coordinates of the pseudoequilibria are $E_9^p=x_9^p,Y_T$ and $E_{10}^p=x_{+},Y_T$. See Figure 1(f).

(v) If $F{x}_{-}<0$ and $F{x}_{+}<0$, then $Fx$ has one positive root. Therefore, the solution is the horizontal coordinate of the pseudoequilibrium, denoted by $x_{11}^p$. So the coordinate of the pseudoequilibrium is $E_{11}^p=x_{11}^p,Y_T$. See Figure 1(g).

[figure(s) omitted; refer to PDF]

Next, we analyze the stability of the pseudoequilibria. We can see from Figure 1 that $F^{\prime }x<0$ when $x=x_i^{p}i=\mathrm{1,2,3,5,6,8,9,11}$, so the pseudoequilibrium $E_i^p$ is stable. $F^{\prime }x>0$ when $x=x_7^p$, so the pseudoequilibrium $E_7^p$ is unstable. By using the proof of Theorem 5 in [28] it can be shown that $E_{10}^p$ is a stable pseudoequilibrium and $E_4^p$ is an unstable pseudoequilibrium.

We summarize the above in Table 1.

Table 1

Summary of pseudoequilibrium.

Condition	The case of pseudoequilibrium
$C_6\ge 0$	$E_1^p$ (stable)
$C_6<0$ and $\Delta \le 0$	$E_2^p$ (stable)
$C_6<0$, $\Delta >0$, and $F{x}_{-}>0$	$E_3^p$ (stable)
$C_6<0$, $\Delta >0$, and $F{x}_{-}=0$	$E_4^p$ (unstable) and $E_5^p$ (stable)
$C_6<0$, $\Delta >0,F{x}_{-}<0$, and $F{x}_{+}>0$	$E_6^p$ (stable), $E_7^p$ (unstable), and $E_8^p$ (stable)
$C_6<0$, $\Delta >0,F{x}_{-}<0$, and $F{x}_{+}=0$	$E_9^p$ (stable) and $E_{10}^p$ (unstable)
$C_6<0$, $\Delta >0,F{x}_{-}<0$, and $F{x}_{+}<0$	$E_{11}^p$ (stable)

4.3. Regular Equilibria

For the subsystem $S_1$, $E_{01}=\mathrm{0,0}$, and $E_{11}=r-d/a,0$ are regular equilibria. Furthermore, for the equilibrium $E_1^{\ast }=x_1^{\ast },y_1^{\ast }$ there are the following conclusions. Before discussing whether $E_1^{\ast }$ is a regular or virtual equilibrium, it is necessary to ensure that it exists, in other words, to ensure that the condition $r>d+am/up$ holds. If $y_1^{\ast }<Y_T$, then the equilibrium $E_1^{\ast }$ is a regular equilibrium for the subsystem $S_1$, denoted by $E_1^R$. If $y_1^{\ast }>Y_T$, then the equilibrium $E_1^{\ast }$ is a virtual equilibrium, denoted by $E_1^V$.

Moreover, about regular equilibria for the subsystem $S_2$, $E_{02}=\mathrm{0,0}$, and $E_{12}=r-d/a,0$ are virtual equilibria. In addition, we have the following results based on the size of $Y_T$ and $y_2^{\ast }$. If $y_2^{\ast }>Y_T$ and $r>d+am/up-mq$ with $up>mq$, then the equilibrium $E_2^{\ast }$ is a regular equilibrium for the subsystem $S_2$, denoted by $E_2^R$. If $y_2^{\ast }>Y_T$ and $r>d+am/up-mq$ with $up>mq$, then the equilibrium $E_2^{\ast }$ is a virtual equilibrium for the subsystem $S_2$, denoted by $E_2^V$.

The expressions of $y_1^{\ast }$ and $y_2^{\ast }$ are complex, so their sizes cannot be directly compared. First, we assume that $y_1^{\ast }>y_2^{\ast }$. Then, we obtain an inequality between the parameter $k$ and the other parameters$$\begin{array}{}\text{(39)}& k>\frac{u^3{p}^{4}r-dup-mq-ma-u{p}^2{up-mq}^2}{{up-mq}^2{upr-d-ma}^2+u^2{p}^{2}am+dup-mq}=k^{\ast }\text{.}\end{array}$$

Similarly, when assuming that $y_1^{\ast }<y_2^{\ast }$ it is obtained that $k<k^{\ast }$.

Based on the above analysis, when both $E_1^{\ast }$ and $E_2^{\ast }$ exist and $k>k^{\ast }$, if $y_2^{\ast }<Y_T<y_1^{\ast }$, then the equilibria $E_1^{\ast }$ and $E_2^{\ast }$ are virtual equilibria. However, when both $E_1^{\ast }$ and $E_2^{\ast }$ exist and $k<k^{\ast }$, if $y_1^{\ast }<Y_T<y_2^{\ast }$, then both $E_1^{\ast }$ and $E_2^{\ast }$ are regular equilibria.

We summarize the regular and virtual states of the internal equilibria under different conditions in the following Table 2.

Table 2

Summary of regular equilibria.

Prerequisite	Size relation between $Y_T$ and $y_i^{\ast }i=\mathrm{1,2}$	Regular and virtual state of equilibria
$k>k^{\ast }$	$y_2^{\ast }>Y_T$	$E_1^V$, $E_2^R$
	$y_2^{\ast }<Y_T<y_1^{\ast }$	$E_1^V$, $E_2^V$
	$y_1^{\ast }<Y_T$	$E_1^R$, $E_2^V$
$k<k^{\ast }$	$y_1^{\ast }>Y_T$	$E_1^V$, $E_2^R$
	$y_1^{\ast }<Y_T<y_2^{\ast }$	$E_1^R$, $E_2^R$
	$y_2^{\ast }<Y_T$	$E_1^R$, $E_2^V$

4.4. Boundary Equilibrium

The boundary equilibrium of the Filippov system (8) satisfies the following equation:$$\begin{array}{}\text{(40)}& \begin{array}{l}\frac{rx}{1+\epsilon ky}-dx-a{x}^2-\frac{pxy}{1+\epsilon qx}=0\text{,}\\ \frac{upxy}{1+\epsilon qx}-my=0\text{,}\\ y=Y_T\text{.}\end{array}\end{array}$$

We denote $E_1^B$ provided (40) as $\epsilon =0$, i.e., $Y_T=y_1^{\ast }$ (if $y_1^{\ast }$ exists) and $E_2^B$ provided (40) as $\epsilon =1$, i.e., $Y_T=y_2^{\ast }$ (if $y_2^{\ast }$ exists).

So, the coordinates of the boundary equilibria are given in the following equation:$$\begin{array}{c}\text{(41)}& E_1^B=\frac{m}{up},Y_T,Y_T=y_1^{\ast }=\frac{upr-d-ma}{u{p}^2}\text{,}{E}_2^B=\frac{m}{up-mq},Y_T,Y_T=y_2^{\ast }=\frac{-C_2+\sqrt{C_2^2-4C_1{C}_3}}{2C_1}\text{.}\end{array}$$

Therefore, the boundary equilibria $E_1^B$ and $E_2^B$ exist at the critical values $y_1^{\ast }$ and $y_2^{\ast }$ of $Y_T$, respectively.

4.5. Tangent Points

According to Definition 4, a point $E_T=x_T,Y_T$ will be a tangent point on sliding segment ${\mathrm{\Sigma }}_s$ if$$\begin{array}{}\text{(42)}& F_{S_1}HZ=0\mathrm{o}\mathrm{r}{F}_{S_2}HZ=0\text{.}\end{array}$$Accordingly, the following equations are obtained:$$\begin{array}{}\text{(43)}& upxy-my=0\text{,}\end{array}$$or$$\begin{array}{}\text{(44)}& \frac{upxy}{1+qx}-my=0\text{.}\end{array}$$

The coordinates of the two tangent points can be obtained by solving the above equations (43) and (44): $E_1^T=m/up,Y_T$ and $E_2^T=m/up-mq,Y_T$. When $Y_T$ is the critical value $y_1^{\ast }$$y_2^{\ast }$, the boundary equilibria $E_1^B$$E_2^B$ collide with the tangent points $E_1^T$$E_2^T$.

5. Sliding Bifurcation Analysis of the Filippov System (8)

5.1. Regular/Virtual Equilibrium Bifurcation

Based on the above analysis of the different equilibria it is known that $r$ and $Y_T$ are the key factors that determine the existence of different types of equilibria. Therefore, we define four curves with respect to the parameters $r$ and $Y_T$ as follows:$$\begin{array}{c}\text{(45)}& L_1=r,Y_{T}r=d+\frac{am}{up}\text{,}{L}_2=r,Y_{T}r=d+\frac{am}{up-mq}\text{,}\\ L_3=r,Y_T{Y}_T=\frac{upr-d-ma}{u{p}^2}\text{,}\\ L_4=r,Y_T{Y}_T=\frac{-C_2+\sqrt{C_2^2-4C_1{C}_3}}{2C_1}\text{,}\end{array}$$where $C_1$, $C_2$, and $C_3$ have the same expressions as in (15).

The four curves divide the parameter space into seven regions and mark whether the equilibrium in each region is regular equilibrium or virtual equilibrium. Similarly, the boundary equilibria $E_1^B$ and $E_2^B$ appear on the corresponding curves $L_3$ and $L_4$. In particular, in Figure 2, not only the two virtual equilibria $E_1^V$ and $E_2^V$ can coexist, but also the two regular equilibria $E_1^R$ and $E_2^R$. Therefore, when other parameters are fixed, the regular and virtual states of the internal equilibria of different thresholds will be different according to the size of the birth rate, and the number threshold of prey’s fear effect on predators will affect the future development trend of the prey population.

[figure(s) omitted; refer to PDF]

5.2. Boundary Equilibrium Bifurcation

Boundary equilibrium bifurcation may occur in the Filippov system (8) once $E_P$, $E_T$, and $E_R$ or $E_T$ and $E_R$ collide at the same time as $Y_T$ passes a critical value, where $E_P$, $E_T$, and $E_R$ represent the pseudoequilibrium, tangent point, and regular equilibria of the system, respectively. In this part, $Y_T$ is chosen as the bifurcation parameter, and all other parameters are fixed as those in Figures 3 and 4.

[figure(s) omitted; refer to PDF]

5.2.1. Boundary Node Bifurcation

A simple calculation gives the critical value of $Y_{T1}$ of 0.2864. When $Y_T=0.1$, the stable node $E_2^R$ exists. Also, the visible tangent point $E_2^T$ exists on the boundary, as shown in Figure 3(a). It can be seen from the Figure 3(a) that the trajectories starting from region $S_2$ and remaining at $S_2$ without colliding with the boundary will converge to $E_2^R$. The trajectories from region $S_2$ will also collide with the sliding region or the crossing region. The trajectories which hit the sliding region will pass through the visible tangent point $E_2^T$ and finally converge to $E_2^R$. At the same time, the trajectories of collision with the crossing region converges to $E_2^R$. The trajectories from region $S_1$ will cross the crossing region and eventually converge to $E_2^R$.

The regular equilibrium $E_2^R$ becomes boundary equilibrium $E_2^B$ which collides with the tangent point $E_2^T$, when $Y_T=0.2864$. All the trajectories hitting sliding region or crossing region will converge to $E_2^B$, see Figure 3(b) for more details.

For $Y_T=0.5$, the regular equilibrium $E_2^R$ will turn into a virtual equilibrium $E_2^V$. At the same time, the visible tangent point $E_2^T$ becomes the invisible tangent point $E_2^T$. In addition, the pseudoequilibrium $E_P$ also exists. In this case, the trajectories will remain on the sliding segment and converge to pseudoequilibrium $E_P$, as shown in Figure 3(c).

5.2.2. Boundary Focus Bifurcation

The critical value $Y_{T2}$ is 0.98. From Figure 4(a), the stable focus $E_1^R$ and the virtual equilibrium $E_2^V$ coexist when $Y_T=1.2$. The visible tangent point $E_1^T$ exist. In this case, the trajectories will converge to the regular equilibrium $E_1^R$.

For $Y_T=0.98$, the stable focus $E_1^R$ becomes boundary equilibrium $E_1^B$ which collides with the tangent point $E_1^T$. In this case, all the trajectories hitting discontinuous boundary will converge to $E_1^B$, see Figure 4(b) for more details.

The boundary equilibrium $E_1^B$ will turn into the virtual equilibrium $E_1^V$, and the visible tangent point $E_1^T$ becomes invisible tangent point $E_1^T$. In this case, the trajectories will converge to the pseudoequilibrium $E_P$, as shown in Figure 4(c).

5.3. Global Sliding Bifurcation

It can be seen from [6] that there may be standard periodic solutions that lie entirely in the region $S_2$ when the subsystem $S_2$ through a Hopf bifurcation. At the same time, as described in [29], the Filippov system (8) may add a new periodic solution that slides in the sliding segment, that is, the sliding periodic solutions. The orbits corresponding to the sliding cycle solutions are called sliding cycles. In this section, we focus on grazing bifurcation (i.e., touching bifurcation).

It is known from [29] that standard periodic solutions can collide with sliding segments, and this bifurcation is called a grazing (touching) bifurcation. When $Y_T=6$, the region $S_2$ of the Filippov system (8) will have a stable periodic solution, see Figure 5(a) for more details. At this point, there are two tangent points $E_1^T$ and $E_2^T$ lying on the sliding boundary. The subsystem $S_2$ has an unstable regular equilibrium $E_2^R$ and subsystem $S_1$ has an unstable virtual equilibrium $E_1^V$. It can be seen from Figure 5(b) that when the threshold $Y_T$ increases to about 6.843, a touching bifurcation occurs, which means that the standard periodic solution of the Filippov system (8) collides with its tangent point $E_2^T$. As $Y_T$ continues to increase, the cycle becomes a sliding cycle where sliding occurs on the sliding segment, as shown in Figure 5(c) with $Y_T=6.9$. When $Y_T$ exceeds 6.98, the virtual equilibrium $E_1^V$ becomes the regular equilibrium $E_1^R$. In other words, the Filippov system (8) has an unstable regular equilibrium $E_2^R$ and a stable regular equilibrium $E_1^R$, as shown in Figure 5(d) with $Y_T=7.2$. At this moment, the sliding cycle remains. In addition, some solution trajectories converge to the stable regular equilibrium $E_1^R$, as shown by the black curves in Figure 5(d).

[figure(s) omitted; refer to PDF]

Especially when the bifurcation parameter $Y_T$ increases to 8.8, the stable periodic cycle disappears and the regular equilibrium $E_2^R$ becomes the virtual equilibrium $E_2^V$, as shown in Figure 5(e). Meanwhile, Figure 5(e) also shows that all solution trajectories converge to $E_1^R$. Figure 5(f) takes the sliding cycle out separately to observe the changes in the sliding cycles. It can be seen from Figure 5(f) that the length of sliding on the slide segment decreases as the sliding bifurcation parameter $Y_T$ increases.

5.4. Impact of $k$ on the System

In this section, we analyze the effect of the fear effect $k$ on the dynamics of the subsystem $S_2$ and the Filippov system (8). By fixing parameters and initial conditions other than $k$, we plot the time series diagram of $x$ and $y$ for the subsystem $S_2$ and the Filippov system (8) for $k=0,k=0.3$, and $k=20$.

It can be seen from Figures 6(a) and 7(a) that the number of predator and prey in the two systems is in a state of periodic fluctuation (i.e. the solution is a periodic solution), but the peak value of fluctuation in Figure 6(a) is higher than that in Figure 7(a). The low-level fear effect ultimately reduce the range of fluctuations in the solutions of the subsystem $S_2$ and Filippov system (8), and the number of prey is higher than the number of predator, as detailed in Figures 6(b) and 7(b). By observing Figures 6(c) and 7(c), we can see that the quantitative relationship of subsystem $S_2$ under the high level of fear effect is consistent with that low-level fear effect but finally fluctuates in a smaller range. The high level of fear effect stabilizes the number of predator and prey in the Filippov system (8), but in contrast to the number relationship in the subsystem $S_2$.

[figure(s) omitted; refer to PDF]

6. Discussion

In recent years, the Filippov system has been widely used in integrated pest management, epidemic control, and predator-prey research [15–18, 21, 30]. The Filippov system demonstrates complex dynamic behaviors on sliding segments based on the inheritance of the dynamic behaviors of atomic systems. New equilibria such as pseudoequilibrium, boundary equilibrium and tangent points can appear on the sliding segment. In addition, new bifurcations such as boundary bifurcations and global sliding bifurcations will also appear on the sliding segment.

In this study, we establish and analyze a Filippov predator-prey model with the fear effect. The threshold of the model is the number of predators, and when the number of predators is below the threshold, the prey does not have a fear effect; when the number of predators exceeds the threshold, the prey has fear effect. Also, the functional responses that respond to the predator-prey relationship are different in the two subsystems.

First, we analyze the existence and stability of the equilibrium of the two subsystems, and discuss whether the equilibria of the two subsystems are regular equilibria or virtual equilibria. The analysis results are summarized in Table 2. In addition, the numerical simulation of the regular/virtual equilibrium bifurcation is also carried out. It is found that with the change of threshold and birth rate, the regular and virtuality of the equilibria also change, and there will be two regular equilibria or virtual equilibria, as shown in Figure 2. Second, by using the theory of the Filippov system, we give the conditions for the existence of sliding segments and various equilibria. We focus on the existence and stability of the pseudoequilibrium and summarize the results in Table 1.

Finally, we study the boundary equilibrium bifurcation and global sliding bifurcation. From the above numerical simulation, we find that in the case of $y_1^{\ast }>y_2^{\ast }$, two virtual equilibria coexist in the system, and a stable pseudoequilibrium appears on the sliding segment. In the case of $y_1^{\ast }<y_2^{\ast }$, both equilibria of the system are regular equilibria. With the appropriate parameters, the system has a standard periodic solution, as shown in Figure 5(a). With the increase of $Y_T$, the system exhibits touching bifurcation (Figure 5(b)), sliding cycle (Figure 5(c)), sliding cycle, stable regular equilibrium $E_1^R$ coexistence (Figure 5(d)), and periodic solution disappears (Figure 5(e)).

This paper has obtained some meaningful results, of course, there are some shortcomings. For example, we only analyze the possible existence of pseudoequilibrium using a combination of theoretical and graphical methods and discuss the local stability of the pseudoequilibrium, without examining their global stability, which will provide a direction for our next study.

Acknowledgments

This study was partially supported by the National Natural Science Foundation of China (Grant nos. 11961024, 12261104, and 12201196).

References

[1] A. J. Lotka, Elements of Physical Biology, 1925.

[2] V. Volterra, "Fluctuations in the abundance of a species considered mathematically," Nature, vol. 118 no. 2972, pp. 558-560, DOI: 10.1038/118558a0, 1926.

[3] S. L. Lima, "Nonlethal effects in the ecology of predator-prey interactions," BioScience, vol. 48 no. 1, pp. 25-34, DOI: 10.2307/1313225, 1998.

[4] S. Creel, D. Christianson, "Relationships between direct predation and risk effects," Trends in Ecology & Evolution, vol. 23 no. 4, pp. 194-201, DOI: 10.1016/j.tree.2007.12.004, 2008.

[5] W. Cresswell, "Predation in bird populations," Journal of Ornithology, vol. 152 no. S1, pp. 251-263, DOI: 10.1007/s10336-010-0638-1, 2011.

[6] X. Wang, L. Zanette, X. Zou, "Modelling the fear effect in predator–prey interactions," Journal of Mathematical Biology, vol. 73 no. 5, pp. 1179-1204, DOI: 10.1007/s00285-016-0989-1, 2016.

[7] J. Wang, Y. H. Chen, L. Lv, J. Li, "Overlapping group screening for binary cancer classification with TCGA high-dimensional genomic data," Journal of Bioinformatics and Computational Biology, vol. 21 no. 3,DOI: 10.1142/S0219720023500130, 2023.

[8] X. Zhang, H. Zhu, Q. An, "Dynamics analysis of a diffusive predator-prey model with spatial memory and nonlocal fear effect," Journal of Mathematical Analysis and Applications, vol. 525 no. 1,DOI: 10.1016/j.jmaa.2023.127123, 2023.

[9] M. Chen, Y. Takeuchi, J. Zhang, "Dynamic complexity of a modified Leslie–Gower predator–prey system with fear effect," Communications in Nonlinear Science and Numerical Simulation, vol. 119,DOI: 10.1016/j.cnsns.2023.107109, 2023.

[10] S. K. Sasmal, "Population dynamics with multiple allee effects induced by fear factors–a mathematical study on prey-predator interactions," Applied Mathematical Modelling, vol. 64,DOI: 10.1016/j.apm.2018.07.021, 2018.

[11] H. Zhang, Y. Cai, S. Fu, W. Wang, "Impact of the fear effect in a prey-predator model incorporating a prey refuge," Applied Mathematics and Computation, vol. 356, pp. 328-337, DOI: 10.1016/j.amc.2019.03.034, 2019.

[12] S. K. Sasmal, Y. Takeuchi, "Dynamics of a predator-prey system with fear and group defense," Journal of Mathematical Analysis and Applications, vol. 481 no. 1,DOI: 10.1016/j.jmaa.2019.123471, 2020.

[13] B. Chakraborty, H. Baek, N. Bairagi, "Diffusion-induced regular and chaotic patterns in a ratio-dependent predator–prey model with fear factor and prey refuge, Chaos," An Interdisciplinary Journal of Nonlinear Science, vol. 31 no. 3, 2021.

[14] R. Yang, D. Jin, "Dynamics in a predator-prey model with memory effect in predator and fear effect in prey," Electronic Research Archive, vol. 30 no. 4, pp. 1322-1339, DOI: 10.3934/era.2022069, 2022.

[15] S. Tang, J. Liang, Y. Xiao, R. A. Cheke, "Sliding bifurcations of Filippov two stage pest control models with economic thresholds," SIAM Journal on Applied Mathematics, vol. 72 no. 4, pp. 1061-1080, DOI: 10.1137/110847020, 2012.

[16] W. Qin, X. Tan, M. Tosato, X. Liu, "Threshold control strategy for a non-smooth Filippov ecosystem with group defense," Applied Mathematics and Computation, vol. 362,DOI: 10.1016/j.amc.2019.06.046, 2019.

[17] W. Li, Y. Chen, L. Huang, J. Wang, "Global dynamics of a Filippov predator-prey model with two thresholds for integrated pest management," Chaos, Solitons & Fractals, vol. 157,DOI: 10.1016/j.chaos.2022.111881, 2022.

[18] W. Li, L. Huang, J. Wang, "Global asymptotical stability and sliding bifurcation analysis of a general Filippov-type predator-prey model with a refuge," Applied Mathematics and Computation, vol. 405,DOI: 10.1016/j.amc.2021.126263, 2021.

[19] J. Wang, F. Zhang, L. Wang, "Equilibrium, pseudoequilibrium and sliding-mode heteroclinic orbit in a Filippov-type plant disease model," Nonlinear Analysis: Real World Applications, vol. 31, pp. 308-324, DOI: 10.1016/j.nonrwa.2016.01.017, 2016.

[20] C. Chen, Y. Kang, R. Smith, "Sliding motion and global dynamics of a Filippov fire-blight model with economic thresholds," Nonlinear Analysis: Real World Applications, vol. 39, pp. 492-519, DOI: 10.1016/j.nonrwa.2017.08.002, 2018.

[21] S. Tang, Y. Tian, R. A. Cheke, "Dynamic complexity of a predator-prey model for ipm with nonlinear impulsive control incorporating a regulatory factor for predator releases," Mathematical Modelling and Analysis, vol. 24 no. 1, pp. 134-154, DOI: 10.3846/mma.2019.010, 2019.

[22] M. Biák, T. Hanus, D. Janovská, "Some applications of Filippov’s dynamical systems," Journal of Computational and Applied Mathematics, vol. 254, pp. 132-143, DOI: 10.1016/j.cam.2013.03.034, 2013.

[23] Y. V. Tyutyunov, L. Titova, "From lotka–volterra to arditi–ginzburg: 90 years of evolving trophic functions," Biology Bulletin Reviews, vol. 10 no. 3, pp. 167-185, DOI: 10.1134/s207908642003007x, 2020.

[24] J. Yang, S. Tang, R. A. Cheke, "Global stability and sliding bifurcations of a non-smooth Gause predator–prey system," Applied Mathematics and Computation, vol. 224,DOI: 10.1016/j.amc.2013.08.024, 2013.

[25] X. Zhang, S. Tang, "Existence of multiple sliding segments and bifurcation analysis of Filippov prey–predator model," Applied Mathematics and Computation, vol. 239, pp. 265-284, DOI: 10.1016/j.amc.2014.04.098, 2014.

[26] X. Zhang, S. Tang, "Filippov ratio-dependent prey-predator model with threshold policy control," Abstract and Applied Analysis, vol. 2013,DOI: 10.1155/2013/280945, 2013.

[27] R. Cristiano, T. Carvalho, D. J. Tonon, D. J. Pagano, "Hopf and homoclinic bifurcations on the sliding vector field of switching systems in R3: a case study in power electronics," Physica D: Nonlinear Phenomena, vol. 347, pp. 12-20, DOI: 10.1016/j.physd.2017.02.005, 2017.

[28] B. Tang, W. Zhao, "Sliding dynamics and bifurcations of a Filippov system with nonlinear threshold control," International Journal of Bifurcation and Chaos, vol. 31 no. 14,DOI: 10.1142/s021812742150214x, 2021.

[29] Y. A. Kuznetsov, S. Rinaldi, A. Gragnani, "One-parameter bifurcations in planar Filippov systems," International Journal of Bifurcation and Chaos, vol. 13 no. 08, pp. 2157-2188, DOI: 10.1142/s0218127403007874, 2003.

[30] A. Wang, Y. Xiao, "Sliding bifurcation and global dynamics of a Filippov epidemic model with vaccination," International Journal of Bifurcation and Chaos, vol. 23 no. 8,DOI: 10.1142/s0218127413501447, 2013.