Dynamics of a delayed predator-prey system with fear effect, herd behavior and disease in the susceptible prey

In this article, a delayed predator-prey system with fear effect, disease and herd behavior in prey incorporating refuge is established. Firstly, the positiveness and boundedness of the solutions is proved, and the basic reproduction number R0 is calculated. Secondly, by analyzing the characteristic equations of the system, the local asymptotic stability of the equilibria is discussed. Then taking time delay as the bifurcation parameters, the existence of Hopf bifurcation of the system at the positive equilibrium is given. Thirdly, the global asymptotic stability of the equilibria is discussed by constructing a suitable Lyapunov function. Next, the direction of Hopf bifurcation and the stability of the periodic solution are analyzed based on the center manifold theorem and normal form theory. What’s more, the impact of the prey refuge, fear effect and capture rate on system is given. Finally, some numerical simulations are performed to verify the correctness of the theoretical results.


Introduction
The relationship between predator and prey is the main field in natural ecology, and the mathematical models describing the interaction of two populations and the corresponding behaviors of those models are also topics which lots of scholars focus on, for example, the classic predator-prey model was first proposed by Lotka [1] and Volterra [2]. Due to different ecological and environmental factors, many species are affected by diseases [3]. Thus, eco-epidemiology is also one topic of concerns for many researchers, which includes ecology and epidemiology with the main purpose of studying how diseases spread. Chattopadhyay and Arino [4] first proposed a predator-prey model with disease in prey population as follows: where s and i are the number of sound prey (also named the susceptible prey) and infected prey, respectively, y is the number of predator population, γ(i) and ηγ 1 (s) are the predator response functions. They mainly studied the persistence and extinction of the population, the existence and stability of the equilibrium and Hopf bifurcation occurring at the positive equilibrium of system (1.1). Many scholars studied the predator-prey models with disease and time delay based on system (1.1) [5][6][7][8][9][10][11][12][13].
In the natural world, in order to avoid being caught by predators, prey often likes to form a herd for life. This herd can be used as a protective shield against predators, because predators cannot easily attack the prey in the herd. For the predator population, it is obviously easy to catch the prey who does not form a herd. Let M(t) be the density of the population that forms the herd at time t and occupies area A, the population density of the individuals staying outside the herd is proportional to the length of the herd at time t, that is, its length is proportional to √ A. Because M is distributed on a flat two-dimensional area, √ M is an individual at the edge of the patch. At present, many scholars have studied the population with the herd behavior [14][15][16][17][18]. Recently, Saha and Samanta [19] proposed a predator-prey system with herd behavior and disease and incorporating prey refuge, where X S (T ) and X I (T ) are the density of susceptible prey and infected prey at time T , respectively, Y(T ) is the density of predator at time T . K 0 is the Allee threshold of susceptible prey population in absence of predator, B is the susceptible prey feels the population pressure of both susceptible and infected prey, (1 − m)X I (m ∈ [0, 1)) is the number of infected prey that can be caught by predator. The authors obtained some dynamical behaviors, such as the boundedness and positivity of the solution, the local stability of the equilibrium. In addition, in the biological system, the fear effect is a universal phenomenon. Almost every creature can feel the crisis of being captured for prey and make a difference. For example, when the prey feel the risk of predation, they may abandon their previous high-risk predation habits, and re-choose a relatively low-risk predation habit, which deplete the energy of the prey [20]. And more and more scholars have done a lot of research on the fear effect between predator and prey [21][22][23][24][25][26]. Furthermore, due to the inherent fear of prey on predators, we need to take certain protection measures for the prey species to prevent all prey from being caught by the predators. Recently, some scholars have studied the predator-prey model with prey refuge, and showed that the use of refuge has a positive effect on the prey [27][28][29][30][31][32]. For example, Zhang et al. [30] studied a predator-prey model with fear effect and incorporating prey refuge (1.3) Where x and y are the density of prey and predator, respectively, K is the level of fear. (1 − m)x(m ∈ [0, 1)) is the number of prey that can be caught by predator. They found that the fear effect can reduce the population density of predators and change stability of system by eliminating the existence of periodic solutions. However, some behaviors of species do not appear immediately in natural ecology, that is, there is a time delay [33][34][35][36]. Delay differential equations can reflect reality more truly and its dynamic behavior is more complicated. It is very suitable for scholars to combine time delay with stage structure, disease and other factors to study the influence of time delay [37][38][39][40]. Furthermore, the functional response of predators to prey density is essential in predator-prey system, and it can enrich the dynamics of predator-prey systems. In ecology, many factors could affect functional responses, such as prey escape ability, predator hunting ability and the structure of the prey habitat [41]. Generally, functional responses can be divided into the following types: prey-dependent (such as Holling I-III [42]) and predator-dependent (such as Beddington-DeAngelis [43], Crowley-Martin [44]). Recently, more and more scholars have studied predator-prey system with functional response functions and time delays [31,39,45].
We assume that the number of susceptible prey is proportional to the number of existing susceptible prey, and the number of infected prey is proportional to the number of existing infected prey. Similarly, the number of predator is directly proportional to the number of existing predator. Motivated by the literatures [4,19,30], we will take fear effect, herd behavior, disease, Holling type II functional response, time delay and refuge for prey into model (1.2) as follows with the initial conditions where X S (t) and X I (t) are the density of susceptible prey and infected prey at time t, respectively, Y(t) is the density of predator at time t. r is the intrinsic growth rate of susceptible prey, k 2 is the level of fear of predators by the susceptible prey, b is the susceptible prey feels the population pressure of both susceptible and infected prey, β is the infection rate from susceptible prey to infected prey. α is the predator's capture rate for susceptible prey, k 1 is the half-saturation constant, (1 − m 1 )X S (m 1 ∈ [0, 1)) is the number of susceptible prey that can be caught by predator. δ is the predator's capture rate for infected prey, (1 − m)X I (m ∈ [0, 1)) is the number of infected prey that can be caught by predator. c and ξ are the conversion rates of the predator, d 1 , d 2 and d 3 are the natural death rates of susceptible prey, infected prey and predator, respectively. τ is the reaction time of the susceptible prey when they feel the predation crisis. All parameters are positive constants. The organization of this article is as follows: In Section 2, the positiveness and boundedness of system (1.4) without time delay, the basic reproduction number R 0 and the existence and stability of the equilibria are given. In Section 3, the stability of the positive equilibrium and the existence of Hopf bifurcation are studied, respectively. In Section 4, we give the global stability of the diseasefree equilibrium and the positive equilibrium. The Lyapunov exponent of system (2.1) is calculated in Section 5. In Section 6, the direction and the stability of Hopf bifurcation are studied based on the center manifold theorem and the normal form theory. To support our theoretical results, some numerical simulations are given in last section.

