Stability of a fear effect predator–prey model with mutual interference or group defense

In this paper, we consider a fear effect predator–prey model with mutual interference or group defense. For the model with mutual interference, we show the interior equilibrium is globally stable, and the mutual interference can stabilize the predator–prey system. For the model with group defense, we discuss the singular dynamics around the origin and the occurrence of Hopf bifurcation, and find that there is a separatrix curve near the origin such that the orbits above which tend to the origin and the orbits below which tend to limit cycle or the interior equilibrium.


Introduction
Most of predator-prey systems take into account only the direct killing by predators [2,7,11,12,19,20,22,29,30,32,34,35,38,40], but not the cost of fear effect, which is considered to be more powerful than direct predation. Zanette et al. [36] showed that the bird reduce 40% less offspring by predation fears. Therefore, Wang et al. [25] investigated the following predator-prey system with fear effect They considered the cost of fear into prey reproduction. a is the birth rate of prey; d is the natural death rate of prey; f is the level of fear; F(f , y) is the cost of anti-predator defence due to fear.
Here F(f , y) can be reasonably assumed They showed the fear effect has no impact on the stability of system (1). But for the Holling type II of system (1), they showed that the fear can make the system become stable. Cong et al. [5] investigated a three-species food chain model with fear effect and showed that the fear effect can change the model from a chaotic state to a stable state. Based on system (1), Pal et al. [14] discussed the stability, Hopf-bifurcation and Bogdanov-Takens bifurcation of a predator-prey system with fear effect and hunting cooperation. Zhang et al. [37] investigated the influence of fear effect and prey refuge on the stability of a predator-prey system and showed that the fear effect can stabilize the system. Tiwari et al. [23] considered a predator-prey model with Beddington-DeAngelis functional response and the fear effect into prey, and discussed the stability and bifurcations of the system. By incorporating spatial memory and pregnancy period, Wang et al. [27] performed the detailed bifurcation analysis, and their results indicated that the memory ability and pregnancy cycle of prey play significant roles in the spatiotemporal dynamics of prey-predator models in an intimidatory environment. Sasmal [18] studied the following predator-prey model with fear effect and Allee effect where 0 < δ < k expresses a strong Allee effect; 1 1+fv is one type of fear effect, which satisfies (2). Due to the increase of the fear effect, the author showed the decrease of per-capita growth rate and the multistability of the system. From more results about fear effect, see [16,21,26].
In 1975, Hassell's research of the capturing behaviour between larvae of Plodia interpunctella and Colorado potato beetles showed that the predators will leave each other when they encountered. Thus this mutual interference phenomenon can reduce the search efficiency of predators. This functional response is particularly applicable to systems in patchy environments where parasites tend to cluster in some patches rather than others. Hassell [10] introduced a mutual interference constant θ (0 < θ < 1) into a Volterra model. Although models considering predator interference have different mathematical expressions, the same qualitative results are given [3]. That is when predator interference is low, increasing predator interference has a positive effect on the asymptotic stability of system. Further, Freedman [9] analysed the stability of the mutual interference system. Wang and Zhu [28] considered the stability of a delayed impulsive prey-predator system with mutual interference. Xiao and Li [31] considered a predator-prey model (1) with mutual interference where 0 < θ < 1. When θ = 1, system (4) is reduced to model (1). Notice that the positive equilibrium of system (1) is unstable, but when considering system (1) with mutual interference, that is system (4), this unstable positive equilibrium becomes stable. Hence, the mutual interference can stabilize the predator-prey system.
In nature, prey gathers together in herds to protect themselves from predator. That is prey allows the weakest individuals to occupy the inside of the herd for defensive purposes, leaving the healthier and stronger prey around it. Then in the predator-prey model, in which prey exhibits herd behaviour, predator interacts with prey along the outer corridor of the prey group. By assuming the prey individuals at the herd boundary is proportional to the square root of the area, Braza [4] and Ajraldi et . al [1] studied the models with square root of the prey population, which represent the predator interacts with the prey along the outer corridor of the herd of prey. Due to the square root term, Braza [4] showed the origin is more subtle. That is, the first quadrant is divided into two parts by a separatrix curve near the origin. The orbits above the separatrix curve will tend to origin, but the orbits below the separatrix curve will leave the origin. This phenomenon also appears in the predatorprey with ratio-dependent functional responses [7,32]. Salman [17] considered the stability of a discrete-time system with square root functional responses. When considering the interaction between the prey and the predator in both cases 2D and 3D herd shapes, a new functional response in the θ power of prey has been proposed. The authors [24] and [33] proposed a functional response in terms of the θ (0 < θ < 1) power of prey, which is more general than [1,4]. More precisely, they studied the following model: where 0 < θ < 1. When θ = 1 2 , the functional responses of the above system are reduced to square root functional responses [1,4]. The authors [24,33] found that the solution behaviour near the origin is singularity and can occur Hopf bifurcation. Djilali [8] extended the model in [24,33] to functional response with Holling-II.
In this paper, we study the following model: where all parameters of system (5) are positive constants; r is the birth rate of prey; K is the carrying capacity; p is the capture rate; h is the conversion efficiency; e is the death rate of predator; 1 1+fy is the cost of anti-predator defence due to fear [18,25]; f is the level of fear; the fear effect can lead to the decrease of per-capita growth rate. Let Dropping the bar, system (5) is reduced to Based on the above discussion on mutual interference and group defense, we expect to discuss the effect of mutual interference and group defense on the dynamical behaviours of the system (6). Therefore, we also respectively study where θ 1 (0 < θ 1 < 1) is the mutual interference constant; and where θ 2 (0 < θ 2 < 1) is the group defense constant. The organization of this paper is as follows: In Section 2, the stability of system (6) is investigated. In Sections 3 and 4, the dynamic behaviours of systems (7) and (8) are respectively investigated. In Section 5, the influences of fear effect, mutual interference and group defense on population density are respectively discussed. Finally, a brief discussion is given.

