Bifurcations of a delayed predator-prey system with fear, refuge for prey and additional food for predator

1.   Introduction

Predator-prey systems are almost ubiquitous in the real world. Usually, there are two kinds of interspecific relationships between prey and predator ^[1,2]. One is the familiar and lethal predation relationship, i.e., the prey is killed directly by the predator for its survival ^[3,4,5], and the other is the non-consumptive relationship of the prey's fear of the predator ^[6,7,8].

In face of the threat of predators, the prey will be circumspect or vigilant to the predator's sounds or smells and reduce foraging activity. After a long time, the physiological characteristics of the prey change, and the reproduction rate declines, which has been demonstrated by several examples, such as elks and wolves, hares and dogs, mule deer and mountain lions, etc. ^[9,10,11]. Recently, Wang et al. ^[12] proposed a predator-prey model with the prey's fear of the predator. They observed that the growth rate of the prey was decreased by fear, while it could stabilize the system. Through subcritical Hopf bifurcation, an inferior level of fear could cause a multi-stability scenario. Currently, many ecological models with fear effect have been investigated, and nice results have been obtained; see ^[13,14,15] and references cited therein.

Due to fear of the predator, the prey often searches for a refuge area to avoid the predator's killing ^[16,17,18]. Actually, almost every prey population has refuges to maintain higher biomass. For example, rats have tall grass as their refuge area to hide from predators, such as owls and cats and to keep a higher population density. In a lake ecosystem, the phytoplankton often use the benthic sediments as their refuge area to lay eggs and avoid the zooplankton's predation. What's more, they can use water stratification as temporary refuge to escape from predation by zooplankton ^[19]. By empirical study, Mullin et al. ^[20] found that with the help of a refuge area, prey was protected by decreasing the oscillatory tendency of species interactions and stabilizing the community's equilibrium state. In addition to the protection of prey, it can affect the predator's survival and destabilize the long-term behaviors of the ecosystem; the refuge has two sides for the survival of species ^[21].

Taking into account the anti-predation behavior of prey, sometimes the predator will lack sufficient food to maintain normal survival, so additional food supplement is very urgent to avoid the oscillation and extinction of the species ^{[22,23,24,25,26]}. Additional food can distract the predator's attention from prey and decrease the predation rate. Simultaneously, because of sufficient additional food, the predator's reproduction rate may increase, resulting in the reduction of the prey's density or even its extinction. For example, in a biological control strategy, usually the additional food is provided for the predator (natural enemies), so that the predation rate of the predator on the prey (pests) is increased. Thus the prey density is reduced, the attack of prey on seeds declines largely, and the loss of economy scales down. In addition, a field experiment was executed by the Moorland Working Group ^[27] on Lang-holm Moor in Scotland in the spring and summer of 1998 and 1999. They released sufficient food for the hen harriers (predator), and then the prey (chicks) was predated less than before. The predation rate declined from 3.7 chicks per 100 hours to 0.5 chicks per 100 hours. Additional food supplement has both positive and negative functions in system balance. It can not only help to eliminate the prey, but also eliminate the predator itself because of too much reproduction and interspecific competition. By manipulating the quality and quantity of additional food to the predator, a desired level of species within an acceptable range can be driven, and the dynamics of an ecosystem can be controlled in reasonable ways ^[28,29]. Now, the research of additional food supplement has been one of the hottest topics.

On the other hand, time delay is almost everywhere in any biological proceeding. For example, there always exist time delays for predator's gestation or digestion and so on ^{[30,31,32,33]}. The time delay in fear is quite common. Generally, when the prey sense the risk of predation, it will take some time to assess the risk level and perform some anti-predation behaviors. There is a time delay in the reduction of the growth rate of the prey ^[34,35,36].

The objective of our research is to construct a predator-prey model with delayed fear, refuge for prey and additional food for predator, and we investigate how these factors affect system stability.

We structure our research as follows. First, the model is constructed in Section 2, and then some preparations are given in Section 3. Next, we discuss the asymptotical stability of equilibria and Hopf bifurcation of parameters in Section 4. In Section 5, we carry out some numerical simulations and explore the effects of the above factors on the system stability. Finally, in Section 6, we give a concluding remark to end this paper.

2.   Model construction

We start with the following steps to construct our model.

(i) The prey species obeys the usual logistic growth driven by the equation $\frac{dx}{dt} = x(b-cx).$ When the predator appears, the functional response of the predator to the prey follows Holling type II. That is, the basic prey-predator model is

(ii) Due to the prey's fear of the predator and the prey's anti-predation behaviors, we have

where $\varsigma+\frac{\sigma(1-\varsigma)}{\sigma+y }$ is a fear level which affects the growth rate of prey (more details in ^[8]). $\frac{1}{1+a}$ is caused by the anti-predation activities ^[32]. It indicates the reduction of predation rate. Considering the anti-predation effect of prey, the mortality rate of predator is assumed to be nonlinear, denoted by $e+\frac{\omega-e}{1+y}$. For detailed biological explanations see ^[8].

(iii) As an anti-predation behavior, the prey is assumed to search for a refuge area to protect itself. The biomass $\kappa x$ lives in refuge area, and $(1-\kappa)x$ is exposed to predator. Meanwhile, additional food is supplied to the predator to maintain the survival of the predator. Then, we get the following model:

where $\frac{q(1-\kappa)y}{d+(1-\kappa) x+\gamma\delta A}$ represents the functional response when additional food appears. The biological explanation is referred in ^[23].

(iv) Incorporating the delay of fear on the prey's growth rate into (2.1), we have

The novelty of (2.2) is explained as follows. First, due to the powerful predator's threat on the prey, the growth rate of prey is reduced. Second, the prey searches for a refuge area to prevent predation, which leads to the reduction of the predator's predation rate. Third, additional food is supplied to maintain the survival and prevent the oscillation or extinction of the species. Finally, in view of a time delay of the effect of the prey's fear on the growth rate of prey, the delayed model (2.2) is constructed.

The history functions of (2.2) are taken as

where $\chi_i(\vartheta)\geq0\; (i = 1, 2)$ denotes any continuous function mapping from $[-\tau, 0]$ into $R_+ = [0, \infty)$ with a sup norm $||\chi|| = \sup\limits_{[-\tau, 0]}\{|\chi_1(\vartheta)|, |\chi_2(\vartheta)|\}.$ For biological explanations of the parameters in model (2.2), refer to Table 1, where the meaning of the dimension is referred to ^[29].

3.   Preparations

3.1.   Persistence and extinction of species

We start with the positiveness of solutions of system (2.2). Let $\Psi_1(x, y) = b\left(\varsigma+\frac{\sigma(1-\varsigma)}{\sigma+y(t-\tau)}\right) -c x- \frac{1}{1+a}\frac{q(1-\kappa)y}{d+(1-\kappa) x+\gamma\delta A},$

$\Psi_2(x, y) = \frac{p((1-\kappa)x +\delta A)}{d+(1-\kappa) x +\gamma\delta A}-(e+\frac{\omega-e}{1+y}),$ and then we rewrite (2.2) as

Consequently,

which implies the positiveness of the solution of (2.2).

Next, we prove the boundedness of solutions. From (2.2), we have

and then we derive from inequality theory that $x\leq\frac{b(\varsigma+1) }{c}.$ Nevertheless, under the condition that $\frac{q}{\gamma(1+a)}-e < 0$, we can choose a sufficiently small positive number $\xi$ such that $\frac{q}{\gamma(1+a)}+\frac{\xi q}{p(1+a)}-e < 0.$ Define $M(t) = x(t)+\frac{q}{p(1+a)}y(t)$, and then we have

By the differential comparison theorem, we have

Letting $t\rightarrow \infty$ leads to $M(t)\leq\frac{\xi}{\mathfrak{M}}$.

On the other hand, by the monotonicity, we obtain from (2.2) that

If $\frac{cp\delta A}{cd+(1-k)b(\varsigma+1)+c\gamma\delta A}-\omega > 0,$ then by the comparison theorem, we have

Similarly, we have

If $b\varsigma -\frac{1}{1+a}\frac{q(1-\kappa)\mathfrak{M}}{\gamma\delta A} > 0$, then

Consequently, the boundedness of $x(t)$ and $y(t)$ are derived. That is, both prey and predator are persistent from the biological angle.

At the end of this section, we give some criteria of species extinction. From (2.2) we have $\frac{dy }{dt}\leq y(p(1+\frac{1}{\gamma})-e).$ If $p(1+\frac{1}{\gamma})\leq e,$ then $\lim\limits_{t\rightarrow \infty}y(t) = 0,$ i.e., the predator will die out. By the monotonicity, together with the boundedness of $x\leq\frac{b(\varsigma+1) }{c},$ it follows from (2.2) that

If $\frac{p((1-\kappa)\frac{b(\varsigma+1) }{c} +\delta A)} {d+(1-\kappa) \frac{b(\varsigma+1) }{c} +\gamma\delta A} < e$, then $\lim\limits_{t\rightarrow \infty}y(t) = 0$, and the predator is also extinct.

In addition, due to the positiveness of $y(t)$, letting $\Upsilon = \liminf\limits_{t\rightarrow \infty} y(t),$ we obtain from (2.2) that

If $b(\varsigma+1)- \frac{1}{1+a}\frac{cq(1-\kappa)\Upsilon}{cd+(1-\kappa) b(\varsigma+1) +c\gamma\delta A} < 0$, then $\lim\limits_{t\rightarrow \infty}x(t) = 0,$ that is, the prey is extinct.

3.2.   Existence of fixed point

We discuss the trivial equilibrium state of (2.2). It is clear that $F_0^\star(0, 0)$ is the trivial fixed point of (2.2). We compute that the boundary-value fixed points are $F_1^\star(\bar{x}, 0)$ and $F_2^\star(0, \bar{y}),$ where $\bar{x} = \frac{\omega(d+\gamma\delta A)-p\delta A}{(p-\omega)(1-\kappa)}$ and $\bar{y} = \frac{\omega(d+\gamma\delta A)-e\delta A}{c\delta A-e(d+\gamma\delta A)}$. The coexistence fixed point, say $F^\star(x^\star, y^\star)$, should satisfy the following equations:

By (3.2) we have

If $e(d+(1-\kappa)x^\star+\gamma\delta A) < p((1-\kappa)x^\star+\delta A) < w(d+(1-\kappa)x^\star+\gamma\delta A)$, then $y^\star > 0$. Putting $y^\star$ into (3.1), we have

We rewrite this as

where $H_1 = C_1D_1, H_2 = (C_1D_2+C_2D_1), H_3 = (C_1D_3+C_3D_1+C_2D_2+G_1), H_4 = (C_2D_3+C_3D_2+G_2),$ $H_5 = C_3D_3+G_3 = 0,$ and

By choosing $\omega < \sigma e < e$, $p < e$ and suitable $\varsigma$ and $A$ such that $H_i > 0(i = 1, 2, 3, 4)$ and $H_5 < 0$, then by Descartes's rule of signs, we conclude that there is a unique positive solution $x^\star$. In later discussion, we always assume that system (2.2) has a unique fixed point denoted by $F^\star(x^\star, y^\star)$.

4.   LAS and Hopf bifurcation

In this section, we investigate the local asymptotical stability (LAS, for simplicity) of solutions around the fixed point $F^\star(x^\star, y^\star)$ and analyze the existence of Hopf bifurcation. For system (2.2), we make a transformation of $X = x-x^\star, Y = y-y^\star,$ and then it can be linearized as

where

The Jacobian matrix at $F^\star(x^\star, y^\star)$ is

For simplicity, we define

and then it is not difficult to obtain that

By an easy computation, the characteristic equation $|A-\lambda I| = 0$ ($I$ is a unit matrix) reads

Next, we discuss the following two scenarios.

4.1.   No delay case

If $\tau = 0$, then Eq (4.1) becomes

At the equilibrium state $F^\star_0(0, 0)$, both characteristic roots are zero, and the stability cannot be decided directly. At $F_1^\star(\bar{x}, 0)$, there is a positive root $\lambda = c\bar{x} > 0,$ and at $F_2^\star (0, \bar{y})$, there is a root $\lambda = a_4 > 0$, so they are both unstable. Now, we investigate the LAS of the solutions of (2.2) around $F^\star(x^\star, y^\star)$.

If $a_1+a_4 > 0, a_1a_4-a_2a_3-a_3a_5 > 0$, then both the characteristic roots of (4.2) are negative, and (2.1) is LAS around $F^\star(x^\star, y^\star)$. If $a_1+a_4 < 0$ or $a_1+a_4 > 0, a_1a_4-a_2a_3-a_3a_5 < 0$, there exists a positive root at least, and system (2.1) is unstable. Consequently, we have the following finding.

Theorem 4.1. System (2.1) is LAS around the coexistence point $F^\star(x^\star, y^\star)$ provided $a_1+a_4 > 0, a_1a_4-a_2a_3-a_3a_5 > 0$ and is unstable provided $a_1+a_4 < 0$ or $a_1+a_4 > 0, a_1a_4-a_2a_3-a_3a_5 < 0$.

It is easy to obtain that $a_1a_4-a_2a_3-a_3a_5 > 0$ is equivalent to

On the other hand, let

Together with the definition of $F^\star$, we have

We take $\sigma$ as a variable and solve the equation $\wedge(\sigma) = 0,$ and then the solution is

By computation, we have

that is, the transversality condition of Hopf bifurcation holds. Thus, by the Hopf bifurcation theorem, system (2.1) undergoes a Hopf bifurcation accompanied by a limit cycle at $\sigma = \sigma^\star.$ Similarly we can take the refuge $\kappa$, the quality of additional food $\gamma$ and the level of additional food $\delta$ as bifurcation parameters. In the same manner, we obtain the critical value such that $\wedge(\kappa^\star) = 0$, $\wedge(\gamma^\star) = 0$ and $\wedge(\delta^\star) = 0$, and $\left.\frac{d\wedge}{d\sigma}\right|_{\kappa = \kappa^\star}\neq0$, $\left.\frac{d\wedge}{d\gamma}\right|_{\gamma = \gamma^\star}\neq0$, and $\left.\frac{d\wedge}{d\delta}\right|_{\delta = \delta^\star}\neq0$. Hence, they all satisfy the transversality condition, and Hopf bifurcations occur at $\kappa = \kappa^\star$ or $\gamma = \gamma^\star$ or $\delta = \delta^\star$, respectively. We summarize our conclusions in the following theorem.

Theorem 4.2. If $a_1a_4-a_2a_3-a_3a_5 > 0$ holds, then system (2.1) has a Hopf bifurcation around the fixed point $F^\star$ at the parameter thresholds of

Remark 4.1. Theorems 4.1 and 4.2 show that under some conditions, the prey and predator can coexist, and the system is stable. If the conditions fail, then they cannot steadily coexist. The prey and predator can be in a steady state of oscillatory coexistence after some parameter critical value.

4.2.   Delayed-case

For system (2.2), the equilibrium state is unchanged. The stability of $F_0^\star$, $F_1^\star$ and $F_2^\star$ are all identical with that in (2.1), and hence, we only study the stabilities of coexisting equilibrium state $F^\star(x^\star, y^\star)$.

Under the conditions of Theorem 4.1, the fixed point $F^\star(x^\star, y^\star)$ is LAS. We denote the characteristic roots of (4.1) by $\lambda(\tau) = \phi(\tau)\pm\varphi(\tau),$ then by the stability theory, the LAS of (2.1) indicates that $\phi(0) < 0$. By the continuity of $\phi(\tau)$, for some sufficiently small $\tau > 0$, we have $\phi(\tau) < 0$, and then $F^\star(x^\star, y^\star)$ retains LAS still. Now, we try to find a $\tau^\star$ such that $\phi(\tau^\star) = 0$ and $\varphi(\tau^\star)\neq0$, that is, the real part of the eigenvalue is zero, and the sign of the roots of (4.1) is changed when $\tau < \tau^\star$ and $\tau > \tau^\star$, whence the qualitative behaviors change with the varying $\tau$. Then, the coexistence steady state may lose its stability, and a Hopf bifurcation occurs.

Substituting $\lambda(\tau) = \phi(\tau)\pm\varphi(\tau)$ in (4.1) and separating the real and imaginary parts from both sides gives:

For the objective of the existence of purely imaginary roots, we take $\phi(\tau) = 0$, and then (4.4) turns to

Squaring and adding the two sides of (4.5), we get a bi-quadratic equation of $\varphi$ as follows:

Let $W = \varphi^2, \varrho_1 = (a_1^2+a_4^2+2a_2a_3), \varrho_2 = (a_1a_4-a_2a_3)^2-a_3^2a_5^2$, and then (4.5) is equivalent to the following quadratic equation:

If $\varrho_1 < 0, \varrho_2 > 0$, then Eq (4.6) contains no positive real root, and therefore there exists no real number $\phi.$ Hence for any $\tau > 0$, the coexistence equilibrium point $F^\star$ is always LAS.

Condition 1 $\varrho_2 < 0$, $2((\varphi^\star)^2+a_2a_3-a_1a_4)+(a_1+a_4)^2 > 0.$

If $\varrho_2 < 0$, there exists exactly one real root denoted by $W_1^+$ for (4.6). Then, the two real roots of (4.6) are $\varphi = \pm\sqrt{W_1^+}$. Now, we want to find a threshold value $\tau = \tau^\star$ such that (4.6) has two purely imaginary roots, i.e., $\lambda(\tau) = \phi(\tau)\pm\varphi(\tau)i = \pm i\sqrt{W_1^+}.$

From (4.5), we obtain the critical value of $\tau$,

Let $\tau^\star = \min\limits_{j}\tau_j$, and then $\tau^\star$ meets the above requirement. Next, it is only to verify that

Differentiating (4.1) w.r.t $\tau$ leads to

At $\tau = \tau^\star$, if $2((\nu^\star)^2+a_2a_3-a_1a_4)+(a_1+a_4)^2 > 0$, we get

That is, the transversality condition holds. Therefore, system (2.2) has a Hopf bifurcation at $\tau = \tau^\star$.

Under condition (2), Eq (4.6) has two positive roots denoted by $W_1^\star$ and $W_2^\star$. Let $\varphi_1^\star = \pm\sqrt{{W}_1^\star }$ be two real roots and ${\varphi}_2^\star = \pm\sqrt{{W}_2^\star }$ be another two real roots of (4.6), respectively. Then, corresponding to ${\varphi}_1^\star$ and ${\varphi}_2^\star$, we can find two threshold values ${\tau}_1^\star$ and $\tau_2^\star$ such that $\phi({\tau}_1^\star) = 0, \varphi({\tau}_1^\star) = {\varphi}_1^\star$ and $\phi({\tau}_2^\star) = 0, \varphi({\tau}_2^\star) = {\varphi}_2^\star.$ Then, the roots of the characteristic Eq (4.1) are purely imaginary $\pm i{\varphi}_1^\star$ and $\pm i{\varphi}_2^\star.$ From (4.5), we obtain

Let ${\varphi}_k^\star = \min\limits_{j}\tau_j^k, k = 1, 2$, and then they are the critical values satisfying the above requirements.

Next, we verify the transversality condition. Under Condition 2, we have

and

Then, by Hopf bifurcation theorem, the stability of $F^\star$ changes between stability and instability $r$ times, where $r$ is a positive integer. $F^\star$ is LAS when $\tau\in[0, {\tau}_0^1)\cup({\tau}_0^2, {\tau}_1^1)\cup\cdots\cup({\tau}_{r-1}^2, {\tau}_r^1)$ and unstable for all $\tau\in[{\tau}_0^1, {\tau}_0^2)\cup({\tau}_1^1, {\tau}_1^2) \cup\cdots\cup({\tau}_{r-1}^1, {\tau}_{r-1}^2)\cup({\tau}_{r}^1, \infty).$ That is, system (2.2) has multiple Hopf-bifurcations at $\tau = \tau_j^{1}$ and $\tau = \tau_j^{2}, j = 1, 2, \cdots, r,$ respectively. We summarize our findings as follows.

Theorem 4.3. We suppose $F^\star$ of (2.1) is LAS, and then for system (2.2), we have the following results.

(i) If $\varrho_1 < 0$ holds, then the solutions of (2.2) are LAS around $F^\star$ for any $\tau > 0$.

(ii) If Condition 1 holds, then the solutions of (2.2) are LAS around $F^\star$ for any $\tau < \tau^\star$ and unstable for all $\tau > \tau^\star$. System (2.2) undergoes a Hopf-bifurcation at $\tau = \tau^\star.$

(iii) If Condition 2 holds, then the solutions of (2.2) are LAS around $F^\star$ for any $\tau\in[0, {\tau}_0^1)\cup({\tau}_0^2, {\tau}_1^1)\cup\cdots\cup({\tau}_{r-1}^2, {\tau}_r^1)$ and unstable for all $\tau\in[{\tau}_0^1, {\tau}_0^2)\cup({\tau}_1^1, {\tau}_1^2) \cup\cdots\cup({\tau}_{r-1}^1, {\tau}_{r-1}^2)\cup({\tau}_{r}^1, \infty).$ System (2.2) has multiple Hopf-bifurcations at $\tau = \tau_j^{1}$ and $\tau = \tau_j^{2}, j = 1, 2, \cdots, r,$ respectively.

Remark 4.2. Theorem 4.3 indicates that the fear delay brings large influence to the asymptotical stability of the coexistence equilibrium point. When the delay is in some range, the prey and predator can coexist and keep stable, but when the delay changes, they may be in a oscillatory coexistence, and if the delay is continuously increased, then they can coexist again. Different delays will make the system change its stability, and multiple stabilities will occur, which is exhibited graphically in the next section.

5.   Examples and simulations

In this part, we numerically validate the main results and explore how the fear, refuges, additional food and the delay affect the long-term behaviors of all species. Let

By computation, the unique fixed point of (5.1) is $F^\star(1.8793, 1.5447).$ Now, we analyze the impacts of such parameters as $\varsigma, \sigma$ and anti-predation parameter $a$ induced by the prey's fear.

First, the effect of varying parameter values on the equilibrium point is listed in Table 2. From Table 2 we find that when $\varsigma$, $\sigma$ and $a$ get larger, the number of prey is smaller, but the number of predators is larger. This indicates that the fear can decrease the prey's biomass and increase the predator's biomass, that is, it is harmful to the prey species but beneficial to the predator itself.

Next, we analyze how these parameters affect the long-term behaviors of species. Keep the parameters in (5.1) unchanged and take $\omega = 0.2$, and then the fixed point is $F^\star(0.6533, 1.8678).$ By numerical simulation, we find it is stable; see Figure 1.

● Effect of fear

We change the value of fear and take $\sigma = 4$. Then the fixed point is $F^\star(0.6533, 3.0003)$, and simulation results show it is unstable (Figure 2). By Theorem 4.3, we know that there exists a Hopf bifurcation at the critical value of $\sigma$. By computation, the threshold is $\sigma^\star = 5.3945.$ See Figure 3. Similar analysis can be performed for $\varsigma$ and $a$.

● Effect of refuge and additional food

To see the effect of the refuge and additional food, we change the values of $\kappa$, $\gamma$ and $\delta$, and we keep the rest same as with (5.1). The results are visually exhibited in the following Table 3.

From Table 3 we observe that, when the refuge parameter $\kappa$ is larger, the biomass of prey is larger. The biomass of the predator is larger at first, but when $\kappa$ exceeds some value, it gets smaller. That means that the refuge area is helpful to the prey itself and does harm to the predator in some cases. When $\gamma$ is larger, that is, the quality of additional food is better, the prey biomass is larger, and the predator begins with an increasing biomass, but the number of predators decrease instead when $\gamma$ exceeds some critical value. That means appropriate additional food supplement can benefit the growth of predators and decrease the predation rate of predator on prey, which helps to preserve both the prey and predator's survival. Meanwhile, too much good quality food will have negative function to the predator's long-term behavior and as a result, the number of predators is decreased. For the effectual food parameter, Table 3 intuitively shows that $\delta$ is conducive to the growth of predators, but it is adverse to the prey's survival, despite the reduction of being predated by predator. Therefore, it is very important to afford additional food rationally and scientifically. Now, we demonstrate the long-term behaviors of the fixed point of (5.1). Figure 1 shows system (5.1) is stable. Now, we change the value of $\kappa$ and take $\kappa = 0.4$, and then it is unstable (see Figure 4). Figure 5 shows that a Hopf bifurcation occurs at $\kappa = 0.4295$. Simultaneously, we take $\gamma = 4.8$ and $\delta = 0.2$ respectively, and then simulations show that they are both unstable and undergo Hopf bifurcations at $\gamma = 4.8846$ and $\delta = 0.4053$ (see Figures 6 and 7), respectively.

● Effect of delay

Finally, we explore the effect of fear delay. Take parameters the same as in Figure 1, and then it is stable with no delay. If we take $\tau = 1$, the stability is unchanged (Figure 8(a)), but when $\tau = 4$, obviously the stability varies, accompanied by periodic solutions (Figure 8(b)). By simulations, Hopf bifurcation happens at $\tau = 3.3589$ (Figure 8(c)).

For verification of multiple bifurcation of time delay, we take the following parameter values.

By simulation, multiple stability occurs; see Figure 9(a). For system (5.2), when the fear delay $\tau\in(0, 4)$, it is stable, but when $\tau\in(4, 18)$, (5.2) becomes unstable with periodic fluctuation. With the continuous increase of delay as $\tau\in(18, 33)$, it gets stable again, and as $\tau\in(33, 50),$ (5.2) gets unstable once again. That is, the stability of (5.2) can be changed by the values of fear delay. With the continuous increase of delay, the system stability undergoes multiple changes from stable to unstable, to stable, and once again to unstable.

Take another set of parameter values,

We can see from Figure 9(b) that system (5.3) changes its stability from unstable to stable, then to unstable again. The stability of system (5.3) undergoes a similar change as appeared in Figure 9(a). These examples show multiple stability induced by time delay, which can be manipulated by choosing some suitable delay value.

Simulation results indicate that time delay brings significant influence to the system dynamics. It can make system dynamics change from stability to instability or vice versa, and it can even produce multiple stabilities. Therefore we should consider the effect of time delay and make good use of the positive function for the system development and natural balance.

6.   Conclusions

Considering the inner effect of interplay of species, such as fear of predator and anti-predation response of prey, and the effect outside like additional food supplement, we propose a delayed prey-predator model with fear, refuge and additional food. By use of Lyapunov stability theory together with the tool of the Jacobian matrix, we study the basic properties of solutions and the existence of Hopf bifurcations, and we establish some criteria. See Theorems 4.1–4.3, which are all demonstrated by examples numerically (Figures 1–9).

By theoretical study and numerical simulations, we find such factors as fear, additional food, refuge area and time delays play key roles in the system stability, which are summarized as follows.

(i) The prey's fear of the predator not only affects the prey's survival but also hurts the predator itself if the fear is too large.

(ii) Refuge area is useful to preserve the prey, while it will affect the predator's survival, so additional food should be afforded in some cases.

(iii) Additional food supplement has both positive and negative effects on the system stability, and it should be used rationally and scientifically.

(iv) Time delay of fear makes the system exhibit abundant dynamical properties like multiple stabilities and periodic limit cycles, which should be well used in the system control.

In addition, every population system is exposed in an open environment, so the random environmental disturbance should be involved ^[37]. The foraging activities happen in adult predators and the young predator has almost no predation ability, so the stage structure of the predator should be considered ^[38]. Particularly, the main results of ^[39] and ^[40] show that the bifurcation behaviors of integer-order and fractional-order models are much different, so it is necessary to investigate the bifurcation of the corresponding fractional-order models of (2.2). All these are very interesting topics, and we will study them in the near future.

Acknowledgments

This research is supported by the NSF of China (No.11861027). The author is very grateful to the anonymous referees and the editors for their careful reading and valuable comments, which have helped to improve the presentation of this work significantly.

Conflict of interest

There are no competing interests in this paper.