Positivity and boundedness of solutions
In natural ecology, positiveness means that the population can survive, and boundedness means that the resources of the population are limited. Therefore, we will discuss the positiveness and boundedness of system (2.1). Lemma 2.1. All solutions of system (2.1) that start in R 3 + remain positive for all time. Proof. Because system (2.1) is continuous and satisfy the local Lipschitz condition on the continuous function space C, the solution (X S (t), X I (t), Y(t)) of system (2.1) exists and is unique with the positive initial conditions (X S (0), X I (0), Y(0)) on [0, ζ), where 0 < ζ ≤ +∞. From the first equation of system (2.1), we can get that is: Similarly, we have that Lemma 2.2. All solutions of system (2.1) which in Ω are uniformly bounded.
Proof. Let X S (t), X I (t), Y(t) be the solution of system (2.1) under the initial condition.
Next, we construct a function W(t) as follows: where P(b + β) = Qβ, P = cR. By differentiating (2.2) with respect to t, we get Let E = rX S , then E attains its maximum value at E max = 1 P . Hence dW dt That is lim t→∞ W(t) ≤ 1 κ . So, all solutions of system (2.1) will enter into the region:

Existence and local stability of equilibria
In this subsection, we will discuss the existence and local stability of equilibria of system (2.1). By calculations, system (2.1) has four equilibria as follows.

The trivial equilibrium
Theorem 2.1. The trivial equilibrium E 0 of system (2.1) is unstable.
Proof. In order to analyze the local stability of the trivial equilibrium E 0 (0, 0, 0), we redefine the variables X A (t) as X S (t) = X 2 A (t). Then system (2.1) transforms to the following form (2. 3) The characteristic equation of system (2.3) at the trivial equilibrium E 0 is Therefore, Eq (2.4) has three eigenvalues are Thus, when r − d 1 < 0, the trivial equilibrium E 0 is locally asymptotically stable. But the trivial equilibrium E 0 is unstable when r − d 1 > 0. In fact, if r − d 1 < 0, then susceptible prey population not exist, thus E 0 always unstable.