Fear effect predator-prey model (6)
The following lemma shows solutions of system (6) are positive and uniformly bounded.
for all large t.
Proof. Obviously, solutions of system (6) are positive. If x(t) ≥ 1, it follows from the first equation of system (6) that x (t) < 0. Then 0 < x(t) < 1 for all large t.
Define a function V(x, y) = ax + y. Differentiating V along the solution of system (6) for large t, we have for all large t. This completes the proof of Lemma 2.1. Now, we consider the existence and stability of the equilibria of system (6). There always exist a trivial equilibrium E 1 0 (0, 0) and a boundary equilibria E 1 1 (1, 0). By computation, system (6) has a positive equilibrium E * Proof: The eigenvalues of E 1 0 are λ 1 = 1 and λ 2 = −d, so E 1 0 is unstable. Note that the eigenvalues at E 1 When a = d, i.e. λ 2 = 0, translating E 1 1 to the origin by letting x 1 = x − 1, y 1 = y, we obtain the Taylor expansions of system (6) as follows: Further, letting x 2 = x 1 + y 1 , y 2 = y 1 , rewriting x 2 , y 2 as x, y again, we have It follows from Theorem 7.1 in Zhang et al. [39] that E 1 1 is a saddle-node, which includes a stable parabolic sector. Then E 1 1 is stable if a ≤ d. This completes the proof of Lemma 2.2.
is locally asymptotically stable and E * 1 does not exist. Assume that system (6) has a closed orbit, then there must exist an equilibrium in the interior of the closed orbit, which is impossible. Therefore, system (6) does not exist a limit cycle. Hence, E 1 1 is globally asymptotically stable. (2) If a > d, E 1 0 and E 1 1 are unstable. The Jacobin matrix at E * 1 is calculated as Noting that x * 1 < 1, the eigenvalues at E * 1 satisfies for all x > 0 and y > 0. It follows from the Dulac theorem that there is no closed orbit.
Hence, E * 1 is globally asymptotically stable. This completes the proof of Theorem 2.1.
It follows from the above discussion that l 1 and l 2 have a unique interior intersection point in the first quadrant. Therefore, system (7) admits a unique positive equilibrium . This completes the proof of Lemma 3.1. (7) is globally asymptotically stable.
The Jacobin matrix at E * 2 is calculated as Noting that x * 2 < 1, the eigenvalues at E * 2 satisfy λ 1 λ 2 > 0 and λ 1 + λ 2 < 0, that is E * 2 is stable. Define a Dulac function B(x, y) = x −1 y −1 , then for all x > 0 and y > 0. Then there is no closed orbit in system (7). Hence, E * 2 is globally asymptotically stable. This completes the proof of Theorem 3.1.

