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.
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.,where and 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. and are the natural birth rate and natural mortality rate of the prey, respectively, indicates the decay rate of the prey due to intraspecific competition, and denotes the rate of search for the prey by the predator. The rate of conversion of prey biomass to predator biomass is denoted by , and expresses the natural mortality of the predator. In addition, the assumption that 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 . That is, there is no fear effect before the predator reaches the threshold value , and the interactions between the predator and prey is assumed to follow model (1). When , 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]where is the half-saturation constant and is a decreasing function with respect to and . reflects the degree of fear that drives antipredatory behaviours in prey. According to the biological significance of , , and , the function should satisfy the following four assumptions:(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): if there are no predators, there is no fear effect on prey, and therefore the natural birth rate of prey does not decrease(3): the natural birth rate of prey decreases with the increasing fear effect(4): if the density of predators increases, the fear effect increases, and therefore the natural birth rate of prey decreases
A simple form of with the above assumptions satisfied is , proposed by [6]. So, system (2) can be rewritten as follows:
So, systems (2) and (3) can be integrated into the following nonsmooth dynamical system:with
2.2. Preliminaries
To further analyze the Filippov system, here are some symbols [24–26]. Let be the smooth vector fields and be a partition function that depends on . Let 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:
For convenience, we denote
Then, system (4) with (6) can be rewritten as the following Filippov system:where and . be called as the switching manifold, which is the separation boundary of the region and the discontinuity boundary set can be defined as
Furthermore, divides into two distinct regions and . Then, the Filippov system (8) in region is called subsystem , .
Letwhere denotes the standard scalar product, and Lie derivative is used to denote the directional derivative of concerning the vector field at . Therefore, we partition the discontinuity boundary according to the sign case of .
Accordingly, the boundary will be classified as follows:(i) is called the escaping region if on (ii) is called the sliding region if on (iii) is called the crossing region if on
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 makes or , the point is called the regular equilibrium of the Filippov system (8). If a point makes or , the point is called the virtual equilibrium of the Filippov system (8).
Definition 2. If a point is the equilibrium of the sliding region of the Filippov system (8), the point is called pseudoequilibrium and , that is, , and , where
Definition 3. If or , the point is called a boundary equilibrium of the Filippov system (8).
Definition 4. If , , and or , the point 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
When , the dynamical behaviours of system (8) is determined by the subsystems defined by (1). The subsystem has the following three equilibria: , , and where . The equilibria and always exist. The point is in the first quadrant if . 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 is a saddle(II)The boundary equilibrium is locally asymptotically stable if(III)The interior equilibrium 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 , , and are given bySince and , the trivial equilibrium is a saddle.
The trace and determinant corresponding to the matrix are and . Both and can be positive or negative, therefore, from the Routh–Hurwitz criterion, it follows that is locally asymptotically stable when and . Thus, the boundary equilibrium is locally asymptotically stable when . Similarly, the condition of local asymptotic stability of internal equilibrium is .
Theorem 6. When , the boundary equilibrium is globally asymptotically stable. When , the interior equilibrium is globally asymptotically stable.
Proof. The two functions on the right-hand side of system (1) are denoted by and . We consider the Dulac function as . After calculations, we obtainFor . Thus, according to the Dulac–Bendixson theorem, system (1) has no periodic orbits in . Besides, when , the positive equilibrium in has only the interior equilibrium . So all positive solutions will tend to . Combining the proved in Theorem 5 is locally asymptotically stable, it can be concluded that is globally asymptotically stable if condition (12) does not hold.
When , there are only two equilibria and in . In addition, there is no periodic orbit in which means that every positive solution converges to or . It is easy to obtain that is repelling when , so every positive solution converges to . Combining the proved in Theorem 5 is locally asymptotically stable, it can be concluded that is globally asymptotically stable if .
3.2. Dynamics of the Subsystem
For the subsystem , it is easy to calculate the following three equilibria: and where with
Notice that the equilibria and are always present and the equilibrium is feasible if with .
Theorem 7. The trivial equilibrium is a saddle.
Proof. The Jacobian matrix at isObviously, we can get and , so by the Routh–Hurwitz criterion it follows that the trivial equilibrium is a saddle.
The following lemmas [6] give the stability analysis regarding the boundary equilibrium and the internal equilibrium .
Lemma 8. The boundary equilibrium is locally asymptotically stable ifis satisfied and is unstable ifholds.
Lemma 9. The boundary equilibrium is globally asymptotically stable if holds.
Lemma 10. The internal equilibrium is locally asymptotically stable iforIt is unstable if
Lemma 11. The internal equilibrium is globally asymptotically stable if
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:Let the right-hand sides of the above two functions be zero to get and . By definitions, we can obtain the following theorem for the sliding region, escaping region, and crossing region.
Theorem 12. (I)If , the sliding and crossing segment of the Filippov system (8) can be defined as(II)If , the escaping and crossing segment of the Filippov system (8) can be defined as
Proof. (I)If , we can know that . According to (23), (24), and the definition of sliding region, we obtain the sliding region Moreover, there are two possible scenarios that can lead to the presence of transversal intersection and crossing region: or From (28), we calculate that and , i.e., . So, we get the crossing region . Also, from (29), we can get and , i.e., . Thus, the crossing region is . In summary, the crossing region is(II)If , we can know that . By definition, we can get . Thus, the escaping region isSimilarly, (28) and (29) can be obtained for the crossing region. Simplifying (28) can yield and , i.e., . To ensure that the solution must satisfy biological significance, must be a non-negative constant, so this situation does not exist. For (29), the calculation gives and , i.e., . So the crossing region is
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 calculateSo,
Substituting the expressions for and into , we can obtain the following equation:where , , , and .
We know from [27] that we only need to discuss the case of the molecular roots of .
Let
It can be seen that the intersection of the function with the vertical axis is . Its derivative is
From this, it follows that the discriminant of the roots of the derivative equation and the axis of symmetry of the image are
Since all parameters are non-negative numbers, is less than zero and is large than zero. From the uniform persistence of the system it follows that , therefore is large than zero. However, the sign of can be positive or negative.
The next step is to determine the positive roots of the original function based on the derivative.(A) is a non-negative number In this case, and the intersection of with the vertical axis is . It follows that monotone increasing and then decreasing on . intersects with the positive half axis of the longitudinal axis, so has one positive root. This solution is the horizontal coordinate of the pseudoequilibrium, denoted by . So the coordinates of the pseudoequilibrium is . See Figure 1(a).(B) is a negative number In this case, the sign of can be positive or negative. Therefore, we will discuss these categories.(I) In this case, since and , . It follows that is monotone decreasing on and because intersects the positive semiaxis of the vertical axis, has one positive root. This solution is the horizontal coordinate of the pseudoequilibrium, denoted by . So the coordinate of the pseudoequilibrium is . See Figure 1(b).(II) Combining all the above conditions shows that has two distinct positive roots. Suppose the two roots are and . Therefore, is monotone decreasing, then monotone increasing, and finally monotone decreasing on . There are five possible scenarios for the solution of as follows:(i)When , has one positive root. Therefore, the solution is the horizontal coordinate of the pseudoequilibrium, denoted by . So the coordinate of the pseudoequilibrium is . See Figure 1(c).(ii)When , has two positive roots, one of which is and the other is set to . Hence the coordinates of the pseudoequilibria are and . Beyond that, we know that . See Figure 1(d).(iii)When and , has three positive roots. Let these three positive roots be , and . Therefore, the coordinates of the pseudoequilibria are , and . Furthermore, we assume that . See Figure 1(e).(iv)If and , then has two positive roots. One of the roots is and the other root is set to . In addition, we know that . Therefore, the coordinates of the pseudoequilibria are and . See Figure 1(f).(v)If and , then has one positive root. Therefore, the solution is the horizontal coordinate of the pseudoequilibrium, denoted by . So the coordinate of the pseudoequilibrium is . See Figure 1(g).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
Next, we analyze the stability of the pseudoequilibria. We can see from Figure 1 that when , so the pseudoequilibrium is stable. when , so the pseudoequilibrium is unstable. By using the proof of Theorem 5 in [28] it can be shown that is a stable pseudoequilibrium and is an unstable pseudoequilibrium.
We summarize the above in Table 1.
4.3. Regular Equilibria
For the subsystem , , and are regular equilibria. Furthermore, for the equilibrium there are the following conclusions. Before discussing whether is a regular or virtual equilibrium, it is necessary to ensure that it exists, in other words, to ensure that the condition holds. If , then the equilibrium is a regular equilibrium for the subsystem , denoted by . If , then the equilibrium is a virtual equilibrium, denoted by .
Moreover, about regular equilibria for the subsystem , , and are virtual equilibria. In addition, we have the following results based on the size of and . If and with , then the equilibrium is a regular equilibrium for the subsystem , denoted by . If and with , then the equilibrium is a virtual equilibrium for the subsystem , denoted by .
The expressions of and are complex, so their sizes cannot be directly compared. First, we assume that . Then, we obtain an inequality between the parameter and the other parameters
Similarly, when assuming that it is obtained that .
Based on the above analysis, when both and exist and , if , then the equilibria and are virtual equilibria. However, when both and exist and , if , then both and are regular equilibria.
We summarize the regular and virtual states of the internal equilibria under different conditions in the following Table 2.
4.4. Boundary Equilibrium
The boundary equilibrium of the Filippov system (8) satisfies the following equation:
We denote provided (40) as , i.e., (if exists) and provided (40) as , i.e., (if exists).
So, the coordinates of the boundary equilibria are given in the following equation:
Therefore, the boundary equilibria and exist at the critical values and of , respectively.
4.5. Tangent Points
According to Definition 4, a point will be a tangent point on sliding segment ifAccordingly, the following equations are obtained:or
The coordinates of the two tangent points can be obtained by solving the above equations (43) and (44): and . When is the critical value , the boundary equilibria collide with the tangent points .
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 and are the key factors that determine the existence of different types of equilibria. Therefore, we define four curves with respect to the parameters and as follows:where , , and 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 and appear on the corresponding curves and . In particular, in Figure 2, not only the two virtual equilibria and can coexist, but also the two regular equilibria and . 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.
5.2. Boundary Equilibrium Bifurcation
Boundary equilibrium bifurcation may occur in the Filippov system (8) once , , and or and collide at the same time as passes a critical value, where , , and represent the pseudoequilibrium, tangent point, and regular equilibria of the system, respectively. In this part, is chosen as the bifurcation parameter, and all other parameters are fixed as those in Figures 3 and 4.
(a)
(b)
(c)
(a)
(b)
(c)
5.2.1. Boundary Node Bifurcation
A simple calculation gives the critical value of of 0.2864. When , the stable node exists. Also, the visible tangent point 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 and remaining at without colliding with the boundary will converge to . The trajectories from region 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 and finally converge to . At the same time, the trajectories of collision with the crossing region converges to . The trajectories from region will cross the crossing region and eventually converge to .
The regular equilibrium becomes boundary equilibrium which collides with the tangent point , when . All the trajectories hitting sliding region or crossing region will converge to , see Figure 3(b) for more details.
For , the regular equilibrium will turn into a virtual equilibrium . At the same time, the visible tangent point becomes the invisible tangent point . In addition, the pseudoequilibrium also exists. In this case, the trajectories will remain on the sliding segment and converge to pseudoequilibrium , as shown in Figure 3(c).
5.2.2. Boundary Focus Bifurcation
The critical value is 0.98. From Figure 4(a), the stable focus and the virtual equilibrium coexist when . The visible tangent point exist. In this case, the trajectories will converge to the regular equilibrium .
For , the stable focus becomes boundary equilibrium which collides with the tangent point . In this case, all the trajectories hitting discontinuous boundary will converge to , see Figure 4(b) for more details.
The boundary equilibrium will turn into the virtual equilibrium , and the visible tangent point becomes invisible tangent point . In this case, the trajectories will converge to the pseudoequilibrium , 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 when the subsystem 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 , the region 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 and lying on the sliding boundary. The subsystem has an unstable regular equilibrium and subsystem has an unstable virtual equilibrium . It can be seen from Figure 5(b) that when the threshold 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 . As continues to increase, the cycle becomes a sliding cycle where sliding occurs on the sliding segment, as shown in Figure 5(c) with . When exceeds 6.98, the virtual equilibrium becomes the regular equilibrium . In other words, the Filippov system (8) has an unstable regular equilibrium and a stable regular equilibrium , as shown in Figure 5(d) with . At this moment, the sliding cycle remains. In addition, some solution trajectories converge to the stable regular equilibrium , as shown by the black curves in Figure 5(d).
(a)
(b)
(c)
(d)
(e)
(f)
Especially when the bifurcation parameter increases to 8.8, the stable periodic cycle disappears and the regular equilibrium becomes the virtual equilibrium , as shown in Figure 5(e). Meanwhile, Figure 5(e) also shows that all solution trajectories converge to . 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 increases.
5.4. Impact of on the System
In this section, we analyze the effect of the fear effect on the dynamics of the subsystem and the Filippov system (8). By fixing parameters and initial conditions other than , we plot the time series diagram of and for the subsystem and the Filippov system (8) for , and .
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 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 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 .
(a)
(b)
(c)
(a)
(b)
(c)
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 , two virtual equilibria coexist in the system, and a stable pseudoequilibrium appears on the sliding segment. In the case of , 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 , the system exhibits touching bifurcation (Figure 5(b)), sliding cycle (Figure 5(c)), sliding cycle, stable regular equilibrium 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.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was partially supported by the National Natural Science Foundation of China (Grant nos. 11961024, 12261104, and 12201196).