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

1.   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, $\gamma(i)$ and $\eta\gamma_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 $\sqrt{A}$. Because $M$ is distributed on a flat two-dimensional area, $\sqrt{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\in[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

Where $x$ and $y$ are the density of prey and predator, respectively, $K$ is the level of fear. $(1-m)x (m\in[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, $\beta$ is the infection rate from susceptible prey to infected prey. $\alpha$ is the predator's capture rate for susceptible prey, $k_1$ is the half-saturation constant, $(1-m_1)X_S(m_1\in[0, 1))$ is the number of susceptible prey that can be caught by predator. $\delta$ is the predator's capture rate for infected prey, $(1-m)X_I (m\in[0, 1))$ is the number of infected prey that can be caught by predator. $c$ and $\xi$ 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. $\tau$ 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 disease-free 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.

2.   The system without delay

Firstly, we give system (1.4) without time delay as follow:

with the initial conditions

2.1.   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, \zeta)$, where $0 < \zeta\leq +\infty$. 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 $\Omega$ are uniformly bounded.

Proof. Let $X_S(t), X_I(t), Y(t)$ be the solution of system (2.1) under the initial condition.

Case I: If $X_S(0)\leq 1$, then $X_S(t)\leq 1$.

Assuming that it is not true, then $\exists$ $t_1$ and $t_2$ such that $X_S(t_1) = 1$ and $X_S(t) > 1$, $\forall$ $t\in (t_1, t_2)$. So

where $\phi(X_S(s), X_I(s), Y(s)) = \frac{r}{1+k_2Y(s)}-(b+\beta)X_I(s)-\frac{\alpha(1-m_1)Y(s)}{\sqrt{X_S(s)}+k_1(1-m_1)X_S(s)}-d_1$.

Therefore

Since $X_S(t_1) = 1$, $\phi(X_S(s), X_I(s), Y(s)) < 0$. Hence, $X_S(t) < 1$, which is a contradiction.

Case II: If $X_S(0) > 1$, then $\lim\limits_{t\to \infty}$sup$X_S(t)\leq 1$.

Suppose it is not true, then $X_S(t) > 1$, $\forall t > 0$ and $\phi(X_S(s), X_I(s), Y(s)) < 0$. Therefore, we have

Next, we construct a function $W(t)$ as follows:

where $P(b+\beta) = Q\beta, P = c R$. By differentiating (2.2) with respect to $t$, we get

Let $E = rX_S$, then $E$ attains its maximum value at $E_{max} = \frac{1}{P}$. Hence

that is

That is $\lim\limits_{t\to \infty}W(t)\leq \frac{1}{\kappa}$. So, all solutions of system (2.1) will enter into the region:

2.2.   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.

2.2.1.   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_A^2(t)$. Then system (2.1) transforms to the following form

The characteristic equation of system (2.3) at the trivial equilibrium $E_{0}$ is

Therefore, Eq (2.4) has three eigenvalues are $\lambda_1 = \frac{1}{2}(r-d_1), \lambda_2 = -d_2, \lambda_3 = -d_3$.

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.

2.2.2.   The boundary predator-free equilibrium

The boundary predator-free equilibrium $\widetilde{E}_1(\widetilde{X}_S, \widetilde{X}_I, 0)$ exists if $r > d_1$, here $\widetilde{X}_S = \frac{d_2}{\beta}$, $\widetilde{X}_I = \frac{r-d_1}{b+\beta}$. Now we prove the stability of the boundary predator-free equilibrium $\widetilde{E}_1$.

Theorem 2.2. The boundary predator-free equilibrium $\widetilde{E}_1$ of system (2.1) is unstable.

Proof. The Jacobian matrix of system (2.1) is as follows:

where

The characteristic equation of system (2.1) at the boundary predator-free equilibrium $\widetilde{E}_1$ is

Where $\; \; \widetilde{A}_{12} = \frac{-(b+\beta)d_2}{\beta}$, $\widetilde{A}_{21} = \frac{\beta(r-d_1)}{b+\beta}$, $\widetilde{A}_{33} = \frac{c\alpha(1-m_1)d_2}{\sqrt{d_2\beta}+k_1(1-m_1)d_2}-\frac{\xi\delta(1-m)(r-d_1)}{b+\beta}-d_3$.

Therefore, the first eigenvalue of Eq (2.6) is $\lambda_1 = \widetilde{A}_{33}$, and the other two eigenvalues are determined by the following equation

Thus, the boundary predator-free equilibrium $\widetilde{E}_1$ is locally asymptotically stable if $\widetilde{A}_{12}\widetilde{A}_{21} > 0$ $(r-d_1 < 0)$ holds along with $\frac{c\alpha(1-m_1)d_2}{\sqrt{d_2\beta}+k_1(1-m_1)d_2} < \frac{\xi\delta(1-m)(r-d_1)}{b+\beta}+d_3$, but $\widetilde{E}_1$ is unstable if $\widetilde{A}_{12}\widetilde{A}_{21} > 0$ does not hold. In fact, if $r-d_1 < 0$, then susceptible prey does not exist, thus $\widetilde{E}_1$ is unstable.

2.2.3.   The disease-free equilibrium

We can obtain the disease-free equilibrium $\overline{E}_2(\overline{X}_S, 0, \overline{Y})$, here $\overline{X}_S = \frac{d_3^2}{(c\alpha-k_1d_3)^2(1-m_1)^2}$. Further $\overline{Y}$ satisfies the following equation:

where $A_1 = -k_2\alpha(1-m_1), A_2 = -\{\alpha(1-m_1)+k_2d_1[\sqrt{\overline{X}_S}+k_1(1-m_1)\overline{X}_S]\}$, $A_3 = (r-d_1)\cdot[\sqrt{\overline{X}_S}+k_1(1-m_1)\overline{X}_S]$. Let $\Delta = A_2^2-4A_1A_3$.

Theorem 2.3. 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 $\overline{E}_{21}(\overline{X}_S, 0, \overline{Y}_1)$ and $\overline{E}_{22}(\overline{X}_S, 0, \overline{Y}_2)$, here $\overline{Y}_1 = \frac{\sqrt{\Delta}-A_2}{2A_1}, \overline{Y}_2 = \frac{-\sqrt{\Delta}-A_2}{2A_1}$.

(ii) if $\Delta > 0$, $A_3 < 0$, then system (2.1) has one disease-free equilibrium $\overline{E}_{21}(\overline{X}_S, 0, \overline{Y}_1)$, here $\overline{Y}_1 = \frac{\sqrt{\Delta}-A_2}{2A_1}$. \label{Theorem2.2.4}

In addition, the basic reproduction number is given by $R_0 = \frac{\beta\overline{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 $\overline{E}_{21}(\overline{X}_S, 0, \overline{Y}_1)$, and can use the similar methods to obtain the local stability of other disease-free equilibria.

Theorem 2.4. The disease-free equilibrium $\overline{E}_{21}$ is locally asymptotically stable if $(H_1)$ and $(H_2)$ hold along with $\beta\overline{X}_S-\delta(1-m)\overline{Y}-d_2 < 0$, but $\overline{E}_{21}$ is unstable if $(H_1)$ and $(H_2)$ does not hold.

Proof. The characteristic equation of system (2.1) at the disease-free equilibrium $\overline{E}_{21}$ is

where

Therefore, the first eigenvalue of Eq (2.8) is $\lambda_1 = J_{22} = \beta\overline{X}_S-\delta(1-m)\overline{Y}-d_2$, and the other two eigenvalues are determined by the following equation

All eigenvalues of Eq (2.9) have negative real parts if

$(H_1)$: $\frac{r}{1+k_2\overline{Y}} < d_1+\frac{(1-m_1)^2(c\alpha-d_3k_1)^3\overline{Y}}{2d_3c^2\alpha}$,

$(H_2)$: $\frac{c\alpha(1-m_1)\overline{Y}}{2\sqrt{\overline{X}_S}\Big[1+k_1(1-m_1)\sqrt{\overline{X}_S}\Big]^2}(-\frac{d_3}{c}-\frac{rk_2\overline{X}_S}{(1+k_2\overline{Y})^2}) < 0$

holds. Thus the disease-free equilibrium $\overline{E}_{21}$ is locally asymptotically stable if assumptions $(H_1)$ and $(H_2)$ hold along with $\beta\overline{X}_S-\delta(1-m)\overline{Y}-d_2 < 0$, but $\overline{E}_{21}$ is unstable if assumptions $(H_1)$ and $(H_2)$ does not hold.

2.2.4.   The positive equilibrium

Theorem 2.5. The positive equilibrium $E^*(X_S^*, X_I^*, Y^*)$ of system (2.1) exists if the assumption $(H_3)$ is true.

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

According to Eq (2.10), we can obtain that $X_I^* = \frac{1}{\xi\delta(1-m)}[\frac{c\alpha(1-m_1)X_S^*}{\sqrt{X_S^*}+k_1(1-m_1)X_S^*}-d_3], \; Y^* = \frac{\beta X_S^*-d_2}{\delta(1-m)}$, and $X_S^*(s = \sqrt{X_S^*})$ is the positive root of the equation:

where

Thus, the positive equilibrium $E^*(X_S^*, X_I^*, Y^*)$ of system (2.1) exists if the assumption

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.

3.   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 bifurcation of system (1.4). For convenience, let $\bar{X}_S(t) = X_S(t)-X_S^*, \bar{X}_I(t) = X_I(t)-X_I^*, \bar{Y}(t) = Y(t)-Y^*$, then we have the following linearized system

where

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.

Case I: $\tau = 0$

The Eq (3.2) is reduced to

where $p_{12} = p_2+s_2, \; \; p_{11} = p_1+s_1, \; \; p_{10} = p_0+s_0$.

So we know that all roots of Eq (3.3) have negative real parts if the following assumption

$(H_4): p_{12} > 0, \quad p_{10} > 0, \quad 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: $\tau\neq0$

Let $i\omega (\omega > 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

Let $\omega^2 = v$, then Eq (3.5) can be written as

Denote $f(v) = v^3+e_{12}v^2+e_{11}v+e_{10}$, then $f(0) = e_{10}, \lim\limits_{v\to +\infty}f(v) = +\infty$, $f^{'}(v) = 3v^2+2e_{12}v+e_{11}$.

After discussion about the roots of Eq (3.6) by the method in ^[47], there are the following conditions.

$(H_5):\; \; e_{10}\geq0, \quad \vartriangle = e_{12}^2-3e_{11}\leq 0$,

$(H_6):\; \; e_{10}\geq0, \quad \vartriangle = e_{12}^2-3e_{11} > 0, \quad v^* = \frac{-e_{12}+\sqrt{\vartriangle}}{3} > 0, \quad f_1(v^*)\leq 0,$

$(H_7):\; \; e_{10} < 0$.

Then we have the following lemma.

Lemma 3.1. For the polynomial Eq (3.6), we have the following results.

(1) If the assumption $(H_5)$ holds, then Eq (3.6) has no positive root.

(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 $\omega_{k} = \sqrt{v_{k}}$, $k = 1, 2, 3$. According to (3.4), if $v_{k} > 0$, the corresponding critical value of time delay is

where $A_{40} = s_{1}-p_{2}s_{2}, \; A_{20} = p_{2}s_{0}+p_{0}s_{2}-p_{1}s_{1}, \; A_{10} = -p_{0}s_{0}, B_{40} = s_{2}^2, \; B_{20} = s_{1}^2-2s_{0}s_{2}, \; B_{10} = s_{0}^2.$ Therefore, $\pm i\omega_{k}$ is a pair of purely imaginary roots of Eq (3.2) with $\tau = \tau_{k}^{(j)}$. And let $\tau_{0} = min_{k\in\{1, 2, 3\}}\{\tau_{k}^{(0)}\}, \omega = \omega_{k_0}$.

Lemma 3.2. Suppose that $(H_8):f^{'}(\omega^2)\neq0$, then the following transversality condition $\frac{\mathrm{d}(\mathrm{Re}\lambda)}{\mathrm{d}\tau}\Big|_{\lambda = i\omega}\neq0$ holds.

Proof. Differentiating Eq (3.2) with respect to $\tau$, and noticing that $\lambda$ is a function of $\tau$, then we obtain

Thus, we have

From Eq (3.4) and Eq (3.9),

It follows that $\frac{\mathrm{d}(\mathrm{Re}\lambda)}{\mathrm{d}\tau}\Big|_{\lambda = i\omega}\neq0$ and the proof is complete.

By Lemmas 3.1-3.2 and the Hopf bifurcation theorem ^[48,49], we have the following results.

Theorem 3.1. For system (1.4) with $\tau\neq0$, the following results are true.

(1) If the assumption $(H_5)$ holds, then the positive equilibrium $E^*(X_S^*, X_I^*, Y^*)$ is locally asymptotically stable for all $\tau\geq0$.

(2) If the assumption $(H_6)$ or the assumption $(H_7)$ and $(H_8)$ hold, then the positive equilibrium $E^*(X_S^*, X_I^*, Y^*)$ is locally asymptotically stable for all $\tau\in[0, \tau_{0})$ and is unstable for $\tau > \tau_{0}$. Furthermore, system (1.4) undergoes a Hopf bifurcation at the positive equilibrium $E^*(X_S^*, X_I^*, Y^*)$ when $\tau = \tau_{0}$.

Remark 3.1. The biological meaning of delay $\tau$ 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 $\tau$ is $\tau_0$. According to Theorem 3.1, system (1.4) can keep its stability if the minimum reaction time of susceptible prey does not exceed $\tau_0$, otherwise system (1.4) will change its stability.

4.   Global stability of equilibria

In this section, we will discuss the global stability of equilibria which are locally asymptotically stable under some conditions.

Theorem 4.1. The disease-free equilibrium $\overline{E}_2(\overline{X}_S, 0, \overline{Y})$ is globally asymptotically stable if the assumptions $(H_1)$, $(H_2)$ and $(H_9)$ are true, here

$(H_9):\; \; 1-m_1\geq0, 1-m\geq0, (c\alpha-k_1d_3)(1-m_1)\leq\sqrt{2}d_3$

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\geq T(T > 0)$ if $(H_9)$ holds. That is, $d_3\geq\frac{1}{\sqrt{2}(c\alpha-k_1d_3)(1-m_1)}$, then $\frac{\mathrm{d}V_1}{\mathrm{d}t}\leq0$. Let $M_1$ be the largest invariant set in $\left\{(X_S, X_I, Y)|\frac{\mathrm{d}V_1}{\mathrm{d}t} = 0\right\}$, we obtained that $\frac{\mathrm{d}V_1}{\mathrm{d}t} = 0$ if and only if $X_S = \overline{X}_S, X_I = \overline{X}_I = 0, Y = \overline{Y}$. Thus, $M_1 = \{\overline{E}_2\}$. By the LaSalle's invariance principle ^[50], the disease-free equilibrium $\overline{E}_2(\overline{X}_S, 0, \overline{Y})$ is globally asymptotically stable.

Theorem 4.2. The positive equilibrium $E^*(X_S^*, X_I^*, Y^*)$ is globally asymptotically stable if $(H_6)$ or $(H_7)$ and $(H_8)$, $(H_{10})$ hold, here

$(H_{10})$: $\; \; 1-m_1\geq0, 1-m\geq0, \frac{1}{2\sqrt{X_S^*}}-\sqrt{X_S^*}\leq0$ and $\frac{X_S(t-\tau)[1+k_2Y]}{X_S[1+k_2Y(t-\tau)]} = 1$.

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

Define

For Eq (4.5), for all $t\geq T(T > 0)$ and $1-m_1\geq0, 1-m\geq0, \frac{1}{2\sqrt{X_S^*}}-\sqrt{X_S^*}\leq0$, then $\frac{\mathrm{d}V_2^1}{\mathrm{d}t}\leq 0$. Then

Noting that $\frac{X_S}{1+k_2Y} = \frac{X_S^*}{1+k_2Y^*}$ at $E^*$, we have

When $n = 1$ in ^[51], if

then $\frac{\mathrm{d}V_3}{\mathrm{d}t}\leq 0$. Let $M_2$ be the largest invariant set in $\{(X_S, X_I, Y)|\frac{\mathrm{d}V_3}{\mathrm{d}t} = 0\}$. We obtained that $\frac{\mathrm{d}V_3}{\mathrm{d}t} = 0$ if and only if $X_S = X_S^*, X_I = X_I^*, Y = Y^*, \frac{X_S(t-\tau)[1+k_2Y]}{X_S[1+k_2Y(t-\tau)]} = 1$. Thus, $M_2 = \{E^*\}$. By the LaSalle's invariance principle ^[50], the positive equilibrium $E^*(X_S^*, X_I^*, Y^*)$ is globally asymptotically stable.

Remark 4.1. According to Eq (4.5), if $\frac{1}{2\sqrt{X_S^*}}-\sqrt{X_S^*}\leq0$, then $\frac{\mathrm{d}V_2^1}{\mathrm{d}t}\leq 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 $\frac{\mathrm{d}V_3}{\mathrm{d}t}\leq 0$, which is consistent with the numerical simulation in Section 7 (see Figure 6), thus the positive equilibrium $E^*$ is globally asymptotically stable.

5.   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.

6.   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 $\tau = \tau_{0}.$

Without loss of generality, let $\mu = \tau-\tau_{0}, \mu\in R, t = s\tau, X_S(s\tau) = \hat{X}_S(s), X_I(s\tau) = \hat{X}_I(s), Y(s\tau) = \hat{Y}(s).$ We denote $u_1(t) = X_S(t)-X_S^*, u_2(t) = X_I(t)-X_I^*, u_3(t) = Y(t)-Y^*,$ then system (1.4) can be written as a functional differential equation (FDE) in $C = C([-1, 0], R^3):$

where $u(t) = (u_1(t), u_2(t), u_3(t))^T\in C$, $u_t(\theta) = u(t+\theta) = (u_1(t+\theta), u_2(t+\theta), u_3(t+\theta))^T\in C$ and $L_\mu:C\rightarrow R^3$, $F: R\times C\rightarrow R^3$ are given by

and

where

Here

By the Riesz representation theorem ^[53], there exists a $3\times3$ matrix function $\eta(\theta, \mu)$ for $\theta\in[-1, 0)$ such that

In fact, we can choose

for $\phi\in C^1([-1, 0], R^3)$, we define

and

Then, Eq (6.1) can be rewritten as

For $\varphi\in C^1([-1, 0], (R^3)^*)$, where $(R^3)^*$ is the three-dimensional space of row vectors, we further define the adjoint operator $A^*$ of $A(0)$:

for $\phi\in C^1([-1, 0], R^3)$ and $\varphi\in C^1([-1, 0], (R^3)^*)$, we define the bilinear form

where $\eta(\theta) = \eta(\theta, 0), A = A(0)$ and $A^*$ are adjoint operators. By the discussion in Section 3, we know that $\pm i\omega\tau_{0}$ are eigenvalues of $A(0)$. Thus, they are also the eigenvalues of $A^*$.

We suppose that $q(\theta) = (1, q_2, q_3)^Te^{i\omega\tau_{0}\theta}$ is the eigenvector of $A(0)$ corresponding to the eigenvalue $i\omega\tau_{0}$, and $q^*(s) = D(1, q_2^*, q_3^*)e^{-i\omega\tau_{0}s}$ is the eigenvector of $A^*$ corresponding to the eigenvalue $-i\omega\tau_{0}$. By computation, we obtain

Then, from Eq (6.6), we get

Therefore, we choose $\bar{D} = \Big[1+\bar{q}_2^{*}q_2+\bar{q}_3^{*}q_3+\tau_{0}e^{-i\omega\tau_{0}}(b_{13}q_3+b_{11})\Big]^{-1}$, such that $\langle q^*(s), q(\theta)\rangle = 1, \langle q^*(s), \bar{q}(\theta)\rangle = 0$.

Next, let $u_t$ be the solution of Eq (6.5) when $\mu = 0$. We define

On the center manifold $C_0$, it comes to the conclusion that

where $z$ and $\bar{z}$ are local coordinates for $C_0$ in the direction of $q^*$ and $\bar{q}^*$. Note that $W$ is real if $u_t$ is real, we only consider real solutions. From Eq (6.8), we get

For a solution $u_t\in C_0$ of Eqs (6.3), (6.5), (6.6) and $\mu = 0$, we have

What's more, Eq (6.11) can be rewritten as follows,

where

It follows from Eqs (6.8) and (6.9) that

By Eqs (6.13) and (6.14), we obtain the following relevant parameters, which determine the direction and stability of Hopf bifurcation:

and

where $E_1 = \left(E_1^{(1)}, E_1^{(2)}, E_1^{(3)}\right)^T\in R^3$ and $E_2 = \left(E_2^{(1)}, E_2^{(2)}, E_2^{(3)}\right)^T\in R^3$ are also constant vectors and can be determined by the following equations, respectively,

where

Therefore, we can calculate $g_{21}$ and the following values:

which determine the properties of bifurcating periodic solutions at $\tau = \tau_{0}$. From the discussion above, we have the following results.

Theorem 6.1. For system (1.4), the direction of Hopf bifurcation is determined by the sign of $\mu_2$: if $\mu_2 > 0(\mu_2 < 0)$, then the Hopf bifurcation is supercritical (subcritical). The stability of the bifurcating periodic solutions is determined by the sign of $\beta_2$: if $\beta_2 < 0 (\beta_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).

7.   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).

7.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, \delta, \beta, 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, \beta$ 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 $\delta$ has no correlation with the output of system.

7.2.   The stability and Hopf bifurcation

In this subsection, we will analyze the stability of system (2.1) and assume that the initial values of system (2.1) are $X_S(0) = 0.6, X_I(0) = 0.1, Y(0) = 0.2$. First, we choose a set of parameter values: $r = 0.17, k_2 = 10, b = 0.5, \beta = 0.6, \alpha = 0.6, k_1 = 1, d_1 = 0.01, d_2 = 0.9, \delta = 1.5, c = 0.89, \xi = 0.14, d_3 = 0.2, m_1 = 0.2, m = 0.42$. According to the Theorem 2.5, there is a positive equilibrium $E^*(0.6002, 0.0355, 0.0971)$. We find that the positive equilibrium $E^*$ is locally asymptotically stable under this set of parameter values, the susceptible prey, infected prey and predator population have stable dynamics behavior(see Figure 3). Furthermore, according to the discussion in Section 2, we can get $R_0 = \frac{\beta\overline{X}_S}{d_2} = 0.37 < 1$(because $\overline{X}_S = 0.56$ from Section 2.2.3), which indicates that the disease dies out from the system.

Next, for the given parameters, we get $\omega = 0.1486$ and $\tau_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 $\tau = 0.5 < \tau_0 = 1.059$, the susceptible prey, infected prey and predator population have stable dynamics behavior (see Figure 4). Then, as the value of $\tau$ gradually increases, we choose the value of $\tau$ as $\tau = 2 > \tau_0 = 1.059$ and obtain that $C_1(0) = -4.5406-1.7383i, \mu_2 = 5.1018 > 0, \beta_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 $\tau$ exceeds its threshold value $\tau_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 $\tau_0 = 1.059$ and a stable limit cycle appeared in phase portrait.

Then, we choose $r = 0.6, d_2 = 0.5, \delta = 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\times10^{-3}, s = 10, w = 1\times10^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 $\tau = 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 $\tau$ can also determine the global asymptotic stability of system (1.4).

7.3.   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.

7.4.   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.

7.5.   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 $\delta$ 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 $\delta$ 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 $\delta$ increases from Figure 14 (b) and Figure 14 (c). When it reached its critical value $\delta_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.

8.   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 $\widetilde{E}_1$, a disease-free equilibrium $\overline{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 $\tau$ as bifurcation parameters, the stability of this system is discussed. If time delay $\tau$ is greater than its corresponding critical values, the Hopf bifurcation occurs at the positive equilibrium $E^*$. Next, the global stability of the equilibrium $\overline{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 $\delta$ 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 $\delta$ of predators to infected prey has a positive effect on the increase in the number of susceptible prey. The catch rate $\delta$ 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 $\tau = 0.5 < \tau_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 $\tau = 2 > \tau_0 = 1.059$. Thus, Hopf bifurcation takes place when $\tau = \tau_0$ (see Figure 5). In other words, the minimum reaction time of susceptible prey cannot exceed $\tau_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

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant Nos. 11661050 and 11861044), and the HongLiu First-class Disciplines Development Program of Lanzhou University of Technology. We also would like to thank two anonymous reviewers for their valuable comments and suggestions.

Conflict of interest

All authors declare no conflicts of interest in this paper.