Fear effect predator-prey model with group defense
In this section, we consider system (7). There always exist a trivial equilibrium E 3 0 (0, 0) and a boundary equilibria E 3 1 (1, 0). By simple computation, if a > d, system (8) has a positive equilibrium Stimulated by the works of [4,24,33], we have the following theorem about the singular dynamics of system (8) around the origin. From the first equation of system (8), we have Considering the following Bernoulli equation: Hence, , the prey x is extinct. There exists a positive ε such that aε θ 2 − d < 0. Since x is extinct, there exists a t 1 such that x(t) ≤ ε for t > t 1 . Form the first equation of system (8), we have Similarly, we have y(t) ≤ y(t 1 )e (aε θ 2 −d)(t−t 1 ) . Then lim t→+∞ y(t) = 0.

Remark 4.1:
Since the singularity of the origin, then the origin cannot be linearized. If we use the method in the proof of Theorem 3.1, then the origin is a saddle which is incorrect. If the initial condition of system (8) above the separatrix curve L, the orbits will intersect with y-axis, and tend to the origin along the y-axis. If a = d, then λ 2 = 0. Translating E 3 1 to the origin by X = x−1, Y = y, rewriting X and Y as x and y respectively, we obtain the Taylor expansions of system (8) as follows: Letting X = x + y, Y = y, rewriting X and Y as x and y respectively, system (11) can be written as It follows from Theorem 7.1 in Zhang et al, [39] that E 3 1 is a saddle-node, which includes a stable parabolic sector. Then E 3 1 is stable if a ≤ d. This completes the proof of Theorem 4.2.
Hence, system (8) undergoes a Hopf bifurcation at a = d * . With the increase of a, tr(J(E * 3 )) changes sign from negative to positive, which implies E * 3 loses its stability and Hopf bifurcation occurs. Now we determine the stability of the limit cycle by calculating the first Lyapunov coefficient. Since the group defense is defined by the θ 2 (0 < θ 2 < 1) power of prey and the expression of E * 3 is complicated, it is difficult to calculate the first Lyapunov coefficient. For simplicity, stimulated by the technique of [6,13], we make the following scaling for system (8) Note that a(x * 3 ) θ 2 = d and x * 3 < 1, thenā =d and b > 1, respectively. Dropping the bar, system (8) is reduced to (8) is transformed to E * (1, 1) of system (12). Hence, we have Then system (12) can be written as follows: Obviously system (13) has the same topological structure as that of system (8), then we calculate the first Lyapunov coefficient of system (13) to determine the stability of the limit cycle. The Jacobian matrix of system (13) at E * (1, 1) takes the form Let b * = 1 + 1 1−θ 2 , then we obtain tr(J(E * ))| b=b * = 0. We have the transversality condition as follows: Hence, system (13) undergoes a Hopf bifurcation at b = b * (i.e. a = d * of system (8)).
From the above discussion, we have the following theorem. (1) If d < a < d * , then E * 3 is locally asymptotically stable.
If a = d * , then system (4.1) undergoes a supercritical Hopf bifurcation at E * 3 and exists a stable limit cycle around E * 3 .

The influence of fear effect on population density
Obviously, the fear effect has no impact on the stability of systems (6) -(8) and the prey density at the positive equilibrium of systems (6) and (8). Now we discuss the influence of the fear effect on the predator density at the positive equilibrium of systems (6) and (8). Let By computation, we obtain Then the fear effect will change the density of predator population at the positive equilibrium of systems (6) and (8), which decreases with the increase of fear effect. Next we study the influence of the fear effect on the positive equilibrium of system (7). Denote Similarly to the analysis of the above subsection, we can easily obtain .
We can observe from where J is given in Section 5.1. Obviously, x * 2 increases with the increase of θ 1 (Figure 2  a); meanwhile, y * 2 increases with the increase of θ 1 when x * 2 < 1 2 , and y * 2 decreases with the increase of θ 1 when x * 2 > 1 2 (Figure 2 b). Therefore when the prey density is relatively small, y * 2 increases as θ 1 increases, that is the mutual interference constant θ 1 has less influence on the capturing behaviour of the predator; but for the fact that the increase of θ 1 can result in the increase of the prey density, when the prey density becomes relatively large, the effect of mutual interference on the capturing behaviour becomes stronger, which leads to the decrease of the predator density.

The influence of group defense on population density
Assume that a > d, then system (8) has a positive equilibrium E * 3 (x * 3 , y * 3 ). Next we discuss the influence of group defense on population density under a > d. By calculation, we can easily obtain dx *     . Bifurcation diagrams of system (6). In region D 2 , E * 1 is globally asymptotically stable. In region D 1 , E 1 1 is globally asymptotically stable.   . Bifurcation diagrams of system (8). In region D 2 , E * 1 is unstable. In region D 1 , E * 1 is locally asymptotically stable. In region D 0 , E 1 1 is locally asymptotically stable.
The above analysis shows that the group defense is beneficial to the density of the prey at the positive equilibrium of system (8). For the predator, when the capture rate is large enough which can result in a large enough a, or the group defense constant is relatively small, the group defense has less influence on the capturing behaviour of the predator; but if the capture rate is relatively small, the group defense plays an important role on the capturing behaviour of the predator, the density of the predator at the positive equilibrium decreases with the increase of θ 2 .

Conclusion
In this paper, we considered a fear effect predator-prey model with mutual interference or group defense and studied the stability and Hopf bifurcation of the system. The stability property of system (6) is simple. If a ≤ d, E 1 1 is globally asymptotically stable ( Figure 5  a-b). If a > d, E * 1 is globally asymptotically stable (Figure 5 c). Considering the original coefficients of system (5), a ≤ d is equivalent to K ≤ e hp . Form Figure 6, if a ≤ d, that is the carrying capacity K of prey is small enough, predator will be extinct, while prey will survive. If a > d, that is the carrying capacity K of prey is large enough, predator and prey can coexist. However, the cost of fear effect has no impact on the stability of system (6), which is in accord with the predator-prey model with the linear functional response [25]. However, the mutual interference or group defense can change the stability of system (6).
By introducing the parameter θ 1 (0 < θ 1 < 1), we considered system (6) with mutual interference and proved that the positive equilibrium E * 2 of system (7) is globally stable unconditionally. However we note that both the predator and prey of system (6) are globally stable only when a > d, and when a ≤ d the predator is extinct. Hence the introduction of the mutual interference constant strengthens the stability of the system. Figure 7(a) -(c) show the global asymptotic stability of system (7) when θ 1 = 2 3 under the assumptions a < d, a = d and a > d respectively.
Incorporating system (6) with group defense, the stability property of system (8) is significantly different from that of systems (6) and (7). It follows from Theorem 4.1 that there exists a separatrix curve such that the first quadrant is divided into two parts. In one part, from Theorem 4.1, the orbits will tend to the origin. In the other part, it follows from Theorems 4.2 and 4.3 that the orbits maybe tend to the boundary equilibria E 3 1 (Figures 8 a-b), the positive equilibrium E * 3 (Figure 8 c), or the limit cycle (Figure 8 d). Note that the stability of boundary equilibria or positive equilibrium of systems (6) and (7) is global. However, when considering system (6) with group defense, the stability of boundary equilibria or positive equilibrium of system (8) is just local due to the singularity of the origin. Considering the original coefficients of system (5), a < d * is equivalent to K ≤ e hp ( 2−θ 2 1−θ 2 ) θ 2 . Let the initial condition of system (8) below the separatrix curve L. From Figure 9, if the carrying capacity K of prey is small enough (i.e. K ≤ e hp ), predator will be extinct, while prey will survive (see Figures 8 a,b). If the carrying capacity K of prey is an appropriate value (i.e. e hp < K < e hp ( 2−θ 2 1−θ 2 ) θ 2 ), predator and prey can coexist (see Figure 8 c). If the carrying capacity K of prey is large enough (i.e. K > e hp ( 2−θ 2 1−θ 2 ) θ 2 ), predator and prey will oscillate periodically or be extinct (see Figures 8 d,e). System (8) undergoes a supercritical Hopf bifurcation at positive equilibrium E * 3 and exists a stable limit cycle around E * 3 . Compared with the global asymptotic stability of the positive equilibrium of system (6) with a > d, when considering system (6) with the group defense, by numerical simulation, we find that the limit cycle of system (8) disappears when a is large enough, and the origin is globally stable, that is both prey and predator will be extinct (Figure 8 e), which means that the group defense can destabilize the system. Hence, the group defense makes the dynamics behaviour of system (8) become more complicated.

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

Funding
This work was supported by the Natural Science Foundation of Fujian Province (2021J01613, 2021J011032).