Influence of fear effect on a Holling type III prey-predator system with the prey refuge

1.   Introduction

The research of the dynamic relationship between prey and predator has been and will continue to be a hot topic for a long time because of its extensive existence and importance (see [1-9] and the references cited therein). There are a lot of outstanding works about the famous Lotka-Volterra type prey-predator system after it was brought up by Lotka and Volterra [10,11]. In 1959, a Canadian scholar named Holling [12] proposed the corresponding functional response function for different types of species to depict the predation rate of predator population to prey population according to his experimental results, which include three main types Holling type I, II and III, among them, Holling type III functional response function, i.e., $\frac{\alpha x^{2}}{\beta^{2}+x^{2}}$ is applicable to cattle, sheep and other vertebrates. From then on, the research on functional response of Holling type III has gradually become another important direction in the study of predator-prey dynamics (see [13-21]). In the real world, many prey species use shelters to protect themselves from being captured by predators. In order to investigate the impact of refuge on population interaction, it is necessary to establish a mathematical model of predator-prey including refuge. Many scholars have made achievements in this field [22-32]. Yunjin Huang et al. [25] proposed and investigated a prey-predator system incorporating Holling type III response function and a prey refuge, given by Eq (1.1)

where the meaning of all parameters of system (1.1) is shown in Table 1. The authors of this paper obtain the following conclusions: There is only one limit cycle in the system when the positive equilibrium is unstable; when the positive equilibrium is locally asymptotically stable, it is also globally asymptotically stable; Sufficient shelters can improve the stability of the system by eliminating periodic solutions, while less shelters will not change the dynamic stability of the system.

However, in nature, fear of predators has a variety of effects on animals, including habitat use, foraging behavior, reproduction and physiological changes. There are more and more works about the predator-prey system including fear effect in recent years, see [33-39]. Zanette et al. [40] used the playback of predator calls to control fear factors in the study of the effect of fear on free-living songbird population, and eliminated the effect of direct predation on the experiment by isolation. The research indicated that the number of offspring of sparrows would be reduced by 40$\%$ due to the fear of predators alone, and the predation risk itself was enough to affect the changes of wild animal population. In order to establish a model to simulate the impact of fear on species reduction, we use a function $F(n, y)$ to express the fear factor which is used to measure the consumption of anti-predator defense owing to the fear on the system. From the biological viewpoint and experimental results, the fear factor $F(n, y)$ should meet [33,41] $F(0, y) = 1$, $F(n, 0) = 1$, $\lim\limits_{n\rightarrow +\infty}F(n, y) = 0$, $\lim\limits_{y\rightarrow +\infty}F(n, y) = 0$, $\frac{\partial F(n, y)}{\partial n} < 0$ and $\frac{\partial F(n, y)}{\partial y} < 0$. Wang et al. [42] introduced a simple function $F(n, y) = \frac{1}{1+ny}$ as the fear factor. Here, $n > 0$ indicates the level of fear and $y$ is the predator population density at time $t$.

Up to now, no one has studied a Holling type III prey-predator system with fear effect and a prey refuge. Inspired by the above articles, we extend model (1.1) by incorporating the fear factor $F(n, y) = \frac{1}{1+ny}$ to the intrinsic growth by multiplication. Accordingly, the system $(1.1)$ becomes

The meaning of the parameters of system (1.2) is consistent with system (1.1), what's more, $n > 0$ denotes the level of fear. We can find corresponding models in the real world, such as snow leopard which is predator species and Tibetan antelope which is prey species and has been protected in protected areas. This model has a strong biological background and significance.

The rest of this paper consists of the following sections: In Section 2, we provide a qualitative analysis of the system. In Section 3, we analyze bifurcation of the system and demonstrate the occurrence of limit cycle. In Section 4, we consider the influence of fear effect and the refuge on the system. In Section 5, numerical simulation is done to verify the rationality of the conclusion. In Section 6, we finish this paper with a short discussion.

2.   Qualitative analysis

For the convenience of research, we first take the following variable substitution for system (1.2)

then, system (1.2) is simplified as

where

Taking the existence of the equilibria and practicability of system (1.2) into account, we suppose $c < k\alpha < 2c$ throughout this article. Hence, $B_{0}, \; B_{1}, \; B_{2}, \; B_{3}$ and $B_{4}$ are all positive constants.

Firstly, we provide all possible equilibria of system (2.1).

(i) The extinction equilibrium $E_{0}(0, 0)$ always exists.

(ii) There is only one semi-trivial equilibrium $E_{1}(u_{0}, 0)$.

We know $u_{0}$ is the positive real root of the following cubic equation of one variable

After verification, the equation has only one positive real root $u_{0} = \frac{B_{2}}{B_{3}}$. The other two roots are $u = \pm\frac{1}{1-m}\sqrt{\frac{k\alpha-c}{c}}i$.

(iii) There is only one coexistent (positive) equilibrium $E^{*} = (u^{*}, v^{*})$, where $u^{*} = \frac{1}{1-m}$, $v^{*}$ is the positive real root of the following equation

which is equivalent to the following quadratic equation of one variable

When $\frac{k\alpha^{2}}{c^{\frac{3}{2}}(k\alpha-c)^{\frac{3}{2}}}\Big(\frac{b}{(1-m)^{2}}\sqrt{\frac{c}{k\alpha-c}}-\frac{a}{\beta(1-m)}\Big) < 0$ namely $0\leq m < m_{1},$ where

The equation has unique positive real root

where

Secondly, the dynamic behavior of the equilibria $E_{0}, \; E_{1}$ and $E^{*}$ is discussed.

Theorem 2.1. $E_{0}(0, 0)$ is a saddle point.

Proof. The Jacobian matrix of system (2.1) at $E_{0}$ is

its eigenvalues are $\lambda_{1} = B_{0} > 0$ and $\lambda_{2} = -1 < 0$, so the extinction equilibrium $E_{0}(0, 0)$ is a saddle point.

The proof is complete.

Theorem 2.2. $E_{1}(u_{0}, 0)$ is locally asymptotically stable for $m > m_{1}$, and $E_{1}$ is unstable (saddle) for $0\leq m < m_{1}$.

Proof. The Jacobian matrix of system (2.1) at $E_{1}$ is

its eigenvalues are $\lambda_{1} = -\frac{B_{2}^{3}}{B_{3}^{2}}-B_{0} < 0$ and $\lambda_{2} = -1+(1-m)^{2}\frac{B_{2}^{2}}{B_{3}^{2}}$. Hence $E_{1}$ is locally asymptotically stable if $-1+(1-m)^{2}\frac{B_{2}^{2}}{B_{3}^{2}} < 0$, that is $m > m_{1}$, and $E_{1}$ is unstable (saddle) if $0\leq m < m_{1}$.

Remark 2.1. From Theorem 2.2, where the shelter rate $m$ takes $m_{1}$ as the threshold, we can observe that when the shelter rate exceeds $m_{1}$, the boundary equilibrium point $E_{1}$ is locally gradually stable, otherwise the boundary equilibrium point $E_{1}$ is unstable. A reasonable biological explanation is that: (1) A higher prey refuge rate $m$ is obviously beneficial to prey species, and the natural high prey refuge rate $m$ value always helps prey species obtain their biomass; (2) The high prey refuge rate leads to the excessive lack of food source of predator population, which leads to extinction.

Theorem 2.3. $E^{*}(u^{*}, v^{*})$ is locally asymptotically stable if $m_{2}\leq m < m_{1}$ or $0\leq m < m_{2}$ and $n > n_{1}$ hold, where

and $E^{*}$ is unstable if $0\leq m < m_{2}$ and $0 < n < n_{1}$ hold.

Proof. The Jacobian matrix of system (2.1) at $E^{*}$ is

where

The secular equation of matrix $J(E^{*})$ is $\lambda^{2}-A_{11}\lambda-A_{12}A_{21} = 0$, and the two eigenvalues meet $\lambda_{1}\lambda_{2} = det(J(E^{*})) = -A_{12}A_{21} > 0$, $\lambda_{1}+\lambda_{2} = tr(J(E^{*})) = A_{11}$, then the two eigenvalues are of the same sign, thus $E^{*}$ is stable when $A_{11} < 0$, and unstable when $A_{11} > 0$.

By calculation, $A_{11} < 0$ is equivalent to

If $m_{2}\leq m < m_{1}$ holds, the inequality (2.6) clearly holds because the left of the formula is positive sign. When $0\leq m < m_{2}$ holds, by substituting $v^{*}$ into inequality, after calculation, it needs to satisfy $n > n_{1}$. Hence, $E^{*}$ is locally asymptotically stable for $m_{2}\leq m < m_{1}$ or $0\leq m < m_{2}$ and $n > n_{1}$, and $E^{*}$ is unstable for $0\leq m < m_{2}$ and $0 < n < n_{1}$.

Remark 2.2. From Theorem 2.3, when the shelter rate $m$ is in the interval $[m_{2}, \; m_{1})$, the coexistence equilibrium point is locally asymptotically stable, which means that the appropriate shelter rate makes the system reach a stable state. At this time, the fear factor has no effect on the stability of the system; However, when the shelter rate is relatively small, that is, when the shelter rate is less than $m_{2}$, the stability of the system will be affected by the fear of prey population to predator. Here, we observe that when the fear caused by predator is at low level, the system shows unstable system dynamics. When the fear caused by predator is at high level, it shows a stable state. A reasonable biological explanation for this phenomenon is that when prey species are very afraid of predators, they will reduce foraging activities and adapt to different defense mechanisms to avoid predation. Fear factors greatly help predator species increase their biomass, so in the long run, it also helps the persistence of predator species and improves the stability of the whole system.

Next, the sufficient condition of global stability for the coexistent equilibrium $E^{*}(u^{*}, v^{*})$ is obtained.

Theorem 2.4. $E^{*}$ is globally asymptotically stable if $a > c$ and $\max\{m_{2},$ $m^{*}\}\leq m < m_{1}$ hold, where

Proof. From Theorem 2.3, we know for $m_{2}\leq m < m_{1}$, $E^{*}$ is locally asymptotically stable. Let

and construct a Dulac function as:

then

Set $g(u) = -2B_{3}u^{3}+(B_{2}-(1-m)^{2})u^{2}-B_{0}+1,$ and

then we get $u = 0\; or\; u = \frac{B_{2}-(1-m)^{2}}{3B_{3}}.$ If $a > c$, and $c < k\alpha < 2c$, then $a+c > k\alpha$, i.e. $\frac{B_{2}-(1-m)^{2}}{3B_{3}} > 0.$

So

When $g_{max} < 0$, i.e. $m > 1-\frac{3\sqrt{3}b\beta}{(a+c-k\alpha)}\sqrt{\frac{a-c}{a+c-k\alpha}} = m^{*}$, we get

By the Dulac Theorem in reference [43], there is no limit cycle in the positive region of u-v plane. Thus, $E^{*}$ is globally asymptotically stable if $a > c$ and $\max\{m_{2}, m^{*}\}\leq m < m_{1}$ hold.

We use a table to show the existence and stability of all equilibria of system (2.1), as shown in Table 2.

3.   Bifurcation analysis

In this part, we analyze bifurcation of the system and demonstrate the occurrence of limit cycle.

Theorem 3.1. Suppose $0\leq m < m_{2}$, then system (2.1) goes through a Hopf bifurcation around $E^{*}$ at $n = n_{1}$.

Proof. The secular equation of matrix $J(E^{*})$ is $\lambda^{2}-tr(J(E^{*}))\lambda+det(J(E^{*})) = 0$, and $det(J(E^{*})) = -A_{12}A_{21} > 0,$ $tr(J(E^{*})) = A_{11}$, then

Obviously,

Hence, system (2.1) goes through a Hopf bifurcation around $E^{*}$ at $n = n_{1}$.

Theorem 3.2. Suppose $0\leq m < m_{2}$ and $0 < n < n_{1}$, system (2.1) has one limit cycle.

Proof. From Theorem 2.2 and Theorem 2.3, $E^{*}$ is a unstable point and $E_{0}$, $E_{1}$ are saddle points, when $0\leq m < m_{2}$ and $0 < n < n_{1}$ hold. Suppose that

where $u_{B}, v_{D}, v_{E}$ satisfy $u_{B} > u_{0} > u^{*}, \max\limits_{u^{*} < u < u_{B}}\{u(-B_{3}u^{3}+B_{2}u^{2}-B_{1}u+B_{0}+1), v^{*}\} < v_{E} < v_{D},$ and $D(u^{*}, v_{D}), E(u_{B}, v_{E})$ are shown in Figure 1.

Therefore

By the Poincar$\acute{\rm{e}}$-Bendixson theorem [44], system (2.1) has one limit cycle in the domain I as shown in Figure 1. The proof is complete.

From Theorem 3.2, we know that system (2.1) has one limit cycle in the first quadrant, and from numerical simulation results it is possible to observe that there is a unique limit cycle.

4.   The influence of fear effect

In this part, we will study the impact of fear effect on the system.

$\bullet$ The influence of the fear factor on the stability of system (2.1)

From Theorem 2.3 and Theorem 2.4, the coexistent equilibrium $E^{*}(u^{*}, v^{*})$ is locally asymptotically stable if $m_{2}\leq m < m_{1}$. In such instance, regardless of the level of fear, the stability of the system will not be affected. If $0\leq m < m_{2}$, as the level of fear increases, $E^{*}(u^{*}, v^{*})$ changes from unstable state to stable one, and $n = n_{1}$ is the critical value. At this time, the fear factor can stabilize the system by eliminating periodic solutions. The domains of the stability of $E^{*}$ are shown in Figure 2.

$\bullet$ The influence of fear factor on predators

Because the final prey population has nothing to do with fear level $n$, we just talk about the influence of fear factor on predator species. By the calculation, it's easy to draw a conclusion that the predator population $v^{*}$ decreases with the increase of fear level $n$, since $v^{*}$ is a continuous function of $n$. Finding the derivative of $n$ on both sides of the following formula:

we get

$\bullet$ The comprehensive influence of fear factor and prey shelters on predator-prey species

In order to seek out the comprehensive influence of fear factor and prey shelters on predator-prey species, let's first consider the influence of prey shelters on preys and predators without fear factor, i.e., letting $n = 0$ in system (2.1), we have

then $u^{*}$ and $v^{*}$ are derived from $m$ respectively:

where

Hence, we know that the increase of $m$ can increase prey population, if $0\leq m < m_{1}$ holds. For $v^{*}$: if $0\leq m < m_{3}$, i.e. $\frac{dv^{*}}{dm} > 0$, then the increase of $m$ can increase predator population, while if $m_{3} < m < m_{1}$, i.e. $\frac{dv^{*}}{dm} < 0$, then the increase of $m$ can decrease predator population; When $m = m_{3}$, the predator population $v^{*}$ achieves the maximum value, and when $m = m_{1}$ i.e. $v^{*} = 0$, the predators dies out.

Next we will study the influence of prey shelters with fear factor on predator-prey species. i.e. $n > 0$.

From (2.3), the derivatives along $u^{*}$ and $v^{*}$ with respect to $m$ are

When $0\leq m < m_{1}$, $u^{*}$ is strictly monotone increasing with respect to $m$, which completely coincides with the system without fear effect.

On the other hand, by the calculation, $\frac{dv^{*}}{dm} < 0$ is equivalent to

When $m_{3} < m < m_{1}$ the inequality above clearly holds since the right side of the formula is positive sign. When $0\leq m < m_{3}$ holds, by an equivalent deformation, it needs to meet

When $n = f(m)$ holds, $\frac{dv^{*}}{dm} = 0.$ Here, the function $n = f(m)$ satisfies

and we set $n_{2} = f(0) = \frac{a(k\alpha-c)^{\frac{3}{2}}m_{3}}{2kb^{2}c^{\frac{1}{2}}\beta^{3}}.$

Hence, as shown in Figure 3, we know that $n > n_{2}$ and $0\leq m < m_{1}$ hold, $v^{*}$ is strictly monotone decreasing with regard to $m$, that means the increase of $m$ can decrease predator population, and predator population gets its maximum value at $m = 0$, i.e. without refuge; When $0 < n < n_{2}$ and $0\leq m < m_{1}$ hold, $v^{*}$ is strictly monotone increasing with respect to $m$ in the interval $[0, f^{-1}(n)]$ and decreasing with respect to $m$ in the interval $[f^{-1}(n), m_{1}]$, then predator population reaches its maximum value at $m = f^{-1}(n)$, which is influenced by fear effect and different from the situation without fear factor; When $m = m_{1}$, i.e. $v^{*} = 0$, the predator species dies out, which is similar to the situation without fear factor.

5.   Numerical results

In this part, the numerical simulations are done to further verify the validity of the above conclusions. Let's set the following parameters as:

Under these set of parameters, we get:

The simulation results are shown as follow:

In Figure 4:$\; m = 0.1, \; n = 0.01$. By a calculation, we have $n_{1} = 0.0146$. Here $0\leq m < m_{2}, \; 0 < n < n_{1}$, the system (2.1) has two saddle points: a extinction equilibrium $E_{0} = (0, 0)$, and a boundary equilibrium point $E_{1} = (8.1650, 0)$; the only unstable coexist equilibrium point $E^{*} = (1.1111, 33.4783)$. There is a limit cycle in the system, we can clearly observe that the trajectories of an initial value inside and outside the limit cycle approach the limit cycle.

In Figure 5:$\; m = 0.1, \; n = 1$, then $0\leq m < m_{2}, \; n > n_{1}$. The system (2.1) also has two saddle points: $E_{0} = (0, 0)$, $E_{1} = (8.1650, 0)$; and a unique coexist equilibrium point $E^{*} = (1.1111, 5.7810)$, which is a locally asymptotically stable spiral source point. Compared with Figure 4, we know that increase the fear level will decrease the final number of predators $v^{*}$, and change $E^{*}$ from unstable point to stable one. At this time, the fear factor can eliminate the limit cycle oscillation and enhance the stability of the system.

In Figure 6:$\; m = 0.1$, then $n_{1} = 0.0146$. The system (2.1) goes through a Hopf bifurcation around $E^{*}$ at $n = 0.0146$.

In Figure 7:$\; m = 0.5, \; n = 1$, then $m_{2}\leq m < m_{1}$. The system (2.1) also has two saddle points: $E_{0} = (0, 0)$, $E_{1} = (8.1650, 0)$; the only local asymptotical stable coexist equilibrium point $E^{*} = (2, 4.7256)$. Compared with Figures 4 and 5, increase the fear effect can decrease the final number of predators $v^{*}$, but not alter the stability of $E^{*}$.

In Figure 8:$\; m = 0.8, \; n = 1$, then $a > c$ and $\max\{m_{2}, m^{*}\}\leq m < m_{1}$. The system (2.1) has two saddle points: $E_{0} = (0, 0)$, $E_{1} = (8.1650, 0)$; and a unique coexist equilibrium point $E^{*} = (5, 1.2595)$, which is a globally asymptotically stable point. Compared with Figures 4, 5 and 7, we know that increase refuge can increase the number of preys $u^{*}$.

In Figure 9:$\; m = 0.8, \; n = 10$, then $a > c$ and $\max\{m_{2}, m^{*}\}\leq m < m_{1}$. The system (2.1) also has two saddle points: $E_{0} = (0, 0)$, $E_{1} = (8.1650, 0)$; and a unique coexist equilibrium point $E^{*} = (5, 0.12887)$, which is also a globally asymptotically stable point. Compared with Figure 8, the $v^{*}$ is further reduced to near zero by the increase of fear factor $n$, but $v^{*}$ is always greater than zero. In such instance, the fear effect does not cause the extinction of the predator population.

In Figure 10:$\; m = 0.9, \; n = 1$, then $m > m_{1}$. The system (2.1) has one extinction equilibrium $E_{0} = (0, 0)$, which is a saddle point, and one boundary equilibrium point $E_{1} = (8.1650, 0)$, which is locally asymptotically stable, and no coexist equilibrium point. In such instance, there are enough prey refuges to cause the extinction of the predator population.

In Figure 11:$\; m = 0, \; n = 1$, then $0\leq m < m_{2}, \; n > n_{1}$. At this time, the system (2.1) has no refuge. Similar to Figure 7, the system also has two saddle points: $E_{0} = (0, 0)$, $E_{1} = (8.1650, 0)$; and only one coexist equilibrium point $E^{*} = (1, 5.7764)$, which is a locally asymptotically stable point. Obviously, in the system (2.1) let $m = 0$ we can get a new system that only contains fear factor and no prey shelter, and the dynamic properties of the new system can be easily obtained from the dynamic properties of system (2.1).

In Figure 12:$\; m = 0, \; n = 0$. At this time, the system (2.1) becomes a system that has neither a refuge for the prey population nor a fear effect, which is the same as the system in [45]. Similar to Figure 4 the system has two saddle points: a extinction equilibrium $E_{0} = (0, 0)$, and a boundary equilibrium point $E_{1} = (8.1650, 0)$; the only unstable coexist equilibrium point $E^{*} = (1, 36.5636)$ and there is a limit cycle in the system.

In Figure 13:$\; m = 0.9, \; n = 0$. At this time, the system (2.1) becomes a system without fear effect which is the same as the system in [25]. Similar to Figure 10, the system only has one extinction equilibrium $E_{0} = (0, 0)$, which is a saddle point, and one boundary equilibrium point $E_{1} = (8.1650, 0)$, which is locally asymptotically stable. and no coexist equilibrium point. Compared with Figure 10, we know that the extinction of predator species has nothing to do with fear effect.

6.   Conclusions

In this paper, the influence of anti-predator behavior caused by fear of predators in a prey-predator system with Holling type III response function and prey refuge is considered. We analyze the dynamic behavior of the system mathematically, including the stability of the system and the occurrence of Hopf bifurcation around the positive equilibrium point and the existence of limit cycle emerging through Hopf bifurcation. We discover that the fear effect can stabilize the system by eliminating periodic solutions and decrease the final number of predator species at the coexist equilibrium, but not cause the extinction of predators, which is different from the system without fear factor. We also discover that prey shelters has vital role on the permanence of the predators. When $n > n_{2}$ and $0\leq m < m_{1}$ hold, the increase of the quantity of shelters can decrease predator population, and the final number of predator species reaches its maximum value without prey refuge namely $m = 0$; when $0 < n < n_{2}$ and $0\leq m < m_{1}$ hold, $v^{*}$ increases monotonically at first and then decreases monotonically in the interval $[0, \; m_{1})$ with respect to $m$, then predator species reaches its maximum value at $m = f^{-1}(n)$, which is influenced by fear effect and different from the situation without fear effect; when $m = m_{1}$ the predator species dies out, which is similar to the situation without fear effect. The system in this paper has complex dynamic behavior, which enrich the dynamic behavior of predator-prey system. From the real world, we can protect endangered animals and achieve ecosystem balance by setting up appropriate reserves.

Acknowledgments

All of the authors would like to thank Professor Zhengce Zhang for his help in revising the paper, which undoubtedly improved the paper. We are also very grateful to the editors and the anonymous reviewers for their helpful comments and constructive suggestions. We thanks the National Science Foundation of China (No. 11801238) and the horizontal research projects: Study on mathematical modeling and integrated control of diseases and insect pests in Camellia oleifera plantation (No. 204302500023) for supporting our research.

Conflict of interest

All authors declare no conflicts of interest in this paper.