The boundary predator-free equilibrium
The boundary predator-free equilibrium E 1 ( X S , X I , 0) exists if r > d 1 , here X S = d 2 β , X I = r−d 1 b+β . Now we prove the stability of the boundary predator-free equilibrium E 1 .
Theorem 2.2. The boundary predator-free equilibrium E 1 of system (2.1) is unstable.
Proof. The Jacobian matrix of system (2.1) is as follows: The characteristic equation of system (2.1) at the boundary predator-free equilibrium E 1 is Therefore, the first eigenvalue of Eq (2.6) is λ 1 = A 33 , and the other two eigenvalues are determined by the following equation Thus, the boundary predator-free equilibrium E 1 is locally asymptotically stable if A 12 A 21 > 0 (r − d 1 < 0) holds along with cα(1−m 1 )d 2 √ We can obtain the disease-free equilibrium Further Y satisfies the following equation: where The disease-free equilibria of system (2.1) are as follows.
(i) if A 3 > 0, then system (2.1) has two disease-free equilibria E 21 (X S , 0, Y 1 ) and In addition, the basic reproduction number is given by R 0 = βX S d 2 . The R 0 mainly represents the number of new infective generated from a single infected individual [46]. If R 0 < 1, then the disease dies out from the system, but it remains the endemic in the host population if R 0 > 1. Now we prove the stability of the disease-free equilibrium E 21 (X S , 0, Y 1 ), and can use the similar methods to obtain the local stability of other disease-free equilibria.
Proof. The characteristic equation of system (2.1) at the disease-free equilibrium E 21 is Therefore, the first eigenvalue of Eq (2.8) is , and the other two eigenvalues are determined by the following equation All eigenvalues of Eq (2.9) have negative real parts if holds. Thus the disease-free equilibrium E 21 is locally asymptotically stable if assumptions (H 1 ) and Proof. We assume that E * (X * S , X * I , Y * ) is a positive equilibrium of system (2.1), then X * S , X * I and Y * satisfy the following equations (2.10) According to Eq (2.10), we can obtain that X * is the positive root of the equation: Thus, the positive equilibrium is true.
Next, we discuss the dynamical behavior of the positive equilibrium E * (X * S , X * I , Y * ) of system (2.1) with time delay.

The existence of Hopf bifurcation of system with delay
In this section, we discuss the local stability of the system at the positive equilibrium and the existence of Hopf Therefore, the characteristic equation corresponding to system (3.1) can be given where In order to study the distribution of the root of Eq (3.2), we will discuss it in the following cases.
So we know that all roots of Eq (3.3) have negative real parts if the following assumption (H 4 ) : p 12 > 0, p 10 > 0, p 12 p 11 > p 10 , holds. That is, the positive equilibrium E * (X * S , X * I , Y * ) of system (3.1) is locally asymptotically stable if the condition (H 4 ) is satisfied. Case II: τ 0 Let iω(ω > 0) be a root of Eq (3.2). By separating real and imaginary parts, it follows that Adding squares of Eq (3.4), we can get where 11 . After discussion about the roots of Eq (3.6) by the method in [47], there are the following conditions.
Then we have the following lemma.
Lemma 3.1. For the polynomial Eq (3.6), we have the following results.
(2) If the assumption (H 6 ) or the assumption (H 7 ) holds, then Eq (3.6) has at least one positive root.
Without loss of generality, we assume that Eq (3.6) has three positive roots defined as v 1 , v 2 and v 3 . Then Eq (3.5) has three positive roots ω k = √ v k ,k = 1, 2, 3. According to (3.4), if v k > 0, the corresponding critical value of time delay is where Proof. Differentiating Eq (3.2) with respect to τ, and noticing that λ is a function of τ, then we obtain From Eq (3.4) and Eq (3.9), It follows that d(Reλ) dτ λ=iω 0 and the proof is complete.
Remark 3.1. The biological meaning of delay τ is the reaction time of the susceptible prey when they feel the predation crisis in our system (1.4). More detailed discussions on similar time delay can be found in [45]. According to our analysis of system (1.4), we find the minimum critical value of τ is τ 0 . According to Theorem 3.1, system (1.4) can keep its stability if the minimum reaction time of susceptible prey does not exceed τ 0 , otherwise system (1.4) will change its stability.

Global stability of equilibria
In this section, we will discuss the global stability of equilibria which are locally asymptotically stable under some conditions.
Proof. We construct a suitable Lyapunov function Here, V 1 (X S , X I , Y) is a positive definite function for all (X S , X I , Y). Then differentiating (4.1) with respect to t, we get Therefore, for all t ≥ T (T > 0) if (H 9 ) holds. That is, By the LaSalle's invariance principle [50], the disease-free equilibrium E 2 (X S , 0, Y) is globally asymptotically stable.
Proof. We construct a suitable Lyapunov function: Then differentiating (4.2) with respect to t, we get

Now let
Then differentiating (4.3) with respect to t, we get For Eq (4.5), for all t ≥ T (T > 0) and 1 − m Noting that X S 1+k 2 Y = X * S 1+k 2 Y * at E * , we have When n = 1 in [51], if then dV 3 dt ≤ 0. Let M 2 be the largest invariant set in {(X S , X I , Y)| dV 3 dt = 0}. We obtained that dV 3 dt = 0 if and only if X S = X * S , By the LaSalle's invariance principle [50], the positive equilibrium E * (X * S , X * I , Y * ) is globally asymptotically stable. dt ≤ 0, here X * S is the positive root of Eq (2.11). That is, X * has a certain mathematical expression, which is greater than or equal to 0.5. Furthermore, if Eq (4.7) is also true, then dV 3 dt ≤ 0, which is consistent with the numerical simulation in Section 7 (see Figure 6), thus the positive equilibrium E * is globally asymptotically stable.

Numerical calculation of Lyapunov exponents
In a dynamic system, the Lyapunov exponent is an important indicator to measure the dynamics of the system [52]. It indicates the average exponential rate of convergence or divergence of the system between adjacent orbits in phase space.
Theorem 5.1. If the Lyapunov exponents are all negative, then system (2.1) is globally asymptotically stable, but if one is positive, then system (2.1) is unstable.
Proof. The Jacobian matrix corresponding of system (2.1) is Eq (2.5), on the basis of given some parameters values, initialize the Jacobian matrix (2.5), and initialize the three Lyapunov exponents to be calculated, we define a starting time: where s is the number of steps in each evolution, t s is the time step, I 1 is the number of evolutions and w is the total number of cycles. In Eq (2.5), we let: Next, use the Schmidt orthogonalization method to orthogonalize Eq (5.2), we can obtain that By calculation, we make the orthogonalized matrix A = [a 1 , a 2 , a 3 ], and assume that the matrix after orthogonalization is: Then take the modulus length of the three vectors as: Finally, the three Lyapunov exponents of system (1.4) are given: Therefore, if the Lyapunov exponents are all negative, then system (2.1) is globally asymptotically stable, but if one is positive, then system (2.1) is unstable.

Direction and stability of Hopf bifurcation
In this part, we will study the direction of Hopf bifurcation and the stability of bifurcating periodic solutions of system (1.4). The theoretical approach is the normal form theory and center manifold theorem [48]. Throughout this section, we assume that system (1.4) undergoes a Hopf bifurcation at τ = τ 0 .
Theorem 6.1. For system (1.4), the direction of Hopf bifurcation is determined by the sign of µ 2 : if µ 2 > 0(µ 2 < 0), then the Hopf bifurcation is supercritical (subcritical). The stability of the bifurcating periodic solutions is determined by the sign of β 2 : if β 2 < 0(β 2 > 0), then the bifurcating periodic solutions are stable (unstable). The period of the bifurcating periodic solutions is determined by the sign of T 2 : if T 2 > 0(T 2 < 0), then the bifurcating periodic solutions increase (decrease).

Numerical simulations
In this section, we present some numerical simulations of system (1.4) and system (2.1) to support our theoretical results. First, we will conduct the sensitivity analysis of some parameters of system (2.1).

Sensitivity analysis
To determine the parameters that have significant impact on output variables of system (2.1), we conduct global sensitivity analysis of some parameters, many scholars have studied the sensitivity analysis of the parameters in system [54,55]. We calculate partial rank correlation coefficients (PRCCs) between the parameters r, δ, β, m 1 , m and k 2 from system (2.1). Nonlinear and monotone relationships are observed with the input parameters of system (2.1), which is a prerequisite for computing PRCCs. Then, a total of 1000 simulations of the model per Latin hypercube sampling (LHS) were carried out using the baseline values tabulated in Table 1. According to the parameter values in Table 1, we analyze the influence of some parameters in the system on the correlation of infected prey. By sampling these parameters for 1000 times and a scatter plot with a fixed time point of 80, we obtain the sampling results in Figure 1 and scatter plot in Figure  2. Monotonical increasing (decreasing) indicates a positive (negative) correlation of the parameter with the model output. It is known from Figure 2 that several selected parameters exhibit periodic correlation, we can know that the parameters r, β and m 1 show a positive correlation with the output of system, the parameters m, k 2 show a little negative correlation with the output of system, and the parameter δ has no correlation with the output of system.
Next, for the given parameters, we get ω = 0.1486 and τ 0 = 1.059 when system (1.4) with time delay.
According to the Theorem 3.1, we can obtain that the positive equilibrium E * (0.6002, 0.0355, 0.0971) of the system (1.4) is locally asymptotically stable when τ = 0.5 < τ 0 = 1.059, the susceptible prey, infected prey and predator population have stable dynamics behavior (see Figure 4). Then, as the value of τ gradually increases, we choose the value of τ as τ = 2 > τ 0 = 1.059 and obtain that C 1 (0) = −4.5406 − 1.7383i, µ 2 = 5.1018 > 0, β 2 = −9.0812 < 0, T 2 = −2.0684 < 0 by Eq (6.18). Moreover, we know that system (1.4) undergoes a supercritical Hopf bifurcation when the value of τ exceeds its threshold value τ 0 . In addition, the bifurcating periodic solutions of system (1.4) are stable, and the period of bifurcating periodic solutions decreases. We clearly see from Figure 5 that the positive equilibrium E * (0.6002, 0.0355, 0.0971) is destabilized through a Hopf bifurcation in τ 0 = 1.059 and a stable limit cycle appeared in phase portrait.  Then, we choose r = 0.6, d 2 = 0.5, δ = 3, the initial conditions of the variables and other parameter values are the same as those given above. According to the Theorem 4.2, system (1.4) has a positive equilibrium E * (0.9023, 0.1255, 0.1943). By numerical simulations, we find that the positive equilibrium E * is global asymptotically stable under this set of parameter values (see Figure 6). Figure 6 (b) shows that the trajectory of the system (1.4) tends to the positive equilibrium E * for different initial values at any time. Finally, we give some parameter values: t * = 0, t s = 1 × 10 −3 , s = 10, w = 1 × 10 5 , then the Lyapunov exponents have been derived numerically of system (2.1) for different species (see Figure 7 (a)).
According to Theorem 5.1, all Lyapunov exponents are negative (L 1 = −0.0871, L 2 = −1.0740, L 3 = −0.8451), thus system (2.1) is globally asymptotically stable. We draw the maximum Lyapunov exponent of system (1.4) for τ = 2 (see Figure 7 (b)). In the figure, positive values of the maximum Lyapunov exponent indicates that system (1.4) is unstable. Figure 7 (b) corresponding to Figure 6, the value of time delay τ can also determine the global asymptotic stability of system (1.4).

The impact of the susceptible prey refuge
In this subsection, we will discuss how the susceptible refuge parameter m 1 affects on each of the population when the positive equilibrium exists and is locally asymptotically stable. We study the influence of different parameter values of m 1 on the population of prey and predator, which can be seen in Figure 8 and Figure 9. The dynamic response diagram of system (2.1) under different parameter values m 1 is shown in Figure 8, which corresponds to Figure 9. Figure 9 shows that the solution curve of state variables of system (2.1) under different values of m 1 . It can be seen from Figure 9 that as the value of the m 1 increases, the number of susceptible prey and infected prey gradually increases, which results in an immediate reduction in the number of predators, but system (2.1) is locally asymptotically stable all time. This shows that taking refuge measures has a positive effect on the prey population, but has a negative effect on the growth of the predator population. This conclusion is consistent with the laws of nature.

The impact of the infected prey refuge
In this subsection, we will discuss how the infected refuge parameter m affects on each of the population when the positive equilibrium exists and is locally asymptotically stable. We study the influence of different parameter values of m on the population of prey and predator, whose results are shown in Figure 10 and Figure 12.   The dynamic response diagram of system (2.1) under different parameter values m is shown in Figure 10, which corresponds to Figure 12. Figure 12 shows that the solution curve of state variables of system (2.1) under different values of m. It can be obtained from Figure 12 that as the value of the m increases from 0 to 0.8, the number of susceptible prey and predator decreases, but the number of infected prey does not change significantly. This is mainly impact due to disease factors.

The impact of the fear effect
In this subsection, we will discuss how the fear effect parameter k 2 affects on each of the population when the positive equilibrium exists and is locally asymptotically stable. We study the influence of different parameter values of k 2 on the population of prey and predator (seen Figure 11 and Figure 13).  Figure 11 shows that the dynamic response diagram of system (2.1) corresponding to the fear effect parameter k 2 in the time interval [0, 100]. It can be seen that the influence of the fear effect k 2 on the population is negative for prey and predator. That is, when the value of the fear parameter k 2 increases from 0 to 40, the number of prey and predator population decreases (see Figure 13), but the level of fear effect k 2 does not affect the stability of the system (2.1). This conclusion is consistent with the conclusions obtained from system (2.1) state response diagrams in Figure 13.  Finally, we give the effect of the predator's catch rate δ for diseased prey on the number of species in system (see Figure 14). It can be seen from Figure 14 (a) that the number of susceptible prey gradually increase as δ increases, which indicates that the capture rate has a positive effect on the number of susceptible prey. However, we can obtain that the population of infected prey and predator first increases and then decreases when δ increases from Figure 14 (b) and Figure 14 (c). When it reached its critical value δ C = 28, the number of infected prey and predator begins to decrease. At this time, the number of predators tends to become extinct, because disease factors inhibit the growth of the number of predators and infected prey.

Discussion and conclusions
In this paper, a predator-prey model with fear effect and disease in prey incorporating prey refuge is studied. To make such system be more realistic, the effect of delay denoting the reaction time of the susceptible prey who feels the predation crisis on the system is considered, and the herd behavior of prey is also considered. We assume that the predator's capture rate on prey population is in line with Holling-II type functional response function. According to calculation, the basic reproduction number is given that R 0 = 0.37 < 1, which indicates the disease dies out from system. And system (2.1) has a trivial equilibrium E 0 , a boundary predator-free equilibrium E 1 , a disease-free equilibrium E 2 and a unique positive equilibrium E * . By discussing the distribution of the roots of the characteristic equation of system (2.1) and taking the time delay τ as bifurcation parameters, the stability of this system is discussed. If time delay τ is greater than its corresponding critical values, the Hopf bifurcation occurs at the positive equilibrium E * . Next, the global stability of the equilibrium E 2 and E * is proved. Furthermore, we discuss the direction and stability of Hopf bifurcation by using the center manifold theorem and the normal form theory. Finally, the stability and bifurcations affecting the prey and predator population are supported by numerical simulations.
In the absence of time delay, we select m 1 , m, k 2 and δ as the analysis parameters, and discuss the influence of these parameters on stability of system (2.1). We found that the parameter m 1 has a positive effect on the increase of the prey population, and vice versa for the predator population (see Figure 8). The parameter m has a negative effect on the change of the prey population, but does not change significantly for the predator population (see Figure 10). In the literature [19], they studied that if the prey is infected with the disease, then the three species do not coexist. Our paper mainly introduces the influence of the fear effect k 2 on system (2.1). For the parameter k 2 of the level of fear of prey on the predator, as the level of fear increases, the number of prey and predator population decreases (see Figure 11). But the level of fear effect k 2 does not affect the stability of system (2.1). The catch rate δ of predators to infected prey has a positive effect on the increase in the number of susceptible prey. The catch rate δ inhibits the increase in the number of predators and infected prey due to disease factors (see Figure 13). In addition, the stability of system (1.4) with time delay is discussed. The study shows that when τ = 0.5 < τ 0 = 1.059, system (1.4) is locally asymptotically stable at the positive equilibrium E * (see Figure 4), but system (1.4) is unstable when τ = 2 > τ 0 = 1.059. Thus, Hopf bifurcation takes place when τ = τ 0 (see Figure 5). In other words, the minimum reaction time of susceptible prey cannot exceed τ 0 , otherwise the stability of the system (1.4) will change. At the same time, the direction of Hopf bifurcation and the stability of periodic solutions are discussed by the center manifold theorem and normal form theory.
In future, in order to obtain the greatest economic benefits, we will consider the effective harvest of the susceptible prey and predator. In addition, by considering the intrapsychic competition in the susceptible prey, we have the following model