Combining impact of velocity, fear and refuge for the predator–prey dynamics

Abstract

We develop a deterministic predator–prey compartmental model to investigate the impact of their velocities on their interactions. Prey hides in a refuge area and comes out of this area when predation pressure declines. To avoid predation, prey can limit their velocity. For antipredator behaviour, we examined that prey mortality increases when either predator or prey velocity increases while raising antipredator behaviour increases prey density. We proved that predator free equilibrium is globally asymptotically stable and co-existing equilibrium will be globally stable under certain conditions. We find that transcritical bifurcations occur at predator-free equilibrium at the certain value of the death rate of the predator.

Keywords:

• Predator–prey
• refuge
• fear
• velocities
• transcritical bifurcation

1. Introduction

In theory, ecology is the scientific study and aesthetic analysis of interactions among organisms and their environments. It consists of organisms, the communities they create and the nonliving components of their environment that are in dynamic interaction with each other. Historically, predators and their prey have been and will remain an important theme in ecology and mathematical ecology due to their universality. The interaction between them is a key feature of population dynamics.

A prey in the natural world protects itself from predators by hiding in an area where the predator cannot easily discover it. A refuge area is the name given to this location. By using the refuge, some members of the prey population are protected from predators, resulting in a lower predation rate. As a result, refuges are proportionate to the prey density while these are inversely proportional to the predator density. Prey refuge rate has been used as a crucial parameter in a prey–predator model in several studies [7,19,35,36]. Das and Samanta [7] investigated the impact of predator fear on prey species and where both species become infected due to environmental toxicant. A diffusive predator–prey model with prey refuge was examined by [1]. Mangroves, according to Glazner et al. [11] can change predator–prey interactions by increasing the value of prey refuge in a mangrove-marsh ecotone. A stage-organized predator–prey model with prey refuge was examined for uniform persistence and global asymptotic stability by [13]. The global stability of a stage structure predator–prey system with prey refuge was studied in [33]. In [33], the authors found that the model undergoes a Hopf bifurcation when the delay crosses specific critical values. Xie et al. [34] investigated the persistence and stability of a modified Leslie–Gower predator–prey model with prey refuge, and they found that the prey refuge had a favourable effect on the persistence feature. Some mathematical models of biological systems have been developed to capture the effect of prey refuge, which protects some of the prey from predators [2,12].

There is a tendency for predator studies to concentrate only on direct killing by predators since such deaths can easily be observed in nature. However, prey species respond to predation risk in a variety of ways, including habitat change, foraging, and alertness, and different physiological changes. Based on numerous experiments, Zanette et al. [37] have demonstrated that song sparrows reduce by $40\mathrm{\%}$ the offspring by predation fear apart from direct killing of birds. However, many field data show that the indirect impact of predator species on prey species has a major impact on population dynamics [3,4,15,25]. Das and Samanta [8] investigated the impact of predator fear on prey species and where both species become infected due to environmental toxicant. In [38], the authors explored a Holling-II predator–prey model with a prey refuge and the fear effect. They discovered that in a positive equilibrium, the fear effect not only reduces the predator's population density, but also stabilizes the system. The stability of a delayed predator–prey model with prey refuge and fear effect was investigated by [14], and it was discovered that refuge below the threshold level is beneficial to the system. Impact of fear in a delay-induced predator–prey system with intraspecific competition within predator species was studied by Das et al. [5]. Xiao and Li [32] looked at a predator–prey model with mutual interference and a fear effect. They concluded that mutual interference can stabilize the predator–prey system when compared to the corresponding predator–prey model without mutual interference. Sasmal and Takeuchi [26] investigated the multistability and Hopf bifurcation of a predator–prey system with fear effect. Das and Samanta [6] examined a stochastic model where fear is in a prey population and an alternative food is available for predators. A comparison study of predator–prey system in deterministic and stochastic environments influenced by fear and its carry-over effects was examined by [17]. In this regard, [9,16,18,20,21,29,30] have more information.

Antipredator sensitivity refers to a shift in a consumer's behaviour in response to the presence of a predator. To avoid predators, prey can limit their speed, but this reduces their capacity to absorb resources, resulting in a trade-off between hiding and foraging. Sadowski and Grosholz [24] created a basic tritrophic food chain that included predator and prey velocities. Antipredator behaviour, according to their hypothesis, allows prey to adopt to predators by reducing their speed in response to predator density. They studied the impact of both fast-moving (mobile) predators and slow-moving (sit and wait) predators on prey density and predator–prey cycle amplitude in equilibrium. Active search is advantageous when both predator and prey move randomly, but the relative advantage of active search reduces when the predator and prey adopt more directional movement [27]. In [10], the authors found that active searching predators are anticipated to prey on slow-moving species, while ambush predators are predicted to prey on fast-moving animals in a search mode model based on search in a three-dimensional space. Ross and Winterhalder [23] recently created a two-dimensional predator–prey encounter model that includes random prey and predator movement.

In a grassland of small grasses with trees in small numbers, commonly found predators are lions, tigers, etc. They feed upon other animals like zebras, deer, etc. It is very easy for predators to locate their prey in small grassland. Prey animals run fast to reach a safe place and escape from predators. To protect from predators the velocity of the weaker animals helps them to escape from predators. In this way due to their velocity they protect themselves and their chances of survival increase. Some recent studies verified it experimentally and revealed that fast velocity is unrelated to survival. Slow prey has higher chances to escape from predators in a careful and skilful way. Studies found that escape success depends on both velocity and manoeuvrability [24,27]. Ross and Winterhalder [23] examined that the encounter rate of prey and predator increases on increasing prey velocity. Slow prey has better chance to escape predation with high manoeuvrability. Large size prey needs to be faster and manoeuvrable to escape predation than the smaller size prey. It is verified that zebra and impala evade capture if they turn quickly which is possible at slow speed. If prey run very fast, then there are chance for prey to slip and fall and be easily captured by predators. Sadowski and Grosholz [24] also examined and found that animals should run slower with manoeuvrability and fastest velocity is not always best as when animals run fast their control and manoeuvrability decline. Preys with lower velocity than their predators are always in an advantageous position because they can rapidly change directions [27]. The fast running animals are likely to lose controls on their movements. High velocity reduces control and increases the chances of mistakes that goes in favour of predators. Preys of small size like Dick-dik are vulnerable by a greater range of predators. They are unable to defend themselves and to out run predators. They do not move far away their territory and make small feeding groups. They move slowly and feed cautiously. Predators detect prey only when they change their habitats. Prey prefers to hide in one place and remain still to detect ambush predators [24].

We investigate our model by including predator fear, which is based on [28] paradigm. Preys leave the refuge only when they are no longer afraid, otherwise, they return to the safe region. By enabling predator density to alter prey velocity, which specifies how prey behaviour changes in response to predator increases, we included predator and prey velocities with antipredator behaviour in our model as studied in [24]. Our present study explore the combine effects of fear, refuge and velocities on the predator–prey dynamics.

2. Model formulation

We consider three-species model, a predator species (y) and two different habitats of prey species. One is called the prey in the refuge habitat $(x_1)$ and the others are called the prey outside the refuge habitat $(x_2)$. It is assumed that the prey species are available inside the refuge habitat and that their population grows logistically, where $a$ is the growth rate of the prey species inside the refuge, $k$ is the carrying capacity of the prey in the refuge habitat. It is also assumed that the prey species inside the refuge habitat is saved from predation. The prey species moves to a second habitat outside the refuge when the pressure of predation fear is decreased, and in this habitat the prey species can be killed by predators under the law of mass action. μ is the fear parameter, α is the migration rate from the refuge habitat. As predation fear increases in the prey species due to the presence of predators, the prey species migrate to the refuge habitat, where β represents the immigration rate into the refuge habitat and b is the feeding rate of the predator on the prey outside the refuge habitat. c is the conversion rate of prey to predator. In the absence of a prey species, predators die exponentially, where d is the death rate of the predator. We included predator and prey velocities in this model, with w representing the predator velocity. The prey velocity is a decreasing function of y, is $\gamma e^{-\theta y}$. The maximum prey movement rate, defined as the velocity of the prey in the absence of a predator, is specified by $\gamma$. Preys are hampered by increasing their encounter rate with predators as the velocity rises. As a result, we include the antipredator behaviour θ in our model by allowing predator density to affect the prey velocity. Antipredator sensitivity describes how prey behaviour changes in response to predator increases. The following system of differential equations is used to simulate the predator–prey interaction: (1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1)

3. Positivity of solutions

We will prove that all solutions of the model with positive initial data will remain positive for all time t>0. $\begin{array}{rl}\frac{\mathrm{d}{x}_1}{\mathrm{d}t}{\mid }_{x_1=0}=& \beta x_2\ge 0\\ \frac{\mathrm{d}{x}_2}{\mathrm{d}t}{\mid }_{x_2=0}=& \frac{\alpha x_1}{1+\mu y}\ge 0\\ \frac{\mathrm{d}y}{\mathrm{d}t}{\mid }_{y=0}=& 0\ge 0\end{array}$ Hence, all solutions will remain positive for all time [28].

4. Boundedness

Theorem 4.1

All solutions $(x_1(t),x_2(t),y(t))$ of the system 1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) with initial conditions $x_1(0)\ge 0$,$x_2(0)\ge 0$, $y(0)\ge 0$ are bounded for all t>0.

Proof.

Let $w=x_1+x_2+y$ $\begin{array}{rl}\frac{\mathrm{d}w}{\mathrm{d}t}& =\frac{\mathrm{d}{x}_1}{\mathrm{d}t}+\frac{\mathrm{d}{x}_2}{\mathrm{d}t}+\frac{\mathrm{d}y}{\mathrm{d}t}\\ & =a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & +\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & +c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\\ & =a{x}_1(1-\frac{x_1}{k})-(b-c)x_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}$ From system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ), the feeding rate of predator on the prey outside the refuge habitat will always be greater than the conversion rate of prey to predator. Therefore, b>c and so $\begin{array}{rl}& \frac{\mathrm{d}w}{\mathrm{d}t}\le a{x}_1(1-\frac{x_1}{k})-\mathrm{d}y\\ & \frac{\mathrm{d}w}{\mathrm{d}t}\le a{x}_1-\mathrm{d}y\end{array}$ For any positive q, we have $\frac{\mathrm{d}w}{\mathrm{d}t}+qw\le a{x}_1+q{x}_1+q{x}_2+(q-d)y$ Since the carrying capacity k is the maximum population size of prey in the refuge, so $x_1<k$ and $x_2<k$ where $x_2$ is a part of $x_1$. Therefore, $\frac{\mathrm{d}w}{\mathrm{d}t}+qw\le (a+2q)k+(q-d)y$ If $q<d,$ then $\frac{dw}{dt}+qw\le (a+2q)k=L\mathbf{(}\mathbf{s}\mathbf{a}\mathbf{y}\mathbf{)}$ The above differential inequality can be expressed by $\frac{\mathrm{d}}{\mathrm{d}t}(w-\frac{L}{q})\le -q(w-\frac{L}{q})$ Using theory of differential inequality for $w(t)$, we get $0\le w(x_1,x_2,y)\le \frac{L}{q}(1-e^{-qt})+w(0)e^{-qt}$ where $w(0)=w(x_1(0),x_2(0),y(0))$. Moreover, for $t⟶\mathrm{\infty }$, $0<w\le \frac{L}{q}$. Hence all the solutions for the system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ), (initiating in ${\mathbb{R}}_{+}^3$) are confined in the region

$\mathrm{\Delta }=\{(x_1,x_2,y)\in {\mathbb{R}}_{+}^3:w(t)\le \frac{L}{q}+ϵ,\mathit{f}\mathit{o}\mathit{r}\mathit{a}\mathit{n}\mathit{y}ϵ\mathit{>}\mathit{0}\}$. This shows that the solution of the system represented by (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ) is bounded.

5. Equilibrium points of the system

We find three biologically meaningful equilibrium points:

• The extinction of all populations ${\overline{E}}_0=(0,0,0)$.

• Predator Free Equilibrium ${\overline{E}}_1=({\overline{x}}_1,{\overline{x}}_2,0)=(k,\frac{\alpha k}{\beta },0)$, where the prey species only survive and the predator goes to extinction.

• Co-existing equilibrium ${\overline{E}}_2=(x_1^{\ast },x_2^{\ast },y^{\ast })$ where

  $x_2^{\ast }=\frac{d}{c\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}}}$, $x_1^{\ast }=(\frac{1+\mu y^{\ast }}{\alpha })x_2^{\ast }(\beta +b{y}^{\ast }\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}})$ and $y^{\ast }$ is the positive root of the following equation:

(2) $\begin{array}{rl}& \left(\frac{1+\mu y^{\ast }}{\alpha }\right)[\beta +b{y}^{\ast }\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}}][a-\frac{a}{k}\left(\frac{1+\mu y^{\ast }}{\alpha }\right)\left(\frac{d}{c\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}}}\right)\\ & \times \phantom{\left(\frac{d}{c\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}}}\right)}(\beta +b{y}^{\ast }\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}})-\frac{\alpha }{1+\mu y^{\ast }}]+\beta =0\end{array}$(2)

6. Stability analysis

6.1. Local stability analysis

Theorem 6.1

  

1. The trivial equilibrium point ${\overline{E}}_0=(0,0,0)$ is always unstable.

2. Predator free equilibrium ${\overline{E}}_1=(k,\frac{\alpha k}{\beta },0)$ is locally asymptotically stable if $\begin{array}{r}c\overline{x_2}\sqrt{w^2+{\gamma }^2}<d\end{array}$

3. The co-existing equilibrium ${\overline{E}}_2=(x_1^{\ast },x_2^{\ast },y^{\ast })$ will be locally asymptotically stable if $\begin{array}{rl}& (\frac{a{x}_1^{\ast }}{k}+\frac{\beta x_2^{\ast }}{x_1^{\ast }})(\frac{\alpha x_1^{\ast }\mu }{(1+\mu y^{\ast }{)}^2}+\frac{b{x}_2^{\ast }}{\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}}}[w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}(1-\theta y^{\ast })])\\ & >\frac{{\alpha }^2{x}_1^{\ast }\mu }{(1+\mu y^{\ast }{)}^3}\end{array}$

Proof.

The stability matrix for the system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ) linearized about the equilibrium point is (3) $\begin{array}{r}\left(\begin{array}{ccc}a-\frac{2a{x}_1}{k}-\frac{\alpha }{1+\mu y}& \beta & \frac{\alpha x_1\mu }{(1+\mu y{)}^2}\\ \frac{\alpha }{1+\mu y}& -\beta -by\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}& \begin{array}{l}-\frac{\alpha x_1\mu }{(1+\mu y{)}^2}-\frac{b{x}_2}{\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}}\\ \times [w^2+{\gamma }^2{e}^{-2\theta y}(1-\theta y)]\end{array}\\ 0& cy\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}& \begin{array}{l}-d+\frac{c{x}_2}{\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}}\\ \times [w^2+{\gamma }^2{e}^{-2\theta y}(1-\theta y)]\end{array}\end{array}\right)\end{array}$(3)

1. For ${\overline{E}}_0=(0,0,0)$, the stability matrix is $\begin{array}{r}\left(\begin{array}{ccc}a-\alpha -\lambda & \beta & 0\\ \alpha & -\beta -\lambda & 0\\ 0& 0& -\mathrm{d}-\lambda \end{array}\right)\end{array}$ The corresponding characteristic equation is $(d+\lambda )[{\lambda }^2+\lambda (\beta -a+\alpha )-a\beta ]=0$. Therefore, ${\overline{E}}_0=(0,0,0)$ is unstable

2. Now the stability for ${\overline{E}}_1=(k,\frac{\alpha k}{\beta },0)$ is $\begin{array}{r}\left(\begin{array}{ccc}-(a+\alpha )-\lambda & \beta & \alpha \mu k\\ \alpha & -\beta -\lambda & -\alpha {\overline{x}}_1\mu -b{\overline{x}}_2\sqrt{w^2+{\gamma }^2}\\ 0& 0& (c{x}_2\sqrt{w^2+{\gamma }^2}-d)-\lambda \end{array}\right)\end{array}$ The corresponding characteristic equation is $((c{x}_2\sqrt{w^2+{\gamma }^2}-d)-\lambda )[{\lambda }^2+\lambda (a+\alpha +\beta )+a\beta ]=0$. Hence, the equilibrium will be locally asymptotically stable if $c{x}_2\sqrt{w^2+{\gamma }^2}<d$

3. The stability for ${\overline{E}}_2$ is $\begin{array}{r}\left(\begin{array}{ccc}{p}_1-\lambda & p_2& p_3\\ q_1& q_2-\lambda & q_3\\ 0& {\gamma }_2& {\gamma }_3-\lambda \end{array}\right)\end{array}$ where $p_1=\frac{-a{x}_1^{\ast }}{k}-\frac{\beta x_2^{\ast }}{x_1^{\ast }}<0$

   $p_2=\beta >0$

   $p_3=\frac{\alpha x_1^{\ast }\mu }{(1+\mu y^{\ast }{)}^2}>0$

   $q_1=\frac{\alpha }{1+\mu y^{\ast }}>0$

   $q_2=\frac{-\alpha x_1^{\ast }}{x_2^{\ast }(1+\mu y^{\ast })}<0$

   $q_3=-\frac{\alpha x_1^{\ast }\mu }{(1+\mu y^{\ast }{)}^2}-\frac{b{x}_2^{\ast }}{\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}}}[w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}(1-\theta y^{\ast })]<0$ if $\theta <\frac{1}{y}$

   ${\gamma }_2=\frac{y^{\ast }d}{x_2^{\ast }}>0$

   ${\gamma }_3=-d{\gamma }^2\theta y^{\ast }{e}^{-2\theta y^{\ast }}<0$

   The corresponding characteristic equation is (4) $\begin{array}{rl}& {\lambda }^3+{\lambda }^2(-q_2-{\gamma }_3-p_1)+\lambda ((q_2+{\gamma }_3)p_1+q_2{\gamma }_3-{\gamma }_2{q}_3-p_2{q}_1)\\ & +(-p_1{q}_2{\gamma }_3+p_1{\gamma }_2{q}_3+p_2{q}_1{\gamma }_3-p_3{q}_1{\gamma }_2)=0\end{array}$(4) Equation (4(4) $\begin{array}{rl}& {\lambda }^3+{\lambda }^2(-q_2-{\gamma }_3-p_1)+\lambda ((q_2+{\gamma }_3)p_1+q_2{\gamma }_3-{\gamma }_2{q}_3-p_2{q}_1)\\ & +(-p_1{q}_2{\gamma }_3+p_1{\gamma }_2{q}_3+p_2{q}_1{\gamma }_3-p_3{q}_1{\gamma }_2)=0\end{array}$(4) ) can be expressed as (5) $\begin{array}{r}{\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3=0\end{array}$(5) where $a_1=-(q_2+{\gamma }_3+p_1)>0$

   $a_2=(q_2+{\gamma }_3)p_1+q_2{\gamma }_3-{\gamma }_2{q}_3-p_2{q}_1$

   $a_3=(-p_1{q}_2{\gamma }_3+p_1{\gamma }_2{q}_3+p_2{q}_1{\gamma }_3-p_3{q}_1{\gamma }_2)={\gamma }_3(-p_1{q}_2+p_2{q}_1)+{\gamma }_2(p_1{q}_3-p_3{q}_1)$

   where ${\gamma }_3<0$, $-p_1{q}_2+p_2{q}_1<0$ if $p_2{q}_1<p_1{q}_2$

   or $\frac{\beta \alpha }{1+\mu y^{\ast }}<\frac{a\alpha x_1^{\ast 2}}{k{x}_2^{\ast }(1+\mu y^{\ast })}+\frac{\beta \alpha }{1+\mu y^{\ast }}$, $\frac{a\alpha x_1^{\ast 2}}{k{x}_2^{\ast }(1+\mu y^{\ast })}>0$ which is true that ${\gamma }_3(-p_1{q}_2+p_2{q}_1)>0$,

   For ${\gamma }_2(p_1{q}_3-p_3{q}_1)$ where ${\gamma }_2>0$, if $p_1{q}_3>p_3{q}_1$ then $a_3>0$

   or $(\frac{a{x}_1^{\ast }}{k}+\frac{\beta x_2^{\ast }}{x_1^{\ast }})(\frac{\alpha x_1^{\ast }\mu }{(1+\mu y^{\ast }{)}^2}+\frac{b{x}_2^{\ast }}{\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}}}[w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}(1-\theta y^{\ast })])>\frac{{\alpha }^2{x}_1^{\ast }\mu }{(1+\mu y^{\ast }{)}^3}$

   Now, $a_1{a}_2>a_3\Rightarrow (-q_2-{\gamma }_3-p_1)[(q_2+{\gamma }_3)p_1+q_2{\gamma }_3-{\gamma }_2{q}_3-p_2{q}_1]>-p_1{q}_2{\gamma }_3+p_1{\gamma }_2{q}_3+p_2{q}_1{\gamma }_3-p_3{q}_1{\gamma }_2$ or (6) $\begin{array}{rl}& p_1{q}_2^2+2p_1{q}_2{\gamma }_3+q_2^2{\gamma }_3-q_2{\gamma }_2{q}_3-p_2{q}_1{q}_2+p_1{\gamma }_3^2+q_2{\gamma }_3^2-{\gamma }_2{q}_3{\gamma }_3+q_2{p}_1^2+{\gamma }_3{p}_1^2\\ & -p_2{p}_1{q}_1-p_3{q}_1{\gamma }_2<0\end{array}$(6) sign of various terms of Equation (6(6) $\begin{array}{rl}& p_1{q}_2^2+2p_1{q}_2{\gamma }_3+q_2^2{\gamma }_3-q_2{\gamma }_2{q}_3-p_2{q}_1{q}_2+p_1{\gamma }_3^2+q_2{\gamma }_3^2-{\gamma }_2{q}_3{\gamma }_3+q_2{p}_1^2+{\gamma }_3{p}_1^2\\ & -p_2{p}_1{q}_1-p_3{q}_1{\gamma }_2<0\end{array}$(6) ) are mentioned

   Now, $p_1{q}_2^2-p_2{p}_1{q}_1+q_2{p}_1^2-p_2{q}_1{q}_2=q_2(p_1{q}_2-p_2{q}_1)+p_1(p_1{q}_2-p_2{q}_1)=(p_1{q}_2-p_2{q}_1)(q_2+p_1)$

   we proved $p_1{q}_2-p_2{q}_1>0$ and $q_2+p_1<0$. So, $a_1{a}_2>a_3$. Hence non-zero equilibrium will be asymptotically stable if $p_1{q}_3>p_3{q}_1$, or $\begin{array}{rl}& (\frac{a{x}_1^{\ast }}{k}+\frac{\beta x_2^{\ast }}{x_1^{\ast }})(\frac{\alpha x_1^{\ast }\mu }{(1+\mu y^{\ast }{)}^2}+\frac{b{x}_2^{\ast }}{\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}}}[w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}(1-\theta y^{\ast })])\\ & >\frac{{\alpha }^2{x}_1^{\ast }\mu }{(1+\mu y^{\ast }{)}^3}\end{array}$

Theorem 6.2

If the equilibrium point $y^{\ast }$ exists, be the positive real root of the Equation (2(2) $\begin{array}{rl}& \left(\frac{1+\mu y^{\ast }}{\alpha }\right)[\beta +b{y}^{\ast }\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}}][a-\frac{a}{k}\left(\frac{1+\mu y^{\ast }}{\alpha }\right)\left(\frac{d}{c\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}}}\right)\\ & \times \phantom{\left(\frac{d}{c\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}}}\right)}(\beta +b{y}^{\ast }\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}})-\frac{\alpha }{1+\mu y^{\ast }}]+\beta =0\end{array}$(2) ), and the inequality $p_1{q}_3>p_3{q}_1$ is true then coexisting equilibrium ${\overline{E}}_2=(x_1^{\ast },x_2^{\ast },y^{\ast })$ will be asymptotically stable.

6.2. Global stability analysis

Theorem 6.3

  

1. Predator free system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ) equilibrium $E=({\overline{x}}_1,{\overline{x}}_2)$ is globally asymptotically stable.

2. Co-existing equilibrium ${\overline{E}}_2=(x_1^{\ast },x_2^{\ast },y^{\ast })$ is globally asymptotically stable if $\begin{array}{rl}& -a{x}_1^{\ast }+\frac{a}{4k}(k+x_1^{\ast }{)}^2+\frac{\alpha x_1^{\ast }}{1+\mu i_3}-\frac{\beta i_2{x}_1^{\ast }}{s_1}-\frac{\alpha i_1{x}_2^{\ast }}{s_2(1+\mu s_3)}\\ & +\beta x_2^{\ast }+b{x}_2^{\ast }{s}_3\sqrt{w^2+{\gamma }^2}+\mathrm{d}{y}^{\ast }\le 0\end{array}$

Proof.

1. We construct a Lyapunov function as follows:

   $L=x_1-{\overline{x}}_1-{\overline{x}}_{1}ln(\frac{x_1}{{\overline{x}}_1})+x_2-{\overline{x}}_2-{\overline{x}}_{2}ln(\frac{x_2}{{\overline{x}}_2})$ (7) $\begin{array}{rl}\frac{\mathrm{d}L}{\mathrm{d}t}& =(1-\frac{{\overline{x}}_1}{x_1})\frac{\mathrm{d}{x}_1}{\mathrm{d}t}+(1-\frac{{\overline{x}}_2}{x_2})\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\\ & =\left(\frac{x_1-{\overline{x}}_1}{x_1}\right)[a{x}_1(1-\frac{x_1}{k})-\alpha x_1+\beta x_2]\\ & +\left(\frac{x_2-{\overline{x}}_2}{x_2}\right)(\alpha x_1-\beta x_2)\end{array}$(7) at equilibrium, $a{\overline{x}}_1(1-\frac{{\overline{x}}_1}{k})-\alpha {\overline{x}}_1+\beta {\overline{x}}_2=0$

   i.e. $\alpha =a(1-\frac{{\overline{x}}_1}{k})+\beta \frac{{\overline{x}}_2}{{\overline{x}}_1}$ and $\alpha {\overline{x}}_1=\beta {\overline{x}}_2\Rightarrow \beta =\frac{\alpha {\overline{x}}_1}{{\overline{x}}_2}$

   using above values of α and β into Equation (7(7) $\begin{array}{rl}\frac{\mathrm{d}L}{\mathrm{d}t}& =(1-\frac{{\overline{x}}_1}{x_1})\frac{\mathrm{d}{x}_1}{\mathrm{d}t}+(1-\frac{{\overline{x}}_2}{x_2})\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\\ & =\left(\frac{x_1-{\overline{x}}_1}{x_1}\right)[a{x}_1(1-\frac{x_1}{k})-\alpha x_1+\beta x_2]\\ & +\left(\frac{x_2-{\overline{x}}_2}{x_2}\right)(\alpha x_1-\beta x_2)\end{array}$(7) ) we get $\begin{array}{rl}\frac{\mathrm{d}L}{\mathrm{d}t}& =(x_1-{\overline{x}}_1)[a(1-\frac{x_1}{k})-a(1-\frac{{\overline{x}}_1}{k})-\beta \frac{{\overline{x}}_2}{{\overline{x}}_1}+\beta \frac{x_2}{x_1}]\\ & +(x_2-{\overline{x}}_2)[\alpha \frac{x_1}{x_2}-\alpha \frac{{\overline{x}}_1}{{\overline{x}}_2}]\end{array}$ or $\begin{array}{r}\frac{\mathrm{d}L}{\mathrm{d}t}=\frac{-a}{k}(x_1-{\overline{x}}_1{)}^2+\beta (x_1-{\overline{x}}_1)\left(\frac{{\overline{x}}_1{x}_2-{\overline{x}}_2{x}_1}{x_1{\overline{x}}_1}\right)+\alpha (x_2-{\overline{x}}_2)\left(\frac{{\overline{x}}_2{x}_1-{\overline{x}}_1{x}_2}{x_2{\overline{x}}_2}\right).\end{array}$ Since the carrying capacity k is the maximum population size of prey in the refuge, so, ${\overline{x}}_1<k$ and ${\overline{x}}_2<k$. Therefore, we have $\begin{array}{r}\frac{\mathrm{d}L}{\mathrm{d}t}\le \frac{-a}{k}(x_1-{\overline{x}}_1{)}^2-(x_1{\overline{x}}_2-x_2{\overline{x}}_1)[\frac{\beta }{k}\left(\frac{x_1-{\overline{x}}_1}{x_1}\right)-\frac{\beta }{k}\left(\frac{x_2-{\overline{x}}_2}{x_2}\right)]\end{array}$ or $\begin{array}{rl}\frac{\mathrm{d}L}{\mathrm{d}t}& \le \frac{-a}{k}(x_1-{\overline{x}}_1{)}^2-(x_1{\overline{x}}_2-x_2{\overline{x}}_1)\frac{\beta }{k}(-\frac{{\overline{x}}_1}{x_1}+\frac{{\overline{x}}_2}{x_2})\\ & \le \frac{-a}{k}(x_1-{\overline{x}}_1{)}^2-\frac{(x_1{\overline{x}}_2-x_2{\overline{x}}_1{)}^2}{x_1{x}_2}\frac{\beta }{k}\\ & \le 0.\end{array}$

2. From the second equation of (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ) and using the comparison theorem [31], we can obtain $\begin{array}{r}\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\le \alpha x_1-\beta x_2\end{array}$ At supremum of $x_2$, $\frac{\mathrm{d}{x}_2}{\mathrm{d}t}=0$, then (8) $\begin{array}{r}\alpha x_1\ge \beta x_2\end{array}$(8) Therefore, (9) $\begin{array}{r}{x}_2\le \frac{\alpha x_1}{\beta }.\end{array}$(9) From the first equation of system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ), we have $\begin{array}{r}\frac{\mathrm{d}{x}_1}{\mathrm{d}t}\le a{x}_1(1-\frac{x_1}{k})+\beta x_2\end{array}$ By using Equation (8(8) $\begin{array}{r}\alpha x_1\ge \beta x_2\end{array}$(8) ), we get $\begin{array}{r}\frac{\mathrm{d}{x}_1}{\mathrm{d}t}\le a{x}_1(1-\frac{x_1}{k})+\alpha x_1.\end{array}$ At supremum of $x_1$, $\frac{\mathrm{d}{x}_1}{\mathrm{d}t}=0$, then $\begin{array}{rl}& 0\le a-\frac{a{x}_1}{k}+\alpha \\ & x_1\le \frac{k}{a}(a+\alpha ).\end{array}$ Therefore, the supremum of $x_1=\frac{k}{a}(a+\alpha )=s_1$ and from (9(9) $\begin{array}{r}{x}_2\le \frac{\alpha x_1}{\beta }.\end{array}$(9) ) the supremum of $x_2=\frac{\alpha k}{\beta a}(a+\alpha )=s_2$

   Now using system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) )${}_3$, we have $\begin{array}{r}\frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}$ At supremum of y, $\frac{\mathrm{d}y}{\mathrm{d}t}=0$, therefore $\begin{array}{rl}& c{x}_2\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}=d\\ & c{s}_2\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\ge d\\ & \sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\ge \frac{d}{c{s}_2}\\ & e^{2\theta y}\le \frac{{\gamma }^2}{[{\left(\frac{d}{c{s}_2}\right)}^2-w^2]}\\ & y\le \frac{1}{2\theta }ln\left[\frac{{\gamma }^2}{[{\left(\frac{d}{c{s}_2}\right)}^2-w^2]}\right]\end{array}$ Therefore, the supremum of y = $\frac{1}{2\theta }ln[\frac{{\gamma }^2}{[(\frac{d}{c{s}_2}{)}^2-w^2]}]=s_3$.

   Similarly, from the first equation of system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ) we have $\begin{array}{r}\frac{\mathrm{d}{x}_1}{\mathrm{d}t}\ge a{x}_1(1-\frac{x_1}{k})-\alpha x_1\end{array}$ At infimum of $x_1$, $\frac{\mathrm{d}{x}_1}{\mathrm{d}t}=0$, therefore $\begin{array}{rl}& 0\ge a(1-\frac{x_1}{k})-\alpha \\ & x_1\ge \frac{k}{a}(a-\alpha )\end{array}$ Therefore, infimum of $x_1=\frac{k}{a}(a-\alpha )=i_1$.

   Now, from the third equation of system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) )${}_3$ we have $\begin{array}{r}\frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}$ At infimum of y, $\frac{\mathrm{d}y}{\mathrm{d}t}=0$, therefore $\begin{array}{rl}& c{x}_2\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}=d\\ & c{x}_2\sqrt{w^2+{\gamma }^2}\ge d\\ & x_2\ge \frac{d}{c\sqrt{w^2+{\gamma }^2}}.\end{array}$ Therefore, infimum of $x_2=\frac{d}{c\sqrt{w^2+{\gamma }^2}}=i_2$.

   Now, from the second equation of system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ) we have $\begin{array}{rl}& \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}\ge \frac{\alpha i_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\end{array}$ At infimum of $x_2$, $\frac{\mathrm{d}{x}_2}{\mathrm{d}t}=0$, therefore $\begin{array}{rl}& \beta x_2+b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\ge \frac{\alpha i_1}{1+\mu y}\\ & \beta s_2+b{s}_{2}y\sqrt{w^2+{\gamma }^2}\ge \frac{\alpha i_1}{1+\mu y}\\ & y^2[b{s}_2\mu \sqrt{w^2+{\gamma }^2}]+y[[\mu \beta s_2+b{s}_2\sqrt{w^2+{\gamma }^2}]]+(\beta s_2-\alpha i_1)\ge 0.\end{array}$ Therefore, infimum of $y=i_3$ where

   $i_3=\frac{-[\mu \beta s_2+b{s}_2\sqrt{w^2+{\gamma }^2}]+\sqrt{(\mu \beta s_2+b{s}_2\sqrt{w^2+{\gamma }^2}{)}^2-4(b{s}_2\mu \sqrt{w^2+{\gamma }^2})(\beta s_2-\alpha i_1)}}{2b{s}_2\mu \sqrt{w^2+{\gamma }^2}}$

   provided $\beta <\frac{\alpha i_1}{s_2}$.

   Now, we construct a Lyapunov function as follows: $\begin{array}{rl}L& =x_1-x_1^{\ast }-x_1^{\ast }ln\left(\frac{x_1}{x_1^{\ast }}\right)+x_2-x_2^{\ast }-x_2^{\ast }ln\left(\frac{x_2}{x_2^{\ast }}\right)+y-y^{\ast }-y^{\ast }ln\left(\frac{y}{y^{\ast }}\right)\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =(1-\frac{x_1^{\ast }}{x_1})\frac{\mathrm{d}{x}_1}{\mathrm{d}t}+(1-\frac{x_2^{\ast }}{x_2})\frac{\mathrm{d}{x}_2}{\mathrm{d}t}+(1-\frac{y^{\ast }}{y})\frac{\mathrm{d}y}{\mathrm{d}t}\\ & =\left(\frac{x_1-x_1^{\ast }}{x_1}\right)[a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2]\\ & +\left(\frac{x_2-x_2^{\ast }}{x_2}\right)[\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}]\\ & +\left(\frac{y-y^{\ast }}{y}\right)[c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y]\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =(x_1-x_1^{\ast })[a(1-\frac{x_1}{k})-\frac{\alpha }{1+\mu y}+\frac{\beta x_2}{x_1}]\\ & +(x_2-x_2^{\ast })[\frac{\alpha x_1}{x_2(1+\mu y)}-\beta -by\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}]\\ & +(y-y^{\ast })[c{x}_2\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-d]\\ \frac{\mathrm{d}L}{\mathrm{d}t}& =-\frac{a}{k}{(x_1-\frac{k+x_1^{\ast }}{2})}^2-a{x}_1^{\ast }+\frac{a}{4k}(k+x_1^{\ast }{)}^2+\frac{\alpha x_1}{1+\mu y}-\frac{\beta x_2{x}_1^{\ast }}{x_1}\\ & -\frac{\alpha x_1{x}_2^{\ast }}{x_2(1+\mu y)}+\beta x_2^{\ast }-(b-c)x_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}+b{x}_2^{\ast }y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & -c{x}_2{y}^{\ast }\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y+\mathrm{d}{y}^{\ast }\\ \frac{\mathrm{d}L}{\mathrm{d}t}& \le -a{x}_1^{\ast }+\frac{a}{4k}(k+x_1^{\ast }{)}^2+\frac{\alpha x_1^{\ast }}{1+\mu y}-\frac{\beta x_2{x}_1^{\ast }}{x_1}-\frac{\alpha x_1{x}_2^{\ast }}{x_2(1+\mu y)}\\ & +\beta x_2^{\ast }+b{x}_2^{\ast }y\sqrt{w^2+{\gamma }^2}+\mathrm{d}{y}^{\ast }\\ \frac{\mathrm{d}L}{\mathrm{d}t}& \le -a{x}_1^{\ast }+\frac{a}{4k}(k+x_1^{\ast }{)}^2+\frac{\alpha x_1^{\ast }}{1+\mu i_3}-\frac{\beta i_2{x}_1^{\ast }}{s_1}-\frac{\alpha i_1{x}_2^{\ast }}{s_2(1+\mu s_3)}\\ & +\beta x_2^{\ast }+b{x}_2^{\ast }{s}_3\sqrt{w^2+{\gamma }^2}+\mathrm{d}{y}^{\ast }.\end{array}$ Therefore, $\begin{array}{rl}\frac{\mathrm{d}L}{\mathrm{d}t}& \le 0\mathbf{i}\mathbf{f}-a{x}_1^{\ast }+\frac{a}{4k}(k+x_1^{\ast }{)}^2+\frac{\alpha x_1^{\ast }}{1+\mu i_3}-\frac{\beta i_2{x}_1^{\ast }}{s_1}-\frac{\alpha i_1{x}_2^{\ast }}{s_2(1+\mu s_3)}+\beta x_2^{\ast }\\ & +b{x}_2^{\ast }{s}_3\sqrt{w^2+{\gamma }^2}+\mathrm{d}{y}^{\ast }\le 0.\end{array}$

7. Transcritical bifurcation

We investigate the transcritical bifurcation for the system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ) taking $d=d_0$ as the bifurcation parameter. According to the Sotomayor theorem [22] for local bifurcation, the following conditions are satisfied:

1. $w^T{f}_d({\overline{E}}_1,d_0)=0$

2. $w^T(D{f}_d({\overline{E}}_1,d_0).\upsilon )\ne 0$

3. $w^T{D}^{2}f({\overline{E}}_1,d_0)(\upsilon ,\upsilon )\ne 0.$

Then the system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ) experiences a transcritical bifurcation at the free predator equilibrium point ${\overline{E}}_1=(k,\frac{\alpha k}{\beta },0)$ as the parameter d passes through the bifurcation value $d=d_0=c{\overline{x}}_2\sqrt{w^2+{\gamma }^2}$ where Df is denoted by the matrix of partial derivatives of the components of f with respect to the components of $X=(x_1,x_2,y)$ and $f_d$ is denoted for the vector of partial derivatives of the components of f with respect to the scalar d. $(\upsilon )$ is an eigenvector of $A=Df({\overline{E}}_1,d_0)$ corresponding to the eigenvalue ${\lambda }_3=0$ and ω is an eigenvector of $A^T$ corresponding to the eigenvalue ${\lambda }_3=0$,

$Df(X)=\left(\begin{array}{ccc}\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{x}_1}& \frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{x}_2}& \frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }y}\\ \frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{x}_1}& \frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{x}_2}& \frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }y}\\ \frac{\mathrm{\partial }{f}_3}{\mathrm{\partial }{x}_1}& \frac{\mathrm{\partial }{f}_3}{\mathrm{\partial }{x}_2}& \frac{\mathrm{\partial }{f}_3}{\mathrm{\partial }y}\end{array}\right)$

$D^{2}f(X)=\left(\begin{array}{ccccccccc}\frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }{x}_1\mathrm{\partial }{x}_1}& \frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }{x}_1\mathrm{\partial }{x}_2}& \frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }{x}_1\mathrm{\partial }y}& \frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }{x}_2\mathrm{\partial }{x}_1}& \frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }{x}_2\mathrm{\partial }{x}_2}& \frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }{x}_2\mathrm{\partial }y}& \frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }y\mathrm{\partial }{x}_1}& \frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }y\mathrm{\partial }{x}_2}& \frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }y\mathrm{\partial }y}\\ \frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_1\mathrm{\partial }{x}_1}& \frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_1\mathrm{\partial }{x}_2}& \frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_1\mathrm{\partial }y}& \frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_2\mathrm{\partial }{x}_1}& \frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_2\mathrm{\partial }{x}_2}& \frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_2\mathrm{\partial }y}& \frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }y\mathrm{\partial }{x}_1}& \frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }y\mathrm{\partial }{x}_2}& \frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }y\mathrm{\partial }y}\\ \frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_1\mathrm{\partial }{x}_1}& \frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_1\mathrm{\partial }{x}_2}& \frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_1\mathrm{\partial }y}& \frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_2\mathrm{\partial }{x}_1}& \frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_2\mathrm{\partial }{x}_2}& \frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_2\mathrm{\partial }y}& \frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }y\mathrm{\partial }{x}_1}& \frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }y\mathrm{\partial }{x}_2}& \frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }y\mathrm{\partial }y}\end{array}\right)$

Theorem 7.1

The system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ) undergoes a transcritical bifurcation at the free predator equilibrium point ${\overline{E}}_1=(k,\frac{\alpha k}{\beta },0)$ when $d=d_0=c{\overline{x}}_2\sqrt{w^2+{\gamma }^2}$.

The linearized system around the equilibrium point ${\overline{E}}_1$ is $\begin{array}{r}\left(\begin{array}{ccc}-(a+\alpha )& \beta & \alpha \mu k\\ \alpha & -\beta & -\alpha {\overline{x}}_1\mu -b{\overline{x}}_2\sqrt{w^2+{\gamma }^2}\\ 0& 0& c{\overline{x}}_2\sqrt{w^2+{\gamma }^2}-d\end{array}\right)\end{array}$

Now ${\lambda }_3=0$ at $d=c{\overline{x}}_2\sqrt{w^2+{\gamma }^2}$, So $d_0=c{\overline{x}}_2\sqrt{w^2+{\gamma }^2}$ $\begin{array}{r}J({\overline{E}}_1,d_0)=\left(\begin{array}{ccc}-(a+\alpha )& \beta & \alpha \mu k\\ \alpha & -\beta & -\alpha {\overline{x}}_1\mu -b{\overline{x}}_2\sqrt{w^2+{\gamma }^2}\\ 0& 0& 0\end{array}\right)\end{array}$ Let us define $\upsilon =({\upsilon }_1,{\upsilon }_2,{\upsilon }_3{)}^T$ and $\omega =({\omega }_1,{\omega }_2,{\omega }_3{)}^T$ are respectively the right and left eigenvectors of ${\lambda }_3=0$.

Now solving $J({\overline{E}}_1,d_0)\upsilon =0$ $\begin{array}{r}\left(\begin{array}{c}-(a+\alpha ){\upsilon }_1+\beta {\upsilon }_2+\mu \alpha k{\upsilon }_3\\ \alpha {\upsilon }_1-\beta {\upsilon }_2-(\mu \alpha k+b{\overline{x}}_2\sqrt{w^2+{\gamma }^2}){\upsilon }_3\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right)\end{array}$ implies that

$\upsilon =((\frac{-b}{a}{\overline{x}}_2\sqrt{w^2+{\gamma }^2}){\upsilon }_3,(\frac{-\alpha b}{a\beta }{\overline{x}}_2\sqrt{w^2+{\gamma }^2}-(\frac{\mu \alpha k}{\beta }+\frac{b}{\beta }{\overline{x}}_2\sqrt{w^2+{\gamma }^2})){\upsilon }_3,{\upsilon }_3)$, where ${\upsilon }_3$ is any nonzero real number.

Then solving $J^T({\overline{E}}_1,d_0)\omega =0$ $\begin{array}{r}\left(\begin{array}{c}-(a+\alpha ){\omega }_1+\alpha {\omega }_2\\ \beta {\omega }_1-\beta {\omega }_2\\ \mu \alpha k{\omega }_1-(\mu \alpha k+b{\overline{x}}_2\sqrt{w^2+{\gamma }^2}){\omega }_2\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right)\end{array}$ we get ${\omega }_1={\omega }_2=0$ and ${\omega }_3$ is any nonzero real number, so $\omega =(0,0,{\omega }_3)$.

Now, system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ) can be rewritten as in the following vector form: (10) $\begin{array}{r}\dot{X}=f(X),\end{array}$(10) where $X=(x_1,x_2,y)$ and $\begin{array}{r}f(X)=\left(\begin{array}{c}a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ \frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-dy\end{array}\right)\end{array}$

Taking the derivative of $f(X)$ with respect to d, we get (11) $\begin{array}{r}{f}_d(X)=\left(\begin{array}{c}0\\ 0\\ -y\end{array}\right)\end{array}$(11) then $f_d({\overline{E}}_1,d_0)=\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right)$, Hence ${\omega }^T{f}_d({\overline{E}}_1,d_0)=0$, So, the first condition is satisfied.

Next, taking the derivative of $f_d(X)$ with respect to $X=(x_1,x_2,y{)}^T$, we get $\begin{array}{r}D{f}_d(X)=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& -1\end{array}\right)\end{array}$ then $\begin{array}{r}D{f}_d({\overline{E}}_1,d_0)=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& -1\end{array}\right)\end{array}$ Therefore, we have ${\omega }^{T}D{f}_d({\overline{E}}_1,d_0).\upsilon =$ $\left(\begin{array}{ccc}0& 0& {\omega }_3\end{array}\right)\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& -1\end{array}\right)$ $\left(\begin{array}{c}{\upsilon }_1\\ {\upsilon }_2\\ {\upsilon }_3\end{array}\right)=-{\omega }_3{\upsilon }_3\ne 0$.

Therefore, the second condition is satisfied.

Furthermore, $\begin{array}{r}{D}^{2}f(X)=\left(\begin{array}{ccccccccc}{A}_1& A_2& A_3& A_4& A_5& A_6& A_7& A_8& A_9\\ A_{10}& A_{11}& A_{12}& A_{13}& A_{14}& A_{15}& A_{16}& A_{17}& A_{18}\\ A_{19}& A_{20}& A_{21}& A_{22}& A_{23}& A_{24}& A_{25}& A_{26}& A_{27}\end{array}\right)\end{array}$ where

$A_1=\frac{-2a}{k}$, $A_2=0$, $A_3=\frac{\alpha \mu }{(1+\mu y{)}^2}$, $A_4=0$, $A_5=0$, $A_6=0$, $A_7=\frac{\alpha \mu }{(1+\mu y{)}^2}$, $A_8=0$, $A_9=\frac{-2\alpha x_1{\mu }^2}{(1+\mu y{)}^3}$, $A_{10}=0$, $A_{11}=\frac{-\alpha \mu }{(1+\mu y{)}^2}$, $A_{12}=0$, $A_{13}=0$, $A_{14}=0$, $A_{15}=-b\sqrt{w^2+{\gamma }^2}$, $A_{16}=\frac{-\alpha \mu }{(1+\mu y{)}^2}$, $A_{17}=-b\sqrt{w^2+{\gamma }^2}$, $A_{18}=\frac{2\alpha x_1{\mu }^2}{(1+\mu y{)}^3}$, $A_{19}=0$, $A_{20}=0$, $A_{21}=0$, $A_{22}=0$, $A_{23}=0$, $A_{24}=c\sqrt{w^2+{\gamma }^2}$,$A_{25}=0$, $A_{26}=c\sqrt{w^2+{\gamma }^2}$, $A_{27}=0.$

The tensor product $\begin{array}{r}(\upsilon ,\upsilon )=\left(\begin{array}{c}{\upsilon }_1^2\\ {\upsilon }_1{\upsilon }_2\\ {\upsilon }_1{\upsilon }_3\\ {\upsilon }_2{\upsilon }_1\\ {\upsilon }_2^2\\ {\upsilon }_2{\upsilon }_3\\ {\upsilon }_3{\upsilon }_1\\ {\upsilon }_3{\upsilon }_2\\ {\upsilon }_3^2\end{array}\right)\end{array}$ so, $\begin{array}{rl}& {\omega }^T{D}^{2}f(E_1,d_0)(\upsilon ,\upsilon )\\ & =\left(\begin{array}{ccc}0& 0& {\omega }_3\end{array}\right)\left(\begin{array}{ccccccccc}{A}_1& A_2& A_3& A_4& A_5& A_6& A_7& A_8& A_9\\ A_{10}& A_{11}& A_{12}& A_{13}& A_{14}& A_{15}& A_{16}& A_{17}& A_{18}\\ A_{19}& A_{20}& A_{21}& A_{22}& A_{23}& A_{24}& A_{25}& A_{26}& A_{27}\end{array}\right)\left(\begin{array}{c}{\upsilon }_1^2\\ {\upsilon }_1{\upsilon }_2\\ {\upsilon }_1{\upsilon }_3\\ {\upsilon }_2{\upsilon }_1\\ {\upsilon }_2^2\\ {\upsilon }_2{\upsilon }_3\\ {\upsilon }_3{\upsilon }_1\\ {\upsilon }_3{\upsilon }_2\\ {\upsilon }_3^2\end{array}\right)\\ & =2c\sqrt{w^2+{\gamma }^2}{\omega }_3({\upsilon }_2{\upsilon }_3)\ne 0.\end{array}$ Therefore, according to the Sotomayor theorem [22] for local bifurcation, system (1(1) $\begin{array}{r}\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=a{x}_1(1-\frac{x_1}{k})-\frac{\alpha x_1}{1+\mu y}+\beta x_2\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\alpha x_1}{1+\mu y}-\beta x_2-b{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c{x}_{2}y\sqrt{w^2+{\gamma }^2{e}^{-2\theta y}}-\mathrm{d}y\end{array}\end{array}$(1) ) has a transcritical bifurcation at steady state at the free predator equilibrium point ${\overline{E}}_1=(k,\frac{\alpha k}{\beta },0)$ as the parameter d passes through the bifurcation value $d=d_0=c{\overline{x}}_2\sqrt{w^2+{\gamma }^2}$ (Figure 1).

Figure 1. The transcritical bifurcation at d = 0.266204905858465, the parameters are $a=0.7,k=0.8,\alpha =0.035,\mu =5,\beta =0.0119,b=0.0112,c=0.04,d=0.07,\theta =0.69,\gamma =2,w=2$. (The red colour signifies the stable equilibrium and the black color shows the unstable equilibrium).

8. Numerical analysis

In this section, numerical simulation results are given to confirm the above theoretical results. Table 1 contains the parameter values that are used, some of them might be changed to study their effect.

Table 1. Parameter values used for simulations.

Parameter	Value	Units
a	0.7	Time${}^{-1}$
k	0.8	$x_1$ density × length${}^{-2}$
α	0.035	Time${}^{-1}$
μ	5	Time${}^{-1}$
β	0.0119	Time${}^{-1}$
b	0.0112	Time${}^{-1}$
w	2	Length × time${}^{-1}$
γ	2	Length × time ${}^{-1}$
θ	0.69	y density × length${}^{-2}$
c	0.04	Time${}^{-1}$
d	0.07	Time${}^{-1}$

Figures 2 and 3 show the convergence of the free predator equilibrium point $(k,\frac{\alpha k}{\beta },0)$ and the coexistence equilibrium point $(x_1^{\ast },x_2^{\ast },y^{\ast })$ respectively.

Figure 2. The convergence of the free predator equilibrium point $(k,\frac{\alpha k}{\beta },0)$, the parameters are $a=0.7,k=0.8,\alpha =0.035,\mu =5,\beta =0.0119,b=0.0112,c=0.04,d=0.266205,\theta =0.69,\gamma =2,\omega =2$, ${\overline{E}}_1=(0.8,2.352941176,0)$.

Figure 3. The convergence of the coexistence equilibrium point ${\overline{E}}_2$, the parameters are a = 0.7, k = 0.8, $\alpha =0.035$, $\mu =5,\beta =0.0119,b=0.0112,c=0.04,d=0.07,\theta =0.69,\gamma =2,\omega =2$.

8.1. Impact of predator and prey velocities on their populations

From figures 4 we find that increasing the predator velocity has no effect on the prey population in the refuge but this will decrease the prey population outside the refuge habitat. On the other side, when the predator velocity rises, so does the predator population.

Figure 4. The effect of predator velocity on the population densities, parameters are $a=0.7,k=0.8,\alpha =0.035,\mu =5,b=0.0112,\beta =0.0119,c=0.04,d=0.07,\theta =0.69,\gamma =2$.

Figure 5 represents that as the prey velocity increases, the predator population grows while the prey population outside the refuge decreases. Furthermore, the prey population in the refuge habitat is unaffected by the increase in the prey velocity. Hence the high speed of prey becomes more vulnerable to predators. Therefore, any increase in predator or prey velocity raises the prey mortality and hence the prey population declines outside the refuge. Refuge area is safe from predation so, the prey population in the refuge reaches its equilibrium value. Numerical results confirm the claim that the slow speed of prey is useful for the survival of the prey population and at the same time the high speed of the predator is good for the predator population.

Figure 5. The effect of prey velocity on the population densities, parameters are $a=0.7,k=0.8,\alpha =0.035,\mu =5,b=0.0112,\beta =0.0119,c=0.04,d=0.07,\theta =0.69,w=2$.

8.2. Effect of antipredator sensitivity (θ)

As demonstrated in figures 6, raising antipredator sensitivity increases prey density outside of the refuge habitat and as a result of this decreases the predator population density. The prey density inside the refuge is uneffected by antipredator sensitivity and will remain constant.

Figure 6. The effect of antipredator sensitivity on the population densities, parameters are $a=0.7,k=0.8,\alpha =0.035,\mu =5,b=0.0112,\beta =0.0119,c=0.04,d=0.07,w=2,\gamma =2$.

8.3. Effect of refuge parameter $(\beta )$

Figure 7 indicates that as the refuge parameter is increased, then more and more prey will be able to take shelter inside the refuge habitat and so prey population inside the refuge increases until it meet carrying capacity. Due to this process the prey population outside the refuge declines and so the predator population which depends on the prey population will also decline and eventually go extinct.

Figure 7. The effect of refuge parameter on the population densities, parameters are $a=0.7,k=0.8,\alpha =0.035,\mu =5,b=0.0112,c=0.04,d=0.07,\theta =0.69,\gamma =2,w=2$.

8.4. Effect of fear parameter $(\mu )$

Figure 8 represents that as the fear of the predator increases the prey density inside the refuge habitat increases because of predator fear the prey will prefer to stay safe inside the refuge. The population of the prey outside the refuge and the predator population decreases. Since predators will get less food due to less number of prey available outside the refuge. So, by heavy predation on less number of prey the population of prey declines and ultimately the predator population will also decline.

Figure 8. The effect of fear parameter on the population densities, parameters are $a=0.7,k=0.8,\alpha =0.035,\beta =0.0119,b=0.0112,c=0.04,d=0.07,\theta =0.69,\gamma =2,w=2$.

9. Conclusion

Predators provide a threat to all creatures, and they are all at risk of being devoured. Fear of predators causes prey populations to change their behaviour. We investigated the dynamics of predator–prey interactions in two settings, namely, refuge and out-of refuge, in this article. The prey is protected from predatory slaughter in the refuge, and it has enough supplies to survive, and the sanctuary's population rises logistically. Predators engage with prey outside of their safe haven and may kill them. Prey lives in fear of predators, but as that fear fades, prey emerges from their hiding places. Prey seeks refuge in the refuge as their fear of predation grows. By allowing predator density to influence prey velocity, we included predator and prey velocities with antipredator behaviour in our model.

Three biologically viable equilibria are obtained, and their stability is discussed. The equilibrium without prey and predator populations will always exist, but it will be unstable, implying that prey and predator populations will never go extinct. There will always be an equilibrium with zero predator population and a nonzero prey population, and it is globally stable. The nonzero equilibrium will be asymptotically stable if $(\frac{a{x}_1^{\ast }}{k}+\frac{\beta x_2^{\ast }}{x_1^{\ast }})(\frac{\alpha x_1^{\ast }\mu }{(1+\mu y^{\ast }{)}^2}+\frac{b{x}_2^{\ast }}{\sqrt{w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}}}[w^2+{\gamma }^2{e}^{-2\theta y^{\ast }}(1-\theta y^{\ast })])>\frac{{\alpha }^2{x}_1^{\ast }\mu }{(1+\mu y^{\ast }{)}^3}$. At a specific value of the predator's mortality rate, a transcritical bifurcation occurs at predator free equilibrium.

We conducted a numerical analysis to see how the factors affect our system, and we observed that prey outside the refuge habitat and predator population decline and become extinct as the refuge parameter rises, whereas prey in the refuge increases until it meets the carrying capacity. Moreover, the prey density inside the refuge habitat grows as the fear parameter increases, but the prey density outside the refuge and predator density decreases. Furthermore, we discovered that when the antipredator sensitivity is low, predator and prey velocities have the same influence on the final equilibrium prey density. Any increase in predator or prey velocity causes prey mortality to rise and the equilibrium prey density to decline. However, boosting antipredator sensitivity increases the density of prey outside the refuge habitat while reducing the density of the predator population.

Disclosure statement

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

Additional information

Funding

We are very thankful to the reviewers for their valuable comments and suggestions which led to an improvement of our original manuscript. Research is funded by Sultan Qaboos University, Oman.