Ecological system with fear induced group defence and prey refuge

Shivam Yadava, Jai Prakash Tripathia 111Corresponding author: Jai Prakash
Tripathi, Department of Mathematics, Indian Institute of Technology Central
University of Rajasthan, Bandar Sindri, Kishangarh, Ajmer 305817, Rajasthan,
India, Email<EMAIL_ADDRESS>Shrichand Bhuriaa, Satish Kumar
Tiwari${}^{b},$ Deepak Tripathia, Vandana Tiwari${}^{c},$ Ranjit Kumar
Upadhyayd, Yun Kange
a Department of Mathematics, Central University of Rajasthan, Bandar Sindri,
Kishangarh
Ajmer 305817, Rajasthan, India
b Assam Energy Institute (Centre of Rajiv Gandhi Institute of Petroleum
Technology, Jais, Amethi),
Sivasagar, Assam - 785697, India
c Gandhi Smarak PG College, Samodhpur, Jaunpur 223102, U.P., India
d Department of Mathematics & Computing, Indian Institute of Technology
(Indian School of Mines),
Dhanbad 826004, Jharkhand, India
$e$ College of Integrative Sciences and Arts, Arizona State University, Mesa,
AZ85212, USA.

###### Abstract

In this study, we investigate the dynamics of a spatial and non spatial prey-
predator interaction model that includes the following: (i) fear effect
incorporated in prey birth rate; (ii) group defence of prey against predators;
and (iii) prey refuge. We provide comprehensive mathematical analysis of
extinction and persistence scenarios for both prey and predator species. We
investigate how the prey and predator equilibrium densities are influenced by
the prey birth rate and fear level. To better explore the dynamics of the
system, a thorough investigation of bifurcation analysis has been performed
using fear level, prey birth rate, and prey’s death rate caused by intra-prey
competition as bifurcation parameter. All potential occurrences of bi-
stability dynamics have also been investigated for some relevant sets of
parametric values. Our numerical evaluations show that high levels of fear can
stabilize the prey-predator system by ruling out the possibility of periodic
solutions. Also, our model’s Hopf bifurcation is subcritical in contrast to
traditional prey-predator models, which ignore the cost of fear and have
supercritical Hopf bifurcations in general. In contrast to the general trend,
predator species go extinct at higher values of prey birth rates. We have also
found that, contrary to the typical tendency for prey species to go extinct,
both prey and predator populations may coexist in the system as intra-prey
competition level grows noticeably. We have also been obtained that both prey
and predator equilibrium densities increase (decrease) as the prey birth rate
(fear level in prey) increases. The stability and Turing instability of
associated spatial model have also been investigated analytically. We also
perform the numerical simulation to observe the effect of different parameters
on the density distribution of species. Different types of spatiotemporal
patterns like spot, mixture of spots and stripes have been observed via
variation of time evolution, diffusion coefficient of predator population,
level of fear factor and prey refuge. The fear level parameter (k) has a great
impact on the spatial dynamics of model system.

Keywords: Prey-predator interaction; Fear effect; Prey refuge; Bi-stability;
Turing instability, Pattern formation

##  1 Introduction

Traditionally, predators are thought to directly kill prey, reducing prey
survival and numbers, and that this is the end of their ecological role in
prey-predator interactions. However, numerous elegant experiments [1, 2, 3, 4]
have demonstrated that the role of predators in ecology is more profound than
just the direct killing of prey, as it may cause behavioral, psychological,
and neurobiological changes in the prey population, ultimately affecting prey
demography. In reality, when there is a confrontation between prey and
predator, the prey perceives a predation risk. This perceived predation risk
leads to fear or stress in the prey. Until relatively recently, it has been
considered that predator-induced stress (fear) is necessarily acute,
transitory, and does not affect the demographical process much. However, many
experiments were conducted at the beginning of the $1990$s to determine
whether predator-induced fear in response to predator could truly affect prey
population even when direct killing was purposefully avoided [5, 6]. In such
studies, direct killing of prey by predators has often been purposefully
prevented by completely shutting (e.g. [7, 8]) or partially dissecting (e.g.
[1]) the predator population’s mouths (e.g. spiders, stoneflies, and damsel
bugs). These risk predators [8] are then forcibly enclosed in artificial
enclosures with their prey (such as grasshopper nymphs, mayfly larvae, and pea
aphids), as they can not directly kill but only ‘intimidate’ any prey [6]. The
influence of intimidation alone on prey demography is then measured by
comparing prey kept in enclosures with ‘risk’ predators to prey housed by
themselves.

In some other investigations, prey have been ‘intimidated’ using caged
predators or the smell of predators [9]. In a meta analysis of such studies,
Preisser, Bolnick, and Benard [6] found that, on average, only predator
intimidation has an equal or greater impact on prey birth rate, survival, and
growth than direct death. Travers et al.[2] conducted an experiment to show
that female song sparrows who frequently had their nests preyed upon produced
fewer eggs in subsequent clutches and, as a result, had a lower birth rate
than song sparrows who didn’t. Eggers et al.[3] also showed that Siberian jays
produced fewer eggs in their first clutch of the season as a result of
predator call playbacks only. These studies demonstrate that prey reproductive
success and long-term survival may be negatively impacted by exposure to
predators or predator cues. However, none of these studies quantified the
impact of risk predators (fear) on prey demography. However, in 2011, Zanette
et al. [4] did revolutionary work by performing a prolific field experiment in
which direct killing was actively prevented to quantify the impact of risk
predators or predator’s fear on the prey population. In this study, it was
claimed that ‘intimidation’ (perceived predation risk) alone, which was
controlled by predator call playbacks, caused a $40\%$ reduction in the number
of song sparrow chicks produced annually. Female song sparrows who were afraid
of predators produced fewer eggs, more of which failed to hatch due to
interruptions in incubation, and more nestlings who starved to death as a
result of intimidation hampering their foraging activities. These experiments
show that as a consequence of the fear of predators (predator induced stress),
fecundity and survival may be reduced, and the total reduction in the prey
number because of exposure to predators may thus far exceed that due to direct
killing alone. The most convincing reason to believe that fear will affect
prey populations universally is that fearful prey eat less because they can’t
have their heads down to forage and their heads up for vigilance of predators
at the same time[10, 11, 4, 12].

The fear may result in the process of avoiding predators for which prey shows
some anti-predator responses. A simple but effective anti-predator response
involves the use of refuges. Refuge is a concept in which prey protects
themselves from predation by hiding in an area or being in a situation where
it is less accessible or difficult to find. Spatial refuges (burrows, dense
foliage), temporal refuges (temporarily migrating to a low-risk habitat),
forming groups, or simply reducing foraging activity are examples of refuges
that reduce the confrontation of prey with predators[13, 14]. Generally, prey
uses spatial refuges as these are easily available because of environmental
heterogeneity. By using refuges, the prey population is not entirely exposed
to the predators. A certain proportion of the prey population is always
shielded from predators. The benefit of refuge is that it reduces the
predation risk. However, the longer the prey stays in the refuge, the fewer
feeding and mating opportunities it has. For example, if the prey uses the
refuge in the form of a group formation, this can lead to enhanced competition
among the prey population [15]; where the refuge is due to reduced foraging
activity, this probably also cuts down on the prey’s encounter with feeding
opportunity; and where the prey uses the refuge in the form of space or time,
it can cause a shift away from high-energy feeding habitats or activity times.
As a result, spending more time in refuges not only promotes survival, but it
also reduces prey feeding rates. There are two consequences of the cost of
refuge use in lowering prey feeding rate: it reduces the possibility that prey
will escape predator control and works against any stabilising effects of
refuge utilisation. The presence of refuges has a substantial impact on the
coexistence of predator and prey in the system[16]. The only two types of prey
refuge have been considered in extant models; (i) a constant number of the
prey population in refuges[17, 18, 19, 20] (ii) a constant proportion [21, 14,
22, 17, 19] of the prey population in refuges. The investigation of these
refuges has persisted in an elementary stage because incorporating these
‘pseudo-refuges’ into the prey-predator models has proven fairly challenging.
It has been shown that use of refuge enhances the coexistence of prey and
predator systems [23].

Other than using refuge, a number of anti-predator responses, including toxin
release, forming groups, playing dead, mimicry, herd behaviour, camouflage,
and apostatic selection are also triggered in prey by fear of encroaching
predators[24, 25, 26]. Prey grouping is well established as an anti-predator
behaviour in many prey species[27]. Additionally, some researchers have looked
into how group size changes in response to predator risk [28, 29]. Predation
is decreased or even avoided as a result of group defence by prey against
predators due to the greater inclination of the prey to better defend or hide
themselves when their density is high enough. Predators can sometimes find
larger groups easily, though. It is commonly acknowledged that enhancing a
prey’s protection against a predator is costly and causes a decrease in the
amount of food intake of the prey[30]. When the prey is aware of the threat,
it displays a particular sensitivity to the predator for physiological
reasons. As their sensitivity to predators develops, prey often decreases the
amount of time they spend outside and strengthens their group defence, which
decreases their chances of becoming prey. In other words, a prey individual
has a higher probability of surviving due to their greater sensitivity to the
predator. The prey will reduce their reproduction rate as a result of putting
more time and effort into group defence. As a result, the prey’s sensitivity
to the predator may promote the protection of the individual prey while also
hastening the reduction in predator population density. Therefore, many prey
(such as ungulates, ungulata) form huge groups in order to increase their
chances of survival through ‘safety in number’ and lower the cost of
individual defence by using ‘many eyes’ effects[31]. The selection process
will favour those who strike the optimal balance between risk reduction’s
benefits and costs.

The question now is how to incorporate the increased instant survival of prey
and the cost due to group defence simultaneously through functional response.
All the initial functional responses Holling I, II, and III [32], are
monotonically increasing, which means that more the prey is in the
environment, the better for the predators[33]. However, this cannot be true
for functional responses that represent group defense. Numerous studies and
observations also show that monotonicity is not always true and that when the
number of prey is high, a non-monotonic response could emerge[34, 35].
Therefore, taking care above critical issues, Andrews[36] proposed a Holling
type IV functional response $\left(i.e.~{}~{}p(x)=\frac{\beta
x}{a+bx+x^{2}}\right)$, which is the most frequent and straightforward
technique to model the group defence behaviour of prey. When prey density is
high, the rate of predation per predator $p(x)$ decreases until it reaches
zero. Therefore, when functional response of Holling type IV is taken as group
defence function, predators may not survive above a certain prey density
threshold because of the increased prey group defence. Another facet of group
defence is that when prey populations congregate into herds, the weakest
individuals are allowed to inhabit the centre of the herd, leaving healthier
and stronger individuals on the herd’s periphery. In any event, it’s crucial
that the more individualistic population only interacts with the more
socialised one at the herd’s perimeter. When attacks occur, especially in
prey-predator or competitive models, it is typically the individuals at the
herd’s edge who suffer the consequences of the predator’s deeds, which is
reflected by the square root term[37]. In present study, we take into account
the group defence of prey in the Holling type IV functional response. Freedman
and Wolkowicz [38] were the first to present a mathematical model of group
defence.

In $2016$,Wang et al. [39] analysed the prey-predator interaction,
incorporating the cost of fear for the first time, and opening a wide scope in
population ecology. He demonstrated that a high level of fear can have a
stabilizing impact on the prey-predator system with Holling type II by
eliminating the periodic solution. The author of [40] showed that mutual
interference can stabilise the system by altering the extinction of predators
after incorporating it into the Wang et al. model[39]. According to study by
Kar et al.[41] on a prey-predator model that incorporates a Holling-type-II
functional response and prey refuge, an increase in refuge may cause
population outbreaks. In [42], the author extended the model of Kar et al.[41]
by incorporating the fear factor into the model and found that fear can
enhance system’s stability by eliminating the occurrence of limit cycle
oscillation or periodic solution. In [43], the author analysed the prey-
predator model, incorporating the fear effect and the prey refuge along with
the Holling type III functional response. It was discovered that the fear
effect can stabilise the system by removing periodic solutions and reducing
the final number of predator species at the coexist equilibrium, but it cannot
cause predator extinction. Yuxin dong et al.[44] demonstrated that even when
there is a high level of fear, the prey population remains stable rather than
extinction. The high level of predator taxis sensitivity thwarts predation.
Sasmal and Takeuchi[30] studied the Hopf bifurcation and multi-stability in a
prey-predator system, incorporating group defence and the fear effect.
According to Sasmal, the high predator-taxis sensitivity may cause instability
in the system and the extinction of predator populations.

On the basis of the assumption that predators and preys move randomly in the
habitat is modeled as diffusive prey-predator model. In reality, the spatial
movement describes that how predator pursuing prey and in what way prey
escapes from predator [45]. Biodiversity of the environment is governed by the
spatio-temporal dynamics and Pattern formation presumes as an attractive
proposition of the environment [46]. In ecosystem, it the general tendency of
prey is that it tries to stay away from the predators and therefore the escape
velocity of prey may be considered proportional to the diffusive velocity of
the predator population [47]. The significance of spatio-temporal models is
one of the leading issues in both ecology and mathematical ecology[48]. The
pattern formation in reaction–diffusion system has long been the subject of
interest to the researchers in the domain of mathematical ecology because of
its universal existence and importance [49]. Spatiotemporal dynamics of a
prey-predator model was studied by Lui and Kang [50] to investigate the impact
of fear effect and found that with increasing the cost of fear, the density of
predator population decreases at the positive steady state. The effect of fear
factor on system dynamics of diffusive prey-predator model was investigated by
Chen et al. [51] and authors suggested the proper range of parameters at which
the spatiotemporal pattern can occur. Han et al. [52] found that the fear
factor has great influence on spatially inhomogeneous distribution of the two
species under certain cross-diffusivity in a diffusive modified Leslie-Gower
prey-predator model.

Wang et al. [39] and Zhang et al. [43], in their studies, have considered
Holling type II and Holling type III functional responses, respectively. They
obtained that a high level of fear has a stabilising effect on the system by
removing the periodic solutions. As a result, it could be interesting to
investigate whether the fear effect still has a stabilizing effect on the
system in the presence of a Holling type IV functional response (group
defence). In this paper, we study a spatial and its associated non spatial
fear effect induced prey-predator interaction system in the presence of prey
refuge and group defence. Our study objectives include:

  * •

What are the dynamic repercussions of fear effect in presence of group
defence?

  * •

In the presence of group defence, how the prey birth rate and intra-prey
competition affect the dynamics of the model system ?

  * •

To investigate different type of bifurcations (e.g., saddle-node,
transcritical, Hopf) and stability results

  * •

What is the impact of prey refuge and fear level parameters on the spatial
prey-predator system?

The paper is organised in the following manner: The construction of our model
is covered in Section 2, which explains how fear and group defence are related
through predator-taxis sensitivity while incorporating fear effect, refuge,
and non-monotone group defence in the model. Section 3 provides a
comprehensive equilibrium analysis and associated stability analysis. In this
section, we also go through a few additional aspects of the suggested model,
such as species persistence and extinction. Using numerical simulation and
one- and two-parameter bifurcation diagrams, in Section 4, we discuss all
possible local bifurcations, including transcritical, saddle-node, and Hopf
bifurcations. Section 5 discusses the effect of prey birth rate and fear level
on prey and predator equilibrium densities. In section, 6 spatial model system
and its dynamics has been discusses. In section 7, we numerically verify the
conclusions of our analysis. The results and conclusions are discussed in
Section 8.

##  2 Model Formulation:

Assume that we have a closed habitat with two species, namely prey and
predator with population densities of $x$ and $y$, respectively. We consider
intraspecies competition in both prey and predator as it is a very common
phenomenon in real life and bound to happen. The following ordinary
differential equation can be used to describe the dynamics of prey populations
in the absence of predator species:

$\dfrac{dx}{dt}=rx-d_{1}x-d_{2}x^{2},$ (1)

where $r$ is the prey’s birth rate, $d_{1}$ is prey’s natural death rate, and
$d_{2}$ is prey’s death rate caused by intra-prey competition. From Eq.1 It is
straightforward to demonstrate that $x(t)$ tends to $0$ whenever $r<d_{1}$ for
all positive initial conditions, i.e. the prey population goes to extinct
under this assumption. Therefore, we suppose that $r>d_{1}$.

In the presence of predator species, the population of prey species decreases
due to predation by predators. When predators graze on prey species, the
Holling type IV functional response $\left(i.e.~{}~{}p(x)=\frac{\beta
x}{a+bx+x^{2}}\right)$ serves as the group defence mechanism or anti-predation
response. Thus, the following system of differential equations describes the
dynamics of the prey and predator interaction:

$\displaystyle\dfrac{dx}{dt}$ $\displaystyle=rx-
d_{1}x-d_{2}x^{2}-p\left(x\right)y,$ (2) $\displaystyle\dfrac{dy}{dt}$
$\displaystyle=-dy-ey^{2}+cp\left(x\right)y,$

where, $d$ is predator’s natural death rate, $e$ is predator’s death rate
caused by intra-predator competition, $c$ is the conversion efficiency of prey
biomass, $p(x)$ is the group defense function, $\beta$ is the predation rate,
$a$ is the half saturation constant, and $b$ is the tolerance limit (defence
level) of predator.

As previously stated, predators not only kill but also induce fear in the prey
population, which can be more lethal than direct killing in the sense of a
reduction in the prey numbers. Therefore, in order to model this reduction in
prey numbers, we introduce a fear factor $f\left(k,y\right)$, which represents
the cost of fear on prey’s birth rate. The fear of prey also results in
antipredator strategies such as vigilance, group defense, etc. But we also
know that the time and energy of a system are limited and constant, which
means if prey spends more time in group defence or other antipredator
strategies, then it may lead to a lower production rate of prey due to reduced
time for production. Because of this, it’s crucial to interconnect these two
aspects (fear factor and group defense) using a term called “predator taxis
sensitivity” (which represents the sensitiveness of prey to predation) and it
is termed as ‘$\alpha$’. If “predator-taxis sensitivity” increases, then group
defence among prey species will increase, which eventually decreases the
predation rate. Also, as “predator-taxis sensitivity” increases, the value of
the fear factor decreases as more sensitive prey will be more fearful (because
their level of fear will increase) and that eventually results in a reduction
in reproduction rate. In our model, we assume that both the tolerance limit of
predators ‘b’ and the fear level ‘k’ increase linearly with respect to
‘$\alpha$’ $~{}(i.e.,~{}b\to b\alpha,k\to k\alpha)$. therefore, the fear
factor mentioned above $\left(i.e.~{}f\left(k,y\right)\right)$ is now modified
into $f\left(k,\alpha,y\right)$ and the group defence function
$p\left(x\right)$ is modified in to $p\left(\alpha,x\right)$, where

$\displaystyle p\left(\alpha,x\right)=\frac{\beta x}{a+b\alpha x+x^{2}}.$ (3)

Ecologically, as Zanette [4]and some other field experiments [3, 2] showed
that fear causes a reduction in production rate, the fear factor and group
defence function should satisfy the following conditions:

  1. (i)

$f\left(0,\alpha,y\right)=1$: in the absence of fear of predation, the prey
population grows at its birth rate;

  2. (ii)

$f\left(k,\alpha,0\right)=1$: in the absence of predators, there is no fear of
predation in the prey population. Therefore, the fear factor does not affect
the birth rate of prey individuals;

  3. (iii)

$f\left(k,0,y\right)=1$: when prey species are not sensitive enough towards
predation, the birth rate of prey individuals remains unaffected by the fear
factor because they don’t feel any kind of fear from predators;

  4. (iv)

$\lim\limits_{k\to\infty}f\left(k,\alpha,y\right)=0$: when the level of fear
is very high, prey species show high antipredator behaviour that leads the
prey’s birth rate towards zero;

  5. (v)

$\lim\limits_{y\to\infty}f\left(k,\alpha,y\right)=0$: a very large predator
population size triggers a high level of fear in the prey population and leads
to a zero prey birth rate;

  6. (vi)

$\lim\limits_{\alpha\to\infty}f\left(k,\alpha,y\right)=0$: highly sensitive
prey population becomes more fearful, and that leads to the fear factor
approaching zero, which ultimately makes the prey birth rate zero;

  7. (vii)

$\frac{\partial f\left(k,\alpha,y\right)}{\partial k}<0$: the rate of prey
production decreases as the level of fear or antipredator behaviour increases;

  8. (viii)

$\frac{\partial f\left(k,\alpha,y\right)}{\partial y}<0$: prey production rate
decreases as predator population size increases;

  9. (ix)

$\frac{\partial f\left(k,\alpha,y\right)}{\partial\alpha}<0$: as the predator
taxis-sensitivity increases, the prey birth rate decreases;

  10. (x)

$\lim\limits_{\alpha\to\infty}p\left(\alpha,x\right)=0$: when prey species are
highly sensitive towards predation, their group defence becomes stronger and
that leads to a zero predation rate;

  11. (xi)

$\frac{\partial p(x,\alpha)}{\partial\alpha}<0$: when the prey population
becomes highly sensitive, they make their group defence a lot stronger, which
leads to a reduction in predation rates.

Therefore $f\left(k,\alpha,y\right)$ can be taken as $\frac{1}{1+k\alpha
y}~{}or~{}e^{-k\alpha y}$ or some other function. Some scholars and authors
[30], in particular, have used the following form for the fear factor term:

$\displaystyle f\left(k,\alpha,y\right)=\frac{1}{1+k\alpha y}.$ (4)

Based on the preceding discussions, after incorporating the fear factor from
Eq. (4) in Eq. (2), our model becomes

$\displaystyle\dfrac{dx}{dt}$
$\displaystyle=rx\left(f\left(k,\alpha,y\right)\right)-d_{1}x-d_{2}x^{2}-p\left(\alpha,x\right)y,$
(5) $\displaystyle\dfrac{dy}{dt}$ $\displaystyle=-dy-
ey^{2}+cp\left(\alpha,x\right)y.$

Now, in order to be more realistic, it is essential to take into account the
system’s prey refuge. In reality, not all the prey population is exposed to
predators as they often have refuges where they can avoid predators. As a
result, we assume that a portion $mx$ (where $m$ ($0\leq m<1$) is the refuge
parameter) of prey species is completely protected from predators (via
refuges) and that only a portion $\left(1-m\right)x$ of prey species is
available to predators. By incorporating the refuge $mx$ and
$p\left(\alpha,x\right)$ from Eq. (3) into our system, we obtain:

$\displaystyle\dfrac{dx}{dt}$ $\displaystyle=\frac{rx}{1+k\alpha
y}-d_{1}x-d_{2}x^{2}-\frac{\left(-m+1\right)\beta
xy}{a+\left(-m+1\right)\alpha bx+\left(-m+1\right)^{2}x^{2}}~{},$ (6)
$\displaystyle\dfrac{dy}{dt}$ $\displaystyle=-dy-
ey^{2}+\frac{\left(-m+1\right)\beta cxy}{a+\left(-m+1\right)\alpha
bx+\left(-m+1\right)^{2}x^{2}}.$

Table 1 represents the system parameters and their biological meaning. Two
sets of values for these parameters have also been provided, which we have
used throughout this paper.

Parameters | Biological meaning | values (set 1) | values (set 2)
---|---|---|---
$r$ | prey’s birth rate | 1.1 | 1.1
$d_{1}$ | prey’s natural death rate | 0.05 | 0.2
$d_{2}$ | prey’s death rate caused by intra-prey competition | 0.05 | 0.1
$\beta$ | predation rate | 0.071 | 0.241
$m$ | prey refuge parameter | 0.3 | 0.3
$k$ | the level of fear | 0.5 | 8
$\alpha$ | predator-taxis sensitive | 1 | 0.25
$d$ | predator’s natural death | 0.05 | 0.8
$e$ | predator’s death rate caused by intra-predator competition | 0.05 | 0.5
$c$ | conversion efficiency of prey biomass | 10 | 17
$a$ | the half saturation constant | 0.1 | 0.9
$b$ | the tolerance limit of predators | 0.125 | 0.5
Table 1: Biological meaning of model parameters for system (6)

##  3 Mathematical Analysis

In this section we shall discuss the positivity and boundedness of the
solutions of the model system (6). We shall also discuss the persistence and
extinction scenarios, local stability and bifurcation analysis.

###  3.1 Positivity and Boundedeness

It is important to show the positivity of solutions in the system as it
implies that the population survives, and once these solutions enter the
positive quadrant, they will remain in it forever. So the solutions to the
system will be biologically feasible if they lie in the positive quadrant. The
boundedness of the system variables defining populations is vital to show. It
can be seen as a natural limit to the growth of these system variables as a
result of limited resources. Because of this, the model’s boundedness shows
that it is biologically appropriate.

######  Theorem 1.

For the proposed model system (6), $\mathbb{R}_{+}^{2}$ is positively
invariant. Furthermore, any solution of the model system (6) is bounded within
this region $\mathbb{R}_{+}^{2}$ and satisfies the following condition:
(i) $\limsup\limits_{t\to+\infty}x(t)\leq\left(\frac{r-d_{1}}{d_{2}}\right)$;
provided $r>d_{1}$, (ii)
$\limsup\limits_{t\to+\infty}y(t)\leq\left(\frac{c\beta-bd\alpha}{b\alpha
e}\right)$; provided $c\beta>bd\alpha$.

###### Proof.

For the proof, refer to Appendix A.1. ∎

###  3.2 Extinction

In biology, extinction refers to the complete annihilation of a species.
Extinction happens when a species population is reduced due to environmental
factors (habitat fission, natural disasters, over-exploitatage of species) or
evolutionary changes in its members (genetic inbreeding, poor reproduction).
In this section, we’ll discuss about certain parametric criteria that, over a
long span of time, guarantee the extinction of a species.

Let $\overline{x}$=$\limsup\limits_{t\to\infty}x(t)$,
$\overline{y}$=$\limsup\limits_{t\to\infty}y(t)$,
$\underline{x}$=$\liminf\limits_{t\to\infty}x(t)$ and
$\underline{y}$=$\liminf\limits_{t\to\infty}y(t).$

######  Theorem 2.

  1. (i)

(Extinction of prey in the absence of predators) If
$r<d_{1},then\lim\limits_{t\to\infty}x\left(t\right)=0.$

  2. (ii)

(Extinction of predator due to prey population) If
$add_{2}>\left(-m+1\right)\beta c\left(r-d_{1}\right)$, then
$\lim\limits_{t\to\infty}y(t)=0.$

###### Proof.

  1. (i)

If possible, let $\lim\limits_{t\to\infty}x(t)=L>0$, then using the first
equation of the model system (6), we get

$\displaystyle\dfrac{dx}{dt}$ $\displaystyle\leq\frac{rx}{1+k\alpha
y}-d_{1}x-d_{2}x^{2}\leq rx-d_{1}x-d_{2}x^{2}\leq x\left(r-d_{1}\right),$
$\displaystyle x(t)$
$\displaystyle=x(0)exp\left(\int_{0}^{t}\left(r-d_{1}\right)\,ds\right).$

Now if $~{}r<d_{1}$ then $\left(r-d_{1}\right)<0$, hence

$\displaystyle x(t)$
$\displaystyle=x(0)exp\left(\int_{0}^{t}\left(r-d_{1}\right)\,ds\right)\to
0~{}~{}as~{}~{}t\to\infty,$

which is a contradiction to our assumption. Hence, our theorem is proved.

Remark: According to this theorem, in the absence of predators, the prey
population goes extinct since natural deaths are sufficient to eliminate the
prey population from the system. In the considered model system, the growth of
the predator population entirely depends on the prey population. Therefore,
when the prey’s natural death rate exceeds the birth rate of the prey, both
the prey and predator populations are wiped out from the system. Henceforth,
we consider that the prey birth rate is greater than their natural death rate
$(r>d_{1})$ unless otherwise stated.

  2. (ii)

From the second equation of the model system(6), we have

$\displaystyle\frac{dy}{dt}$ $\displaystyle=-dy-ey^{2}+\frac{(-m+1)\beta
cy}{a}\left(x-\frac{(-m+1)\alpha bx^{2}+(-m+1)^{2}x^{3}}{a+(-m+1)\alpha
bx+(-m+1)^{2}x^{2}}\right),$ $\displaystyle\leq y\left(-d-ey+\frac{(-m+1)\beta
cx}{a}\right)\leq y\left(-d+\frac{(-m+1)\beta
c(\frac{r-d_{1}}{d_{2}})}{a}\right),$ $\displaystyle\int_{0}^{t}\
\frac{dy}{y}$
$\displaystyle\leq\int_{0}^{t}\left(-d+\frac{\left(-m+1\right)\beta
c\left(\frac{r-d_{1}}{d_{2}}\right)}{a}\right)\,ds,$ $\displaystyle
y\left(t\right)$ $\displaystyle\leq
y\left(0\right)\exp\left(\int_{0}^{t}\left(-d+\frac{\left(-m+1\right)\beta
c\left(\frac{r-d_{1}}{d_{2}}\right)}{a}\right)\,ds\right).$

So, if $add_{2}>\left(-m+1\right)\beta c\left(r-d_{1}\right)$, then
$y\left(t\right)\to 0~{}~{}as~{}~{}t\to 0$.
From the second equation of the model system(6), we have

$\displaystyle\dfrac{dy}{dt}$ $\displaystyle\leq-
dy+\frac{\left(-m+1\right)\beta cxy}{a+\left(-m+1\right)\alpha
bx}=y\left(\frac{-d\left(a+\left(m-1\right)\alpha bx\right)+\beta
cx\left(1-m\right)}{a+(1-m)\alpha bx}\right),$ $\displaystyle\dfrac{dy}{dt}$
$\displaystyle\leq y\left(\frac{-da+(1-m)x(c\beta-d\alpha b)}{a+(1-m)\alpha
bx}\right).$

So, if $c\beta<d\alpha b$, then $y\left(t\right)\to 0$ as $t\to 0$.

∎

###  3.3 Persistence of Species

A species is ecologically persistent if it is never in danger of extinction or
is never entirely eradicated from the system. The ability of a species to
persist means that it can withstand exogenous perturbations. As a result, it
is a significant system behaviour in the long run. Analytically, we can define
the persistence of species through the following theorem:

######  Theorem 3.

  1. (i)

The prey species is strongly persistent under the following condition:
$\frac{rbe}{be+k\left(c\beta-
bd\alpha\right)}>d_{1}+\frac{\beta\left(-m+1\right)\left(c\beta-
bd\alpha\right)}{abe\alpha}.$

  2. (ii)

The predator species is strongly persistent under the following conditions:
$r>d_{1}$ and
$c\beta\left(1-m\right)\left(\frac{r-d_{1}}{d_{2}}\right)>d\left(a+b\alpha\left(1-m\right)\left(\frac{r-d_{1}}{d_{2}}\right)+\left(1-m\right)^{2}\left(\frac{r-d_{1}}{d_{2}}\right)^{2}\right).$

###### Proof.

  1. (i)

Using upper bound of $y$ from Theorem 1 in the first equation of the model
system(6) to obtain

$\displaystyle\frac{dx}{dt}$ $\displaystyle\geq
x\left(\left(\frac{r}{1+k\alpha\left(\frac{c\beta-bd\alpha}{b\alpha
e}\right)}\right)-d_{1}-d_{2}x-\frac{\left(-m+1\right)\beta\left(\frac{c\beta-
bd\alpha}{b\alpha e}\right)}{a+\left(-m+1\right)\alpha
bx+\left(-m+1\right)^{2}x^{2}}\right),$

since $\left(a+\left(-m+1\right)\alpha
bx+\left(-m+1\right)^{2}x^{2}\right)\geq a$, then from the above expression we
get:

$\displaystyle\frac{dx}{dt}$ $\displaystyle\geq
x\left(\left(\frac{r}{1+k\alpha\left(\frac{c\beta-bd\alpha}{b\alpha
e}\right)}\right)-d_{1}-d_{2}x-\frac{\left(-m+1\right)\beta\left(\frac{c\beta-
bd\alpha}{b\alpha e}\right)}{a}\right).$

Now by comparison lemma ([53]), we obtain

$\displaystyle\liminf\limits_{t\to+\infty}x(t)\geq\frac{rbe}{bed_{2}+kd_{2}\left(c\beta-
bd\alpha\right)}-\frac{d_{1}}{d_{2}}-\frac{\beta\left(-m+1\right)\left(c\beta-
bd\alpha\right)}{abed_{2}\alpha};$ $\displaystyle
provided:~{}~{}\frac{rbe}{be+k\left(c\beta-
bd\alpha\right)}>d_{1}+\frac{\beta\left(-m+1\right)\left(c\beta-
bd\alpha\right)}{abe\alpha}.$

  2. (ii)

Using upper bound of $x$ from Theorem 1 in the first equation of the model
system(6) to obtain

$\displaystyle\frac{dy}{dt}$ $\displaystyle\geq
y\left(-d-ey+\frac{c\beta\left(1-m\right)\left(\frac{r-d_{1}}{d_{2}}\right)}{a+b\alpha\left(1-m\right)\left(\frac{r-d_{1}}{d_{2}}\right)+\left(1-m\right)^{2}\left(\frac{r-d_{1}}{d_{2}}\right)^{2}}\right).$

Now by comparison lemma ([53]), we obtain

$\displaystyle\liminf\limits_{t\to+\infty}y(t)\geq\left(\frac{-d}{e}+\frac{c\beta\left(1-m\right)\left(\frac{r-d_{1}}{d_{2}}\right)}{e\left(a+b\alpha\left(1-m\right)\left(\frac{r-d_{1}}{d_{2}}\right)+\left(1-m\right)^{2}\left(\frac{r-d_{1}}{d_{2}}\right)^{2}\right)}\right),$
$\displaystyle
provided:~{}~{}r>d_{1},\,\,c\beta\left(1-m\right)\left(\frac{r-d_{1}}{d_{2}}\right)>d\left(a+b\alpha\left(1-m\right)\left(\frac{r-d_{1}}{d_{2}}\right)+\left(1-m\right)^{2}\left(\frac{r-d_{1}}{d_{2}}\right)^{2}\right).$

∎

###  3.4 Existence of Equilibrium Points

This subsection investigates the parametric conditions for the presence of
several different feasible equilibrium points. By analysing the eigenvalues of
the related Jacobian matrices, we also investigate the local stability of
these equilibrium points. We consider the model system (6) in order to
determine feasible equilibrium points:

$\displaystyle\frac{dx}{dt}$ $\displaystyle=\frac{rx}{1+k\alpha
y}-d_{1}x-d_{2}x^{2}-\frac{\left(-m+1\right)\beta
xy}{a+\left(-m+1\right)\alpha
bx+\left(-m+1\right)^{2}x^{2}}\equiv~{}g\left(x,y\right),$
$\displaystyle\frac{dy}{dt}$ $\displaystyle=-dy-
ey^{2}+\frac{\left(-m+1\right)c\beta xy}{a+\left(-m+1\right)\alpha
bx+\left(-m+1\right)^{2}x^{2}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\equiv~{}h\left(x,y\right).$

Ecologically, equilibrium points with non-negative co-ordinates only, i.e.,
the points of intersection of $x$ and $y$ growth nullclines in
$\mathbb{R}_{+}^{2}\cup\left(0,0\right)$, are feasible. $\left(\
i.e.,\left(x,y\right)\in\mathbb{R}_{+}^{2}\cup\left(0,0\right)~{}:~{}g\left(x,y\right)=0,h\left(x,y\right)=0\right)$.
Therefore, the model system (6) does have the following feasible equilibrium
point:-

  1. (i)

The population extinct equilibrium point $E_{0}(0,0)$, which always exists.

  2. (ii)

The predator extinct equilibrium point $E_{1}(\frac{r-d_{1}}{d_{2}},0),$ which
exists if $r>d_{1}$.

  3. (iii)

An interior (coexisting) equilibrium point $E^{*}=(x^{*},y^{*})$ exists in
addition to the trivial and boundary equilibrium points if $d<\frac{c\beta
x^{*}(1-m)}{a+b\alpha x^{*}(1-m)+(1-m)^{2}(x^{*})^{2}}$. For
$E^{*}=(x^{*},y^{*})$, the predator species’ component is
$y^{*}=\frac{1}{e}\left(-d+\frac{c\beta x^{*}(1-m)}{a+b\alpha
x^{*}(1-m)+(1-m)^{2}(x^{*})^{2}}\right)$ and the prey component $x^{*}$is a
positive root of the following equation:

$\displaystyle
p(x)=A_{1}x^{7}+A_{2}x^{6}+A_{3}x^{5}+A_{4}x^{4}+A_{5}x^{3}+A_{6}x^{2}+A_{7}x+A_{8}.$
(7)

Appendix B.1 contains the coefficients $A_{1}$, $A_{2}$, $A_{3}$, $A_{4}$,
$A_{5}$, $A_{6}$, $A_{7}$, and $A_{8}$. It is difficult to find the analytical
conditions to ascertain the actual number of interior (coexisting) equilibrium
points of the system (6). However, through the relative position and shape of
coexisting nullclines, the possible number of coexisting equilibrium points
can be determined. To find a different number of equilibrium points, we vary
$d_{2}$ and assign values to other model parameters from set (1) of the Table
(1). As shown in the Fig.1, predator nullcline (red curve) intersects prey
nullclines once, twice, and thrice for $d_{2}=0.05$ (black curve),
$d_{2}=0.055$ (green curve), and $d_{2}=0.09$ (cyan curve), respectively. So,
the model system (6) can have a maximum of three coexisting equilibrium points
for the chosen set of parametric values. Table 2 shows the interior
equilibrium points and corresponding local stability for different values of
$d_{2}$.

###  3.5 Local Stability Analysis

######  Theorem 4.

  1. (i)

The population extinct equilibrium point $E_{0}(0,0)$ is locally
asymptotically stable if $(-d_{1}+r)<0$. In particular, under this condition,
$E_{0}(0,0)$ is globally asymptotically stable.

  2. (ii)

The predator extinct equilibrium $E_{1}(\frac{r-d_{1}}{d_{2}},0)$ is locally
asymptotically stable if $r>d_{1}$, and $\left[d(ad_{2}^{2}+\right.\\\
\left.(-m+1)\alpha b(r-d_{1})d_{2}+(-m+1)^{2}(r-d_{1})^{2})>(-m+1)\beta
c(r-d_{1})d_{2}\right]$.

  3. (iii)

The interior equilibrium point $E^{*}(x^{*},y^{*})$ is locally asymptotically
stable if $Tr(J(E^{*})<0$ and $Det(J(E^{*})>0$, unstable otherwise.

(a)  Figure 1: The red curve represents the predator nullcline for all values
of $d_{2}$. The black curve, the green curve, the magenta curve and the cyan
curve represent the prey nullcline for
$d_{2}=0.05,d_{2}=0.055,d_{2}=0.077854268968,d_{2}=0.09$, respectively.

###### Proof.

  1. (i)

For the proof, refer to Appendix A.2.

  2. (ii)

For the proof, refer to Appendix A.2.

  3. (iii)

The Jacobian matrix of model system (6) at $E^{*}(x^{*},y^{*})$ is given by

$J({E^{*}})={(J)_{{E^{*}}}}=\left({\begin{array}[]{*{20}{c}}{{g_{x}}}&{{g_{y}}}\\\
{}\hfil&{}\hfil\\\ {{h_{x}}}&{{h_{y}}}\\\
\end{array}}\right){|_{{E^{*}}}}=\left({\begin{array}[]{*{20}{c}}{{j_{11}}}&{{j_{12}}}\\\
{}\hfil&{}\hfil\\\ {{j_{21}}}&{{j_{22}}}\\\ \end{array}}\right).$ (8)

Therefore, the characteristic equation of $J(E^{*})$ is

$\lambda^{2}-Tr(J(E^{*}))\lambda+Det(J(E^{*}))=0,$

where, $Tr(J(E^{*}))=j_{11}+j_{22}$ and
$Det(J(E^{*}))=j_{11}j_{22}-j_{12}j_{21}$. Thus, $E^{*}(x^{*},y^{*})$ is
locally asymptotically stable if $Tr(J(E^{*}))<0$ and $Det(J(E^{*}))>0$.

∎

Finding the analytical conditions to determine the actual number of coexisting
equilibrium points is difficult for the system (6). Therefore, here in the
Table 2, we have numerically presented various numbers of equilibrium points
and discussed their corresponding local stability for a chosen set of
parametric values.

S.No. | $d_{2}$ | No. of |  |  |  |
---|---|---|---|---|---|---
|  | E.P | E.P(s) | $\lambda_{1}$ | $\lambda_{2}$ | Stability
1 | 0.05 | 2 | $E_{2}^{*}=(2.5356379,6.2586513)$ | -0.4883204 | 0.2169094 | S.P
$E_{3}^{*}=(0.0234372,1.2768451)$ |  -0.024659 + $\iota$ 0.310341 |  -0.024659 - $\iota$ 0.310341 | S.S
2 | 0.055 | 3 | $E_{1}^{*}=(18.0018467,0.1151048)$ | -0.9931931 | -0.0025924 | S.N
$E_{2}^{*}=(2.6405682,6.0037036)$ | -0.4820788 | 0.1867656 | S.P
$E_{3}^{*}=(0.0234354,1.2766738)$ |  -0.024715 + $\iota$ 0.310322 |  -0.024715 - $\iota$ 0.310322 | S.S
3 | 0.07785426896887 | 2 | $E_{12}^{*}=(4.78094,3.05535)$ | -0.49974307 | -0.00000226 | S.N
$E_{3}^{*}=(0.0234272,1.2758912)$ |  -0.024971+ $\iota$ 0.3102333 |  -0.024971- $\iota$ 0.3102333 | S.S
4 | 0.09 | 1 | $E_{3}^{*}=(0.0234228,1.2754756)$ |  -0.025107+ $\iota$ 0.3101861 |  -0.025107- $\iota$ 0.3101861 | S.S
5 | $d_{2}=0.1,d_{1}=1.2$ | 0 | - | - | - | -
Table 2:  Other than $d_{2}$, all other parametric values are taken from the
set (1) of the Table 1. Interior equilibrium point(s) and corresponding local
stability for different values of $d_{2}$. E.P, S.P, S.S, and S.N are
abbreviations for equilibrium point, saddle point, stable spiral, and stable
node, respectively.

##  4 Local bifurcation analysis

In this section, we shell discuss about how the different parameter values
affect the stability of the different equilibrium points. In addition, we
shell also investigate that how the equilibrium points arise/ dissappear.

###  4.1 Saddle-node bifurcation

######  Theorem 5.

The model system (6) undergoes a saddle-node bifurcation around the interior
equilibrium point $E^{*}(x^{*},y^{*})$ as the intra-prey competition parameter
$d_{2}$ passes through the bifurcation value $d_{2}=(d_{2})_{SN}$ if and only
if

$\displaystyle\left[g_{xx}-\frac{g_{x}}{g_{y}}(g_{xy}+g_{yx})+\left(\frac{g_{x}}{g_{y}}\right)^{2}g_{yy}-\left(\frac{g_{x}}{h_{x}}\right)\left(h_{xx}-\frac{g_{x}}{g_{y}}(h_{xy}+h_{yx})+\left(\frac{g_{x}}{g_{y}}\right)^{2}h_{yy}\right)\right]\Bigg{|}_{((d_{2})_{SN},E^{*})}\neq
0.$

###### Proof.

Let $d_{2}=(d_{2})_{SN}$ and $E^{*}(x^{*},y^{*})$ be the interior equilibrium
point. Define:

$f(x,y)=\left(\begin{array}[]{c}g(x,y)\\\ \\\
h(x,y)\end{array}\right)=\left(\begin{array}[]{c}\frac{rx}{1+k\alpha
y}-d_{1}x-d_{2}x^{2}-\frac{\left(-m+1\right)\beta
xy}{a+\left(-m+1\right)\alpha bx+x^{2}\left(-m+1\right)^{2}}\\\ \\\ -dy-
ey^{2}-\frac{c\beta\left(-m+1\right)xy}{a+b\alpha\left(-m+1\right)x+\left(-m+1\right)^{2}x^{2}}\end{array}\right),$

differentiating $f$ with respect to $d_{2}$, we obtain

$\displaystyle f_{d_{2}}(E^{*})=\left(\begin{array}[]{c}g_{d_{2}}\\\ \\\
h_{d_{2}}\end{array}\right)\Bigg{|}_{(x^{*},y^{*})}=\left(\begin{array}[]{c}-(x^{*})^{2}\\\
\\\ 0\end{array}\right).$

Now we evaluate

$\displaystyle
S=Df(E^{*},(d_{2})_{SN})=\left(\begin{array}[]{cc}g_{x}&g_{y}\\\ \\\
h_{x}&h_{y}\end{array}\right)\Bigg{|}_{(E^{*},(d_{2})_{SN}).}$

So, $V=\left(1~{}~{}\frac{-g_{x}}{g_{y}}\right)^{T}$ and
$W=\left(1~{}~{}\frac{-g_{x}}{h_{x}}\right)^{T}$ are the eigenvectors
corresponding to eigenvalue $\lambda=0$ of $S$ and $S^{T}$, respectively. we
also evaluate

$\displaystyle\left[D^{2}f(E^{*},(d_{2})_{SN})(V,V)\right]=\left(\begin{array}[]{c}g_{xx}-\frac{g_{x}}{g_{y}}(g_{xy}+g_{yx})+\frac{g^{2}_{x}}{g^{2}_{y}}g_{yy}\\\
\\\
h_{xx}-\frac{g_{x}}{g_{y}}(h_{xy}+h_{yx})+\frac{g^{2}_{x}}{g^{2}_{y}}h_{yy}\end{array}\right).$

The differentiations of functions ‘f’ and ‘g’ have been provided in Appendix
B.2. Now using the Sotomayor theorem [54] for saddle-node bifurcation, we
obtain

$\displaystyle W^{T}\left[f_{d_{2}}(E^{*})\right]=-(x^{*})^{2}\neq 0,$

and

$\displaystyle W^{T}\left[D^{2}f(E^{*},(d_{2})_{SN})(V,V)\right]=$
$\displaystyle\left(g_{xx}-\frac{g_{x}}{g_{y}}(g_{xy}+g_{yx})+\frac{g^{2}_{x}}{g^{2}_{y}}g_{yy}-\frac{g_{x}}{h_{x}}\left(h_{xx}-\frac{g_{x}}{g_{y}}(h_{xy}+h_{yx})+\right.\right.$
$\displaystyle\left.\left.\frac{g^{2}_{x}}{g^{2}_{y}}h_{yy}\right)\right)\Bigg{|}_{((d_{2})_{SN},E^{*}).}$

Therefore, system (6) undergoes a saddle-node bifurcation at
$E^{*}(x^{*},y^{*})$ if and only if

$\displaystyle\left[g_{xx}-\frac{g_{x}}{g_{y}}(g_{xy}+g_{yx})+\frac{g^{2}_{x}}{g^{2}_{y}}g_{yy}-\frac{g_{x}}{h_{x}}\left(h_{xx}-\frac{g_{x}}{g_{y}}(h_{xy}+h_{yx})+\frac{g^{2}_{x}}{g^{2}_{y}}h_{yy}\right)\right]\Bigg{|}_{((d_{2})_{SN},E^{*})}\neq
0.$

∎

Analytical demonstration of saddle-node bifurcation for interior equilibrium
point is quite challenging. However, we provide a numerical example to verify
the existence of saddle-node bifurcation for the system (6). Here we verify
the saddle-node bifurcation for $d_{2}=(d_{2})_{SN}=0.0778542689688737$. The
other parameter values are taken from set (1) of Table 1. The interior
equilibrium points $E^{*}_{1}$ and $E^{*}_{2}$ coincide at interior
equilibrium point $E^{*}_{12}=(4.7809,3.0553)$ and mutually annihilate via
saddle-node bifurcation as shown in Fig. 3b.

######  Remark 6.

We also observe a saddle-node bifurcation for
$r=(r)_{SN}=3.27311875453723710727160$ at interior equilibrium point
$E^{*}_{12}=(1.247196,2.438599)$ of the system (6) where two interior
equilibrium points $E^{*}_{1}$ and $E^{*}_{2}$ coincide and mutually
annihilate as depicted in Fig. 4.

###  4.2 Transcritical bifurcation

######  Theorem 7.

The model system (6) undergoes a transcritical bifurcation from the population
extinct equilibrium point $E_{0}=(0,0)$ when $r=r_{TC}=d_{1}.$

###### Proof.

The jacobian matrix at $E_{0}=(0,0)$ and $r=d_{1}$ is given by

$J(E_{0},r_{TC})=\left(\begin{array}[]{cc}0&0\\\ 0&-d\\\ \end{array}\right),$

whose eigenvalues are $\lambda_{1}=0$ and $\lambda_{2}=-d$. Now we choose the
eigenvectors V and W associated with eigen value $\lambda_{1}$ of
$J(E_{0},d_{1})$ and $J(E_{0},d_{1})^{T}$ given respectively by
$V=\left(\begin{array}[]{c}v_{1}\\\
v_{2}\end{array}\right)=\left(\begin{array}[]{c}1\\\ 0\end{array}\right)$ and
$W=\left(\begin{array}[]{c}w_{1}\\\
w_{2}\end{array}\right)=\left(\begin{array}[]{c}1\\\ 0\end{array}\right).$
Recall $f$ from the Theorem 5 and calculate:

$\displaystyle\Delta_{1}$
$\displaystyle=W^{T}f_{r}\left(E_{0},r=r_{TC}\right)=\left(\begin{array}[]{cc}1&0\end{array}\right)\cdot\left(\begin{array}[]{c}\frac{\partial
g}{\partial r}\\\ \\\ \frac{\partial h}{\partial
r}\end{array}\right)_{\left({E_{0}};r_{TC}\right)}=\left(\begin{array}[]{cc}1&0\end{array}\right)\cdot\left(\begin{array}[]{c}0\\\
0\end{array}\right)=0,$
$\displaystyle\Delta_{2}=W^{T}\left[Df_{r}\left(E_{0},r=r_{TC}\right)V\right]=$
$\displaystyle\left(\begin{array}[]{cc}1&0\end{array}\right)\cdot\left(\begin{array}[]{cc}\frac{\partial^{2}g}{\partial
x\partial r}&\frac{\partial^{2}g}{\partial y\partial r}\\\ \\\
\frac{\partial^{2}h}{\partial x\partial r}&\frac{\partial^{2}h}{\partial
y\partial
r}\end{array}\right)_{\left(E_{0};r_{TC}\right)}\cdot\left(\begin{array}[]{c}1\\\
0\end{array}\right)$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}1&0\end{array}\right)\cdot\left(\begin{array}[]{cc}1&0\\\
\\\ 0&0\end{array}\right)\cdot\left(\begin{array}[]{c}1\\\
0\end{array}\right)=1\neq 0,$
$\displaystyle\Delta_{3}=W^{T}\left[D^{2}f\left(E_{0},r=r_{TC}\right)(V,V)\right]=$
$\displaystyle\left(\begin{array}[]{cc}1&0\end{array}\right)\left(\begin{array}[]{c}\frac{\partial^{2}g}{\partial
x^{2}}v_{1}^{2}+\frac{2\partial^{2}g}{\partial x\partial
y}v_{1}v_{2}+\frac{\partial^{2}g}{\partial y^{2}}v_{2}^{2}\\\
\frac{\partial^{2}h}{\partial x^{2}}v_{1}^{2}+\frac{2\partial^{2}h}{\partial
x\partial y}v_{1}v_{2}+\frac{\partial^{2}h}{\partial
y^{2}}v_{2}^{2}\end{array}\right)\Bigg{|}_{(E_{0},r_{TC})}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}\frac{\partial^{2}g}{\partial
x^{2}}v_{1}^{2}+\frac{2\partial^{2}g}{\partial x\partial
y}v_{1}v_{2}+\frac{\partial^{2}g}{\partial
y^{2}}v_{2}^{2}\end{array}\right)\Bigg{|}_{(E_{0},r_{TC})}=\left(\begin{array}[]{c}\frac{\partial^{2}g}{\partial
x^{2}}\end{array}\right)\Bigg{|}_{(E_{0},r_{TC})}$ $\displaystyle=$
$\displaystyle-2d_{2}\neq 0.$

Hence, all the transversality condition are satisfied, therefore the model
system (6) undergoes a transcritical bifurcation from population extinct
equilibrium $E_{0}$ at $r=r_{TC}=d_{1}$.
∎

######  Remark 8.

When we choose the parametric values from set (2) of Table 1 then a
transcritical bifurcation for the system (6) occurs for $r=(r)_{TC}=0.2$ at
population extinct equilibrium point $E_{0}$ (refer to Fig. 4).

######  Theorem 9.

The model system (6) undergoes a transcritical bifurcation from predator
extinct equilibrium point $E_{1}(\frac{r-d_{1}}{d_{2}},0)$ at
$\beta=\beta_{TC}=\frac{d(ad_{2}^{2}+b\alpha(1-m)(r-d_{1})d_{2}+(1-m)^{2}(r-d_{1})^{2})}{c(1-m)(r-d_{1})d_{2}}$,
if ($ec(r-d_{1})^{2})\neq S$, where
$S=d((1-m)^{2}(r-d_{1})^{2}-ad_{2}^{2})(dd_{2}^{2}+ckr\alpha(r-d_{1}))$.

###### Proof.

The jacobian matrix at $E_{1}(\frac{r-d_{1}}{d_{2}},0)$ and
$\beta=\beta_{TC}$,

$J(E_{1},\beta_{TC})=\left(\begin{array}[]{cc}d_{1}-r&\frac{kr\alpha\left(d_{1}-r\right)}{d_{2}}-\frac{d}{c}\\\
\\\ 0&0\\\ \end{array}\right),$

whose eigenvalues are $\lambda_{1}=0$ and $\lambda_{2}=-r+d_{1}$. Now we
choose the eigenvectors V and W associated with eigen value $\lambda_{1}$
(zero eigen value) of $J(E_{1},\beta_{TC})$ and $J(E_{1},\beta_{TC})^{T}$
given respectively by $V=\left(\begin{array}[]{c}v_{1}\\\ \\\
v_{2}\end{array}\right)=\left(\begin{array}[]{c}\frac{-ckr^{2}\alpha+ckr\alpha
d_{1}-dd_{2}}{c(r-d_{1})d_{2}}\\\ \\\ 1\end{array}\right)$ and
$W=\left(\begin{array}[]{c}w_{1}\\\ \\\
w_{2}\end{array}\right)=\left(\begin{array}[]{c}0\\\ \\\ 1\end{array}\right)$.
Recall $f$ from Theorem 5 and calculate:

$\displaystyle f_{\beta}({E_{1},\beta_{TC}})$
$\displaystyle=\left(\begin{array}[]{c}\frac{(1-m)xy}{a+b\alpha(1-m)x+(1-m)^{2}x^{2}}\\\
\\\
\frac{c(1-m)xy}{a+b\alpha(1-m)x+(1-m)^{2}x^{2}}\end{array}\right)\Bigg{|}_{(E_{1},\beta_{TC})},~{}f_{\beta}({E_{1},\beta_{TC}})=\left(\begin{array}[]{c}0\\\
0\end{array}\right),$ $\displaystyle W^{T}[f_{\beta}({E_{1},\beta_{TC}})]$
$\displaystyle=\left(\begin{array}[]{cc}0&1\end{array}\right)\cdot\left(\begin{array}[]{c}0\\\
0\end{array}\right)=0.$

Now, we evaluate $W^{T}[Df_{\beta}(E_{1},\beta_{TC})V]$ and find

$\displaystyle W^{T}[Df_{\beta}(E_{1},\beta_{TC})V]$
$\displaystyle=\frac{cd_{2}(1-m)(r-d_{1})}{(ad_{2}^{2}+(1-m)^{2}(r-d_{1})^{2}+d_{2}b\alpha(1-m)(r-d_{1}))}\neq
0.$

Now, also we have

$\displaystyle W^{T}[D^{2}f(E_{1},\beta_{TC})(V,V)]=$
$\displaystyle\left(\begin{array}[]{c}\frac{\partial^{2}h}{\partial
x^{2}}v_{1}^{2}+\frac{2\partial^{2}h}{\partial x\partial
y}v_{1}v_{2}+\frac{\partial^{2}h}{\partial y^{2}}v_{2}^{2}\\\
\end{array}\right)\Bigg{|}_{(E_{1},\beta_{TC})},$

where

$\displaystyle\frac{\partial^{2}h}{\partial
x^{2}}\Bigg{|}_{(E_{1},\beta_{TC})}=$ $\displaystyle
0,\,\,\,\,\frac{\partial^{2}h}{\partial
y^{2}}\Bigg{|}_{(E_{1},\beta_{TC})}=-2e,$
$\displaystyle\frac{\partial^{2}h}{\partial x\partial
y}\Bigg{|}_{(E_{1},\beta_{TC})}=$
$\displaystyle\frac{-c(1-m)(\beta_{TC})d_{2}^{2}((m-1)^{2}(r-d_{1})^{2}-ad_{2}^{2})}{((m-1)^{2}(r-d_{1})^{2}+b\alpha
d_{2}(1-m)(r-d_{1})+ad_{2}^{2})}$ $\displaystyle=$
$\displaystyle\frac{-dd_{2}((m-1)^{2}(r-d_{1})^{2}-ad_{2}^{2})}{(r-d_{1})}.$

So, using the values of $\frac{\partial^{2}h}{\partial
x^{2}}\Bigg{|}_{(E_{1},\beta_{TC})},\frac{\partial^{2}h}{\partial
y^{2}}\Bigg{|}_{(E_{1},\beta_{TC})},\frac{\partial^{2}}{\partial x\partial
y}\Bigg{|}_{(E_{1},\beta_{TC})}$, $v_{1}$, and $v_{2}$, we get

$\displaystyle
W^{T}[D^{2}f(E_{1},\beta_{TC})(V,V)]=2\left(\frac{-ec(r-d_{1})^{2}+d((1-m)^{2}(r-d_{1})^{2}-ad_{2}^{2})(dd_{2}^{2}+ckr\alpha(r-d_{1}))}{c(r-d_{1})^{2}}\right)\neq
0.$

So, using the Sotomayor theorem [54], the theorem (9) is proved. We provide a
numerical example to verify the existence of transcritical bifurcation for the
system (6). We check the transcritical bifurcation for
$d_{2}=(d_{2})_{TC}=0.052246635114$, which corresponds to the transcritical
bifurcation threshold $\beta=\beta_{TC}$. All other parameters are taken from
set (1) of Table 1. The predator extinct equilibrium point, $E_{1}$, exchanges
its stability with interior equilibrium points, $E^{*}_{1}$, through
transcritical bifurcation as depicted in Fig. 3b. ∎

######  Remark 10.

Two transcritical bifurcations for the system (6) occur at
$r=(r)_{{TC}_{1}}=0.226735001383921$ and $r=(r)_{{TC}_{2}}=0.88701499861607$
(refer to Fig. 4), which correspond to the transcritical bifurcation threshold
$\beta=\beta_{TC}$. All other parameters from set (2) of Table (1) .

###  4.3 Hopf bifurcation

######  Theorem 11.

At the bifurcation threshold $k=k_{H}$, the model system (6) endures a hopf
bifurcation at the coexisting equilibrium point $E^{*}$ if the following
conditions are met: (i) $Tr\left(J_{E^{*}}\right)|_{k=k_{H}}=0$, (ii)
$Det\left(J_{E^{*}}\right)|_{k=k_{H}}>0$, (iii)
$\dfrac{d}{dk}\Big{[}Re\left(\lambda_{i}(k)\right)\Big{]}_{k=k_{H}}\neq 0$ for
$i=1,2$; where $\lambda_{i}$ is an eigenvalue of the Jacobian matrix
corresponding to $E^{*}$.

###### Proof.

For the proof, refer to Appendix A.3. ∎

It is difficult to establish the analytic expression for the Hopf-bifurcation
threshold because of the algebraic complexity of our model system. However, we
numerically test the Hopf-bifurcation for a set of parameter values. If we
take all of the parameter values from set (1) of Table (1) and set
$\alpha=0.25$ and $c=3$, the interior equilibrium point of the system (6)
$E_{2}^{*}=(0.101007273887931,1.809506754599451)$ becomes stable through Hopf-
bifurcation at threshold value $k=k_{H}=0.509642226$ (see Fig. 10).

### Nature of Hopf-bifurcation

By determining the first Lyapunov number at the equilibrium point $E_{H}^{*}$,
we can investigate the direction of the Hopf-bifurcation or the stability of
limit cycle. To begin, we use the transformations $x=\hat{X}-x_{H}^{*}$ and
$y=\hat{Y}-y_{H}^{*}$ to move $E_{H}^{*}(x_{H}^{*},y_{H}^{*})$ to the origin.
Then the model system (6) in the neighbourhood of the origin can be written
as:

$\displaystyle\frac{d\hat{X}}{dt}$
$\displaystyle=\alpha_{01}\hat{Y}+\alpha_{10}\hat{X}+\alpha_{02}\hat{Y}^{2}+\alpha_{20}\hat{X}^{2}+\alpha_{11}\hat{X}\hat{Y}+\alpha_{30}\hat{X}^{3}+\alpha_{03}\hat{Y}^{3}+\alpha_{21}\hat{X}^{2}\hat{Y}+\alpha_{12}\hat{X}\hat{Y}^{2}+F_{1}(\hat{X},\hat{Y})$
(9) $\displaystyle\frac{d\hat{Y}}{dt}$
$\displaystyle=\beta_{01}\hat{Y}+\beta_{10}\hat{X}+\beta_{02}\hat{Y}^{2}+\beta_{20}\hat{X}^{2}+\beta_{11}\hat{X}\hat{Y}+\beta_{30}\hat{X}^{3}+\beta_{03}\hat{Y}^{3}+\beta_{21}\hat{X}^{2}\hat{Y}+\beta_{12}\hat{X}\hat{Y}^{2}+F_{2}(\hat{X},\hat{Y}),$

where, the $\alpha_{ij}=\frac{1}{i!j!}\frac{\partial^{i+j}F}{\partial
x^{i}\partial y^{j}}|_{(E_{H}^{*},\theta_{H})}$ ,
$\beta_{ij}=\frac{1}{i!j!}\frac{\partial^{i+j}G}{\partial x^{i}\partial
y^{j}}|_{(E_{H}^{*},\theta_{H})}$ (i, j = 0, 1, 2 and 3) and
$F_{i}(\hat{x},\hat{y})$ ($i=1,2$) are power series in powers of
$\hat{X}^{i}\hat{Y}^{j}$ satisfying $i+j\geq 4$. The expressions of
$\alpha_{ij}$ and $\beta_{ij}$ have been provided in Appendix B.3.

We know that expression for Lyapunov coefficient ($\sigma$) [54] for a general
planar system (9) is given by:

$\displaystyle\sigma=$
$\displaystyle-\frac{3\pi}{2\alpha_{01}\Delta^{\frac{3}{2}}}\\{[\beta_{10}\alpha_{10}(\beta_{02}\alpha_{11}+\alpha_{11}^{2}+\beta_{11}\alpha_{02})+\alpha_{01}\alpha_{10}(\beta_{11}\alpha_{20}+\beta_{11}^{2}+\beta_{02}\alpha_{11})+\beta_{10}^{2}(2\alpha_{02}\beta_{02}+\alpha_{11}\alpha_{02})$
$\displaystyle-2\beta_{10}\alpha_{10}(-\alpha_{02}\alpha_{20}+\beta_{02}^{2})-2\alpha_{01}\alpha_{10}(-\beta_{02}\beta_{20}+\alpha_{20}^{2})-\alpha_{01}^{2}(\beta_{20}\beta_{11}+2\beta_{20}\alpha_{20})+(-2\alpha_{10}^{2}+\alpha_{01}\beta_{10})$
$\displaystyle(-\alpha_{20}\alpha_{11}+\beta_{02}\beta_{11})]-(\beta_{10}\alpha_{01}+\alpha_{10}^{2})[3(-\alpha_{30}\alpha_{01}+\beta_{03}\beta_{10})+2\alpha_{10}(\beta_{12}+\alpha_{21})+(-\beta_{21}\alpha_{01}+\alpha_{12}\beta_{10})]\\},$

where, $\Delta=\alpha_{10}\beta_{01}-\alpha_{01}\beta_{10}.$ Using the model
parameter values provided above in this subsection, we calculate
$\alpha_{ij},\beta_{ij}$, and get $\sigma=120.438,$ which implies that the
Hopf bifurcation is of subcritical type.

##  5 Effect of fear level, and Birth Rate of Prey

Effect of fear level on prey and predator equilibrium densities: Here, we are
interested in observing the changes of the prey population with respect to the
fear level. Since, the prey population at coexisting equilibrium is given by
the positive solution of the Eq. (7). Therefore, we differentiate the Eq. (7)
with respect to parameter $k$ (level of fear) and we obtain

$\dfrac{dx}{dk}\bigg{|}_{x=x^{*}}=\frac{-\left(x^{7}\dfrac{dA_{1}}{dk}+x^{6}\dfrac{dA_{2}}{dk}+x^{5}\dfrac{dA_{3}}{dk}+x^{4}\dfrac{dA_{4}}{dk}+x^{3}\dfrac{dA_{5}}{dk}+x^{2}\dfrac{dA_{6}}{dk}+x\dfrac{dA_{7}}{dk}+\dfrac{dA_{8}}{dk}\right)}{(7A_{1}x^{6}+6A_{2}x^{5}+5A_{3}x^{4}+4A_{4}x^{3}+3A_{5}x^{2}+2A_{6}x+A_{7})}\bigg{|}_{x=x^{*}},$

where $A_{1},A_{2},A_{3},A_{4},A_{5},A_{6},A_{7},$ and $A_{8}$ have same
meaning as mentioned in Eq. (7) in subsection (3.4). Now, we have

$\displaystyle\dfrac{dx}{dk}\bigg{|}_{x=x^{*}}=\frac{\gamma}{\delta}\bigg{|}_{x=x^{*}},$

where

$\displaystyle\gamma=$ $\displaystyle-\alpha(da+x(-1+m)(\beta c-b\alpha
d+x(-1+m)d))\left((-1+m)\beta(da+x(-1+m)(-\alpha bd+\beta c+x(-1+m)d))\right.$
$\displaystyle\left.+e\left(d_{1}+xd_{2}\right)(a+x(-1+m)(-\alpha
b+x(-1+m)))^{2}\right),$ $\displaystyle\delta=$
$\displaystyle(-1+m)\left(-c(-1+m)\beta^{2}(-2k\alpha d+e)(a+x(-2\alpha
b+3x(-1+m))(-1+m))+2d(-1+m)\beta(-\alpha b+2x\right.$
$\displaystyle\left.(-1+m))(-e+k\alpha d)(a+x(-1+m)(-\alpha
b+x(-1+m)))+3re^{2}(-\alpha b+2x(-1+m))(a+(-1+m)x\right.$
$\displaystyle\left.(-\alpha b+x(-1+m)))^{2}+2\alpha
c^{2}(-1+m)^{2}\beta^{3}kx\right)-e(a+x(-1+m)(-\alpha
b+(-1+m)x))\left(d_{2}\left(a^{2}(-k\alpha d+e)\right.\right.$
$\displaystyle\left.\left.+(-1+m)xa((8x(-1+m)-5\alpha b)(-k\alpha
d+e)-2k\alpha\beta c)+x^{2}(-1+m)^{2}(2\alpha\beta ck(2\alpha
b-3(-1+m)x)+\right.\right.$ $\displaystyle\left.\left.(\alpha
b-mx+x)(-7x(-1+m)+4\alpha b)(-k\alpha d+e))\right)+(-1+m)d_{1}(3a(-\alpha
b+2x(-1+m))(-\alpha dk+e)-a\alpha\right.$ $\displaystyle\left.\beta
ck+x(-1+m)(\alpha\beta ck(-5x(-1+m)+3\alpha b)-3(-\alpha b+x(-1+m))(-\alpha
b+2x(-1+m))(-e+k\alpha d)))\right).$

Further, let

$\displaystyle(da+x(-1+m)(-\alpha bd+\beta c+d(-1+m)x))<$ $\displaystyle
0,~{}\delta<0.$ (10)

The Eq. (10) ensures that $\gamma>0$. If Eq.(10) holds, then
$\dfrac{dx}{dk}\bigg{|}_{x=x^{*}}<0$. It indicates that as fear in a prey
population against predation increases, prey equilibrium density decreases
monotonically; hence, an increasing level of fear gradually decreases the prey
population density. Because the increasing level of fear compels prey to
change their foraging tactics and spend more time on vigilance, which
ultimately hampers their time for reproduction. As a result, prey birth rates
decrease as fear levels rise, which eventually decreases the level of prey
population density. Additionally, when $x^{*}$ is feasible, the predator
population at the coexisting equilibrium point is given by
$y^{*}=\frac{1}{e}\left(-d+\frac{c\beta x^{*}(1-m)}{a+b\alpha
x^{*}(1-m)+(1-m)^{2}(x^{*})^{2}}\right),$ whenever it is feasible. Now, we
differentiate $y$ with respect to $k$ as follows to see how $k$ affects the
predator population:

$\dfrac{dy^{*}}{dk}=\frac{c\beta(-m+1)}{e}\left[\frac{a+x^{*}(-m+1)b\alpha+(-m+1)^{2}(x^{*})^{2}-\left(b(-m+1)\alpha+2x^{*}(-m+1)^{2}\right)x^{*}}{\left(a+x^{*}(-m+1)b\alpha+(x^{*})^{2}(-m+1)^{2}\right)^{2}}\right]\dfrac{dx^{*}}{dk},$
$\dfrac{dy^{*}}{dk}=\frac{\beta
c(-m+1)}{e}\left[\frac{\left(a-(x^{*})^{2}(-m+1)^{2}\right)}{\left(a+x^{*}(-m+1)b\alpha+(x^{*})^{2}(-m+1)^{2}\right)^{2}}\right]\dfrac{dx^{*}}{dk}.$

Further, let $\left(a-(x^{*})^{2}(-m+1)^{2}\right)>0$. It could be observed
that the change in the equilibrium density of the predator population depends
on changes in the prey equilibrium density, i.e., the prey equilibrium density
decreases, then the predator equilibrium density also decreases with respect
to level of fear $k$. Hence, the increased level of fear reduces the prey
population directly and the predator population indirectly. In a biological
sense, increased fear not only forces the prey populations to forage less but
also to tighten their group defences, which consequently hampers the food
intake rate of predators. As a result, predator numbers decrease as the fear
level increases.
Effect of prey birth rate on prey equilibrium density: In a similar fashion,
as we have discussed in the above subsection, we can observe the impact the
prey’s birth rate $(r)$ on the prey equilibrium density. The change in prey
population with respect to $r$ is given by

$\dfrac{dx}{dr}\bigg{|}_{x=x^{*}}=\frac{N}{\delta}\bigg{|}_{x=x^{*}},$

where, $N=-e^{2}(a+x(-1+m)(-\alpha b+(-1+m)x))^{3}<0.$ Therefore, if Eq. (10)
holds then we have $\dfrac{dx}{dr}\bigg{|}_{x=x^{*}}>0$. It indicates that as
prey birth rate increases, prey equilibrium density increases monotonically.
It may happens for two reasons: (1) as prey birth rates rise, so do prey
numbers; and (2) as prey numbers rise, so do prey group defenses, which reduce
the chances of prey being predated and indirectly contribute to prey numbers.

##  6 Spatial Model System

Reaction-diffusion equations are widely used as models for spatial effects in
ecology. They support three important types of ecological phenomena: the
existence of a minimal patch size necessary to sustain a population and the
formation of spatial patterns in the distributions of populations in
homogeneous environments [55]. The spatial model between the interaction of
prey $x(u,v,t)$ and predator $y(u,v,t)$ is described by the following set of
partial differential equations:

$\displaystyle\frac{{\partial x}}{{\partial t}}=$
$\displaystyle\frac{{rx}}{{(1+k\alpha
y)}}-{d_{1}}x-{d_{2}}{x^{2}}-\frac{{\left({-m+1}\right)\beta
xy}}{{a+(-m+1)\alpha bx+{{(-m+1)}^{2}}{x^{2}}}}+{D_{1}}{\nabla^{2}}x,$ (11)
$\displaystyle\frac{{\partial y}}{{\partial t}}=$ $\displaystyle-
dy-e{y^{2}}+\frac{{\left({-m+1}\right)\beta cxy}}{{a+(-m+1)\alpha
bx+{{(-m+1)}^{2}}{x^{2}}}}+{D_{2}}{\nabla^{2}}y,$

with initial conditions

$x(u,\,v,\,0)>0,\,\,y(u,\,v,\,0)>0\,\,\,\,\,{\rm for}\,\,\,(u,\,v)\in\Omega,$

and boundary conditions

$\frac{{\partial x}}{{\partial\vartheta}}=\frac{{\partial
y}}{{\partial\vartheta}}=0,\,\,\,\,\,\,(u,\,v)\in\partial\Omega\,\,\,\,{\rm
for\,\,\,all}\,\,\,\,t,$

where $D_{1}$ and $D_{2}$ are diffusion coefficients of prey and predator
population respectively, the laplacian operator
${\nabla^{2}}\equiv\frac{{{\partial^{2}}}}{{\partial{u^{2}}}}+\frac{{{\partial^{2}}}}{{\partial{v^{2}}}}$
and $\vartheta$ is outward normal to the $\partial\Omega$ .

###  6.1 Dynamics of spatial Model system

To understand the spatial dynamics of model system (11), we consider the
linearized form of the system about $E^{*}(x^{*},y^{*})$ as follows:

$\displaystyle\frac{{\partial\bar{x}}}{{\partial t}}=$
$\displaystyle{j_{11}}\bar{x}+{j_{12}}\bar{y}+{D_{1}}{\nabla^{2}}\bar{x},$
(12) $\displaystyle\frac{{\partial\bar{y}}}{{\partial t}}=$
$\displaystyle{j_{21}}\bar{x}+{j_{22}}\bar{y}+{D_{2}}{\nabla^{2}}\bar{y},$

where we introduce small perturbations $\bar{x}=x-{x^{*}}$ and
$\bar{y}=y-{y^{*}}.$ Let us suppose that the solution of model system (11) is
of the form:

$\left(\begin{array}[]{l}{\bar{x}}\\\ {\bar{y}}\\\
\end{array}\right)=\left(\begin{array}[]{l}{l_{1}}\\\ {l_{2}}\\\
\end{array}\right)\exp({\lambda_{\Lambda}}t)\cos({\Lambda_{u}}u)\cos({\Lambda_{v}}v),$

where $l_{1}$, $l_{2}$ are sufficiently small constants, $\Lambda_{u}$ and
$\Lambda_{v}$ are the component of wave number along $u$ and $v$ direction
respectively and $\Lambda_{k}$ is wave length.
The jacobian matrix of linearized model system (12) is

$\bar{J}={\left({\begin{array}[]{*{20}{c}}{{j_{11}}-{D_{1}}{\Lambda^{2}}}&{{j_{12}}}\\\
{{j_{21}}}&{{j_{22}}-{D_{2}}{\Lambda^{2}}}\\\ \end{array}}\right)},$

where $\Lambda$ is wave number given by
$\Lambda^{2}=\Lambda^{2}_{u}+\Lambda^{2}_{v}.$ The characteristic equation of
$\bar{J}$ is given by

$\lambda_{\Lambda}^{2}+{\rho_{1}}{\lambda_{\Lambda}}+{\rho_{2}}=0,$

where,

$\displaystyle{\rho_{1}}$ $\displaystyle=$
$\displaystyle-({j_{11}}+{j_{22}})+\left({{D_{1}}+{D_{2}}}\right){\Lambda^{2}},$
$\displaystyle{\rho_{2}}$ $\displaystyle=$
$\displaystyle{D_{1}}{D_{2}}{\Lambda^{4}}-\left({{j_{11}}{D_{2}}+{j_{22}}{D_{1}}}\right){\Lambda^{2}}+{j_{11}}{j_{22}}-{j_{12}}{j_{21}}.$

The equilibrium point ${E^{*}}\left({{x^{*}},\,\,{y^{*}}}\right)$ is locally
asymptotically stable in the presence of diffusion if $\rho_{1}>0$ and
$\rho_{2}>0$.

###  6.2 Turing Instability

Turing instability occurs when steady state is stable in the absence of
diffusion $(D_{1}=D_{2}=0)$ and unstable in the presence of diffusion
${D_{i}}\neq 0\,\,\,\,\rm{for}\,\,\,{\it i}=1,\,\,2.$ The condition for
occurrence of Turing instability is given by

$H({\Lambda^{2}})={D_{1}}{D_{2}}{\Lambda^{4}}-({j_{11}}{D_{2}}+{j_{22}}{D_{1}})\Lambda^{2}+{j_{11}}{j_{22}}-{j_{12}}{j_{21}}<0.$

The function $H(\Lambda^{2})$ is a quadratic equation in $\Lambda^{2}$. It
represents a parabola opening upward and $H(\Lambda^{2})$ has its minimum
$({H_{\min}})$ at some value $\Lambda_{c}^{2}$ of $\Lambda^{2}$ at vertex of
the parabola where
$\Lambda_{c}^{2}=\frac{{{j_{11}}{D_{2}}+{j_{22}}{D_{1}}}}{{2{D_{1}}{D_{2}}}}.$
Then the condition that $H(\Lambda_{c}^{2})<0$ gives,

${\left({{j_{11}}{D_{2}}+{j_{22}}{D_{1}}}\right)^{2}}>4{D_{1}}{D_{2}}\left({{j_{11}}{j_{22}}-{j_{12}}{j_{21}}}\right).$
(13)

Therefore, $H(\Lambda^{2}_{c})$ will be negative, when condition of Eq. (13)
satisfies. We summarize the above discussion in the form of following theorem.

######  Theorem 12.

At the homogeneous steady state $E^{*}(x^{*},y^{*})$, the model system
undergoes Turing instability if (i) ${D_{1}}{j_{22}}+{D_{2}}{j_{11}}>0,$ (ii)
${\left({{D_{1}}{j_{22}}+{D_{2}}{j_{11}}}\right)^{2}}>4{D_{1}}{D_{2}}\left({{j_{11}}{j_{22}}-{j_{12}}{j_{21}}}\right),$
(iii)
${j_{11}}+{j_{22}}<0,\,\,\,\,\,\,\,\,\,{j_{11}}{j_{22}}-{j_{12}}{j_{21}}>0.$

The diffusion induced instability in the model system (11) has been discussed
by taking the following parameter values:

$\displaystyle r$ $\displaystyle=$ $\displaystyle
1.1,\,\,d_{1}=0.05,\,\,d_{2}=0.055,\,\,\beta=0.071,\,\,m=0.3,\,\,k=0.5,\,\,\alpha=1,\,\,d=0.05,\,\,e=0.05,\,\,c=5,$
$\displaystyle a$ $\displaystyle=$ $\displaystyle
0.1,\,\,b=0.125,D_{1}=0.005,\,\,\,D_{2}=2.$ (14)

At the set of parametric values given in Eq. (14), we obtain
${E^{*}}\left({{x^{*}},\,\,{y^{*}}}\right)=(0.0488345,1.30182)$. Also

$j_{11}+j_{22}=-0.0293104<0\,\,\,\,\,{\rm
and}\,\,\,\,\,j_{11}j_{22}-j_{12}j_{21}=0.0922047>0$

ensure the locally asymptotical stability of the system without diffusion. For
diffusive model (11), $H(\Lambda^{2})<0$ for $1.7001<\Lambda^{2}<5.42347.$ For
the above set of parameter values, the plot of Turing instability is given in
Fig 2.

(a)

(b)

Figure 2: (a) Behavior of $H(\Lambda^{2})$ for the occurrence of Turing
instability and (b) for different values of level of fear factor $k$, the real
part of the dispersion relation of model system (11). For the values of other
parameters, refer to Eq. (14).

##  7 Numerical simulation

Now, we perform the numerical simulation to validate our obtained analytical
results of spatial and non spatial prey-predator systems.

###  7.1 Numerical Simulation of Non-Spatial Model

Influence of Intra-Prey Competition Parameter: Fig. 3a and Fig. 3b show that
how the existence and stability of various feasible equilibrium points are
influenced by the intra-prey competition parameter $d_{2}$. The interior
equilibrium point $E^{*}_{3}$ exist always and is stable. It can be visualized
in Fig. 3b as the corresponding $Max(Re(\lambda_{i}))$ is negative throughout.
The interior equilibrium point $E^{*}_{1}$ emerges at $(d_{2})_{TC}$, whereas
the predator extinct equilibrium point $E_{1}$ becomes unstable due to
stability exchange with $E^{*}_{1}$ via transcritical bifurcation. The
interior equilibrium point $E^{*}_{1}$ exist and remain stable in the region
$R_{2}$ (presented in Fig. 3a and 3b) because the corresponding
$Max(Re(\lambda_{i}))$ in Fig.3b is negative throughout. The interior
equilibrium point $E^{*}_{2}$ exist in the region $R_{1}\cup R_{2}$ and is
unstable as the corresponding $Max(Re(\lambda_{i}))$ is positive throughout
regions (please refer the Fig. 3b). As $d_{2}$ increases, interior equilibrium
points $E^{*}_{1}$ and $E^{*}_{2}$ in region $R_{2}$ move closer to each
other, and finally collide at $d_{2}=(d_{2})_{SN}$ via saddle-node
bifurcation. The collided equilibrium point $E^{*}_{12}=(4.7809,3.0553)$ will
be half stable. The regions $R_{1}$ and $R_{2}$ are regions of bistabilty as
there are two local basins of attraction. Local attractors in region $R_{1}$
are the interior equilibrium point $E^{*}_{3}$ and the predator extinct
equilibrium point $E_{1}$. Therefore, depending on the initial population
size, the predator population may extinct or the system would get stabilised
at a low prey density. The interior equilibrium points $E^{*}_{1}$ and
$E^{*}_{3}$ are local attractors in region $R_{2}$, so the population always
survives but only changes the level of population. Because there is only one
basin of attraction, which is around the interior equilibrium point of
$E^{*}_{3}$, the population survives at a low prey density in region $R_{3}$.
Hence, from Fig. 3b, we conclude that there is always a possibility of the
coexistence of prey and predator species. The scenario of predator extinction
at low values of $d_{2}$ arises because of a relatively higher level of prey
population than of predator population. It leads to robust group defence by
prey against predation, which eventually forces predators to starve to death.
As a result, predators are no longer present in the system. For a large enough
increment in values of $d_{2}$, two possibilities occur for stable coexistence
of populations depending on the initial populations: (a) both the population
coexist with higher levels of prey population and (b) both the population
coexists with lower levels of prey population. Whichever of these two
possibilities is the case, both species somehow find a way to live with each
other and do not affect the biodiversity of the ecosystem. In the case of the
first possibility, the increasing value of $d_{2}$ has two effects: (1) it
decreases prey numbers directly (2) as well as indirectly by causing possible
laxity in the group defence of prey, which allows predators to hunt more prey.
Therefore, the gradual increment in $d_{2}$ removes the possibility (a),
leaving the system with only the possibility (b). The possibility (b) occurs
for all values of $d_{2}$ in which both the populations live together with
relatively lower densities without forcing each other toward extinction.

(a)

(b)

Figure 3: All parameter values, apart from $d_{2}$, are chosen from set (1) in
Table 1. (a) One parameter bifurcation diagram of model system (6),
corresponding to $d_{2}$. The green, black, and red curves represent the prey
density at equilibrium points $E^{*}_{1},E^{*}_{2}$, and $E^{*}_{3},$
respectively. The cyan and violet vertical lines represent the transcritical
and saddle-node bifurcation thresholds $d_{2}=(d_{2})_{TC}=0.0522466351144$
and $d_{2}=(d_{2})_{SN}=0.077854268968873742$, respectively. (b) The green,
black, red, and magenta curves represent the maximum of the real part of the
eigenvalues corresponding to the interior equilibrium points
$E^{*}_{1},E^{*}_{2},E^{*}_{3}$, and $E_{1}$, respectively

.

Influence of Birth Rate of Prey Parameter: Here, we focus on how the birth
rate of prey influences the prey-predator dynamics. Here, Fig. 4a and Fig. 4b
depict the dependence of existence and stability of various feasible
equilibrium points on the prey birth rate parameter $r$. Only the population
extinct equilibrium point $E_{0}$ is feasible and a global attractor when
$r<r_{TC}$. Therefore, both populations are eliminated from the system when
$r<r_{TC}$. As $r$ achieves $r_{TC}$, the population extinct equilibrium point
$E_{0}$ exchanges its stability to become unstable permanently via
transcritical bifurcation, and the predator extinct equilibrium point $E_{1}$
comes into role. In the region $R_{1}$, the $E_{1}$ remains stable. When $r$
reaches $r_{{TC}_{1}}$, the predator extinct equilibrium point $E_{1}$ becomes
unstable by exchanging stability with the interior equilibrium point
$E^{*}_{2}$ via transcritical bifurcation. The predator extinct equilibrium
point $E_{1}$ remains unstable in the region $R_{2}.$ At $r=r_{{TC}_{2}}$, it
exchanges stability with the interior equilibrium point $E^{*}_{1}$ via
transcritical bifurcation and returns to stability. Region $R_{2}\cup R_{3}$
is the region of existence for interior equilibrium point $E^{*}_{2}$ and it
remains stable because the corresponding $Max(Re(\lambda_{i}))$ is remain
negative throughout this region (presented in Fig. 4b). On the other hand,
interior equilibrium point $E^{*}_{1}$ will unstable whenever it exists.
Interior equilibrium point $E^{*}_{1}$ exist only in the region $R_{3}$ (refer
to Fig. 4a), where it is always unstable because the corresponding
$Max(Re(\lambda_{i}))$ in Fig. 4b is positive throughout. On other hand, it
can be observed that interior equilibrium $E^{*}_{1}$ remains unstable
whenever exists. As $r$ moves through region $R_{3}$, the equilibrium points
$E^{*}_{1}$ and $E^{*}_{2}$ move towards each other, eventually collide at
$E^{*}_{12}=(1.247196,2.438599)$ when $r=r_{SN}$. In region $R_{4}$, only the
predator extinct equilibrium point remains stable and no interior equilibrium
points exist. Since, $R_{3}$ is the bistability region for $E_{1}$ and
$E^{*}_{2},$ therefore, it is the region of local basins for $E_{1}$ and
$E^{*}_{2}$. So the predator population may extinct or the system could
stabilise around co-existence equilibrium point depending on the initial
population size. In the region $R_{2}$, the both the population survive
because the interior equilibrium point $E^{*}_{2}$ is the only basin of
attraction for system. $R_{1}$ and $R_{4}$ are the regions in which predator
population would go to extinction, as the predator extinct equilibrium point
$E_{1}$ is the only local attractor for the system.

As a result, predator species do not exist in the system for both extremely
high and extremely low birth rates of prey. When birth rates of prey are
extremely low, predators do not get enough food to survive and vanish from the
system. In the case of extremely high prey birth rates, prey enable themselves
to tighten the group defence as much as they can so that predators are not
able to hunt them down. It triggers the predators’ starvation and eventually
their extinction from the system. When both prey and predators exist in the
system and prey birth rates increases gradually, then a case of bistability
arises between the predator extinct equilibrium and a coexisting equilibrium.
If the initial predator species’ size is not large enough, then high levels of
prey population will cause the extinction of predators by using their group
defence. In the case of predators with a large enough initial size, both
species live together until prey birth rates exceed a threshold, because after
this, prey numbers will be high enough to make predators starve by using their
group defence and cause predator extinction.

(a)

(b)

Figure 4: All parameter values, apart from $r$, are chosen from set (2) in
Table 1. (a) Diagram of a single parameter bifurcation for model system (6),
taking $r$ as a bifurcation parameter. The blue, and black curves represent
the existence of interior equilibrium points $E^{*}_{1},E^{*}_{2}$
respectively. The red, the violet, the cyan, and the green vertical lines
represent the transcritical bifurcation at $r=(r)_{TC}=0.2$, transcritical
bifurcation at $r=(r)_{{TC}_{1}}=0.226735001383921$, transcritical bifurcation
at $r=(r)_{{TC}_{2}}=0.887014998616079$, and saddle-node bifurcation at
$r=(r)_{SN}=3.273118754537237107271607$, respectively. (b) The blue, black,
and magenta curves represent the maximum of real part of eigen values
corresponding to interior equilibrium points $E^{*}_{1},E^{*}_{2}$ and
predator free equilibrium point $E_{1}$, respectively.

Stability analysis of equilibria: Only the population extinct equilibrium
point is feasible and globally stable if $r<d_{1}$, as stated in Theorem 4(i).
Therefore whenever this condition holds, both populations go to extinction as
time goes on, as depicted in the Fig. 5a. If we choose a suitable set of
parametric values which satisfy the local stability condition of the predator
extinct equilibrium point $E_{1}$ as stated in Theorem 4(ii), then the
predator population will become extinct and the prey population will approach
$\left(\frac{r-d_{1}}{d_{2}}\right)$, as shown in Fig. 5b.

(a) The local stability of $E_{0}(0,0)$ with $d_{1}=1.2$

(b) The local stability of the predator extinct equilibrium $E_{1}$ with
$d_{1}=1,d_{2}=2,c=1$

Figure 5: The figure shows the local stability of equilibria when the
numerical values of parameters are taken same as set (1) in Table 1 apart from
stated above.

Further, we vary $d_{2}$ by taking all other parametric values from the set
(1) of the Table 1, the interior equilibrium point (coexisting equilibrium)
$E^{*}_{3}$ always exists and is locally asymptotically stable. If we consider
$d_{2}=-.055$ and other parametric values same as set (1) of the Table 1 then
interior equilibrium point $E^{*}_{23}=(0.0234354,1.2758912)$ exists. In Fig.
6a, we see that both populations initially exhibit damped oscillations, but
over time, they move closer to a positive level. Therefore, Fig. 6a and 6b
represent the local stability of interior equilibrium point. As a result, we
can achieve the conclusion that both populations may survive throughout time
and the ecosystem may enter a stable state when both populations coexist under
specific environmental conditions.

(a) Long term behavior of the system (6).

(b) Phase portrait corresponding to time series plot Fig. 6a.

Figure 6: The figure shows the local stability of coexisting equilibrium
$E^{*}_{3}$. Here, all parametric values, apart from $d_{2}=0.055$, are chosen
from set (1) of Table 1.

Now, we look at how the fear level and the prey birth rate affect the
equilibrium density of prey and predator. Whenever the conditions (10) hold,
the equilibrium density of prey decreases as the level of fear increases,
which is depicted in Fig. 7a. As shown in Fig. 7b, whenever condition
($\left(a-(x^{*})^{2}(-m+1)^{2}\right)>0$) holds and the prey equilibrium
density decreases as the level of fear increases, the predator equilibrium
density decreases as well. It happens because prey modifies its foraging
behaviour and spends more time being vigilant as its level of fear rises,
which leads to decreased food intake and reduced reproduction. The population
density of prey at equilibrium drops as a result. This altered prey foraging
behaviour reduces the amount of food available to predators, hinders their
capacity to consume it, and lowers the predator population density at
equilibrium. Whenever $\delta<0$ (refer to Eq. (10)), the equilibrium density
of prey and predator increases as the birth rate of prey increases, which is
depicted in Fig. 8. Increasing the prey birth rate does not only increase prey
numbers but also the food availability for predators. Hence, both prey and
predator populations increase at equilibrium.

(a) Effect of fear on prey density at coexisting equilibrium.

(b) Effect of fear on predator density at coexisting equilibrium.

Figure 7: The figure depicts the effect of fear in prey on equilibrium state
of species. Here, all parametric values are taken from set (1) in Table 1.

(a) Effect of prey birth rate on prey density at coexisting equilibrium. (b)
Effect of prey birth rate on predator density at coexisting equilibrium.

Figure 8: The figure depicts the effect of prey birth rate on equilibrium
state of species. Here, all parametric values are taken from set (1) in Table
1.

(a) Region plot for extinction criteria of predator species

(b) Long term dynamics of the model system (6).

Figure 9: The figure ensures the extinction of predator species from prey-
predator system whenever values $r$, and $\beta$ have been chosen from blue
shaded region of Fig. 9a, i.e. $(r=0.6,,\beta=0.00005)$. Here, all parameter
values, apart from $r$ and $\beta$, are chosen from set (1) of Table 1.

The extinction of any species may arise the imbalance in the ecosystem.
Therefore, its study is an important aspect of prey-predator interaction. In
Fig. 9a, the blue-shaded region is the region of extinction for predator
species in $(r-\beta)$ space. This blue-shaded region corresponds to the
parametric condition for the extinction of predators stated in Theorem 2(ii).
The green-shaded region of $(r-\beta)$ space in the Fig. 9a is the region of
uncertainty for predator species extinction, because predator species may or
may not go extinct in this region. In Fig. 9b, we observe that the solutions
curves of the model system (6) converge to predator extinct equilibrium. For
this time series plot, we have chosen the numerical values of $r$, and $\beta$
from the blue shaded region which ensure the extinction of predator species.
Therefore, from Fig. 9, we conclude that the predator species goes to
extinction for the small values of $r$ and $\beta.$

In Fig. 10, we observe the occurrence of periodic solution of the system (6)
with a gradual increment in the level of fear (related parameter is $k$). From
Fig. 10a, we notice that at a low level of fear, predator species extinct from
the system. As the level of fear increases, periodic solutions occur, but the
system remains unstable. But the higher level of fear eliminates the periodic
solutions and stabilises the system. Further, we have that the interior
equilibrium point $E^{*}_{2}$ is unstable for k=0.509554, so the solution
spirals out, as shown in the Fig. 10d. This is also demonstrated by
corresponding time series plot as shown in Fig. 10b and 10c (refer to blue
curves). The gradual increment in the fear enforces the occurrence of periodic
solutions for the system (6). When fear reaches the Hopf bifurcation threshold
$k=k_{H}=0.509642226$, an unstable limit cycle occurs around $E^{*}_{2}$ via
hope bifurcation, as shown in Fig. 10e. The presence of periodic oscillations
in the time series corresponding to phase diagram Fig. 10e has also been shown
in Fig. 10b and 10c by magenta curves. With another gradual increase in the
level of fear, the unstable limit cycle surrounding the stable interior
equilibrium point $E^{*}$ disappears and the solution spirals into $E^{*}_{2}$
for the appropriate initial conditions as shown in Fig. 10f and red curves of
Fig. 10b and 10c.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 10: All parameter values, apart from $\alpha$ = 0.25, c = 3 and k, are
chosen from set (1) in Table 1. Here, the initial condition is
$(x_{0},y_{0})=(0.0999994016,1.784623932176383)$. (a) The blue and red colour
curves represent the time series of $y_{max}$ and $y_{min}$, respectively. The
blue, magenta, and red curves in Figs.(10b)-(10c) represent the time series
corresponding to $E^{*}_{2}$ before the Hopf bifurcation threshold, at the
Hopf bifurcation threshold, and after the Hopf bifurcation threshold,
respectively. (d) It represents the instability of $E^{*}_{2}$ before the
occurrence of the limit cycle through Hopf bifurcation. (e) It represents the
existence of a limit cycle. (f) It represents the stability of $E^{*}_{2}$
after the Hopf bifurcation.

The stability region for the predator extinct equilibrium point $E_{1}$ is
represented by the blue shaded region of Fig. 11a in $(\alpha,\beta)$ space.
If $(\alpha,\beta)$ is taken from the blue shaded region, it meets the locally
asymptotically stable condition of the predator extinct equilibrium point
$E_{1}$ stated in Theorem 4(ii), and the predator species will eventually
become extinct in the system. The cyan-shaded region in Fig. 11a is the region
of instability for the predator extinct equilibrium point $E_{1}$. A time
series of a numerical example for the stability of a predator free equilibrium
point is provided in Fig. 11b by choosing the values of $(\alpha,\beta)$ from
the blue shaded region. In Fig. 11b, we notice that the solution curves
approach to predator extinct equilibrium $E_{1}$ as time passes. Therefore, it
ensure the local stability of $E_{1}$ equilibrium.

(a) Region plot for extinction criteria of predator species in
$(\alpha-\beta)$ space.

(b) Long term dynamics of the model system (6).

Figure 11: The figure ensures the local stability of predator extinct
equilibrium $E_{1}$ of system (6) whenever values $\alpha$, and $\beta$ have
been chosen from blue shaded region of Fig. 11a, i.e.
$(\alpha=5,,\beta=0.025)$. Here, all parameter values, apart from $\alpha$ and
$\beta$, are chosen from set (1) of Table 1.

As shown in Fig. 12a, the blue-shaded region of $(\alpha-\beta)$ space is the
region of persistence for prey species. If $(\alpha,\beta)$ is taken from the
blue shaded region, it meets the persistence condition of the prey species
stated in Theorem 3(i) and ensures the survival of prey species in the system.
The red-shaded region in Fig. 12a represents the region of uncertainty for
prey species persistence, because prey species may persist or become extinct
in this region. A time series of a numerical example for the persistence of
prey species is provided in Fig. 12b by selecting values of $(\alpha,\beta)$
from the blue shaded region.

(a)

(b)

Figure 12: All parameter values, apart from $r=5$, $\alpha$ and $\beta$, are
chosen from set (1) in Table 1. (a) The blue-shaded region is the region of
persistence for prey species. (b) The time series for a numerical example with
$(\alpha,\beta)$=(0.5, 0.002) and initial condition (x(0), y(0))=(80, 60).

###  7.2 Numerical Simulation of Spatial Model

In this subsection, we perform the numerical simulation of spatial model. In
Fig. 13, we observe the effect of time evolution on the spatial density
distribution of predator population. For time iteration $2000$, the density
distribution of predator population shows the mixture of spot and stripes like
pattern. As we increase the value of time iteration to $5000$, the stripes and
spot like pattern moves toward the irregular spots. The density distribution
of predator population shows the mixture of stripes and spot like pattern at
time evolution $7000$. Finally, at time iteration $10000$, the spots comes
closure to each other and form the cluster of spots. Thus time evolution plays
an important role on system dynamics of spatial model. Also we can observe
that the upper and lower bound of density of predator population decreases
with decrease in time evolution (refer to Colorbar Fig. 13). Fig. 14 depicts
the impact of diffusion coefficient of predator population $D_{2}$ on the
dynamics of spatial model system. Initially at $D_{2}=4.0$, the density is
distributed as mixture of stripe and spot pattern. We observe that stripes and
spot are very near to each other. The spots in the pattern die out and form
the irregular stripes at $D_{2}=6.0$. The population density of predator
distributed in the form of spots at $D_{2}=8.0$. From Fig. 15, one can observe
that the predator population is spatially distributed in the form of mixture
of spots and stripes at prey fear level $k=0.1$ and $k=0.3$. When we increase
the value of fear level as $k=0.6$ and $k=0.7$, the spots come closure to each
other and form irregular stripe like pattern. Also we have found that the
upper and lower density of population decreases with increase in prey fear
level parameter $k$ (refer colorbar of Fig. 15 ). We can also observe that the
density of predator population decreases with increase in fear level. Thus we
can say that fear factor have much impact on spatial distribution of predator
population and coexistence of species. Also from Fig. 2b, one can observe that
the lower value of fear effect is responsible for stability of model system as
higher value of fear level leads the occurrence of Turing instability. Thus
fear level acts as control parameter for system dynamics of spatial model
system. Fig. 16 describes the impact of prey refuge on the spatial
distribution of predator population. Initially we have obtained the mixture of
irregular spot and stripe like pattern and higher density of predator
population is distributed in almost whole domain of consideration. However,
when we increase the value of prey refuge there is no much change in the
pattern and higher density of predator are not distributed in whole domain.
Also the the upper and lower bound of density of predator increases with
increase in prey refuge (refer Fig. 16a-16c).

(a)

(b)

(c)

(d)

Figure 13: Spatial distribution of predator $y$ in $u-v$ plane at different
values of iteration (a) $2000$, (b) $5000$, (c) $7000$ and (d) $10000$ for
model system (11). For the values of other parameters, refer to Eq. (14).

(a)

(b)

(c)

Figure 14: Spatial distribution of predator $y$ in $u-v$ plane at time
iteration $5000$ and different values of diffusion coefficient of predator
$D_{2}$ (a) $D_{2}=4.0$, (b) $D_{2}=6.0$ and (c) $D_{2}=8.0$ for model system
(11). For the values of other parameters, refer to Eq. (14).

(a)

(b)

(c)

(d)

Figure 15: Spatial distribution of predator $y$ in $u-v$ plane at time
iteration $5000$ and different values of fear level parameter (a) $k=0.1$, (b)
$k=0.3$, (c) $k=0.6$ and (d) $k=0.7$ for model system (11). For the values of
other parameters, refer to Eq. (14).

(a)

(b)

(c)

Figure 16: Spatial distribution of predator $y$ in $u-v$ plane at time
iteration $5000$ and different values of prey refuge parameter (a) $m=0.1$,
(b) $m=0.5$ and (c) $m=0.7$ for model system (11). For the values of other
parameters, refer to Eq. (14).

##  8 Discussion

A driving force behind theoretical ecology has been the study of prey-predator
interactions with a variety of significant mechanisms, which may help us to
comprehend ecological systems in the real world more broadly. Prey-predator
models have become more realistic over the last few years as lots of factors
have been introduced into them, like the fear effect, group defence, density-
dependent death rates of predators and prey density-dependent predation
rates[39, 30]. Wang et al. [39] studied a prey-predator model with fear effect
incorporating functional response of Holling type II. He explored that high
level of fear stabilises the system by removing the periodic solutions. The
author also found the existence of subcritical Hopf bifurcation which was not
possible in traditional prey-predator model without fear effect. Zhang et al.
[43] studied a prey-predator model incorporating fear effect, prey refuge, and
Holling type III functional response and found that the increasing value of
fear stabilise the system by eliminating the oscillatory solutions. These
studies have motivated us to investigate a prey-predator interaction
incorporating the fear effect, prey refuge, and Holling type functional
response IV to observe whether the fear effect has a stabilising effect on the
system even when prey species show group defence. In our proposed model, the
functional response of the Monod-Haldane type is meant to reflect the prey’s
group defence. The fear effect is assumed to result in a decrease in the birth
rate of prey. Furthermore, the parameter of predator-taxis sensitivity is used
to relate the group defence and decreased birth rate caused by the fear
effect. Due to the prey’s limited overall time and/or energy, if they devote
more of it to group defence, their birth rate should decline[30].

In mathematical analysis, we have systematically examined all attainable
dynamics of the model system (6) including the positivity and boundedness of
solutions, some scenarios for species persistence and extinction, and local
stability and bifurcation around all potential equilibrium states. We have
identified the adequate parametric criterion for the extinction of prey
species in Theorem 2(i). It demonstrates that the population of prey could
extinct in spite of the absence of predator species. The inclusion of a death
rate of predators caused by intra-predator competition in the considered model
system promotes the extinction of predator species. Theorem 2(ii) gives a
sufficient parametric criterion for the extinction of predator species only.
It demonstrates that as a result of the higher natural death rate of predators
and the prey’s death rate caused by intra-prey competition, predator species
become extinct in the system. However, only prey species survive if they exist
at all. The adequate criterion for the persistence of prey species has been
discussed in Theorem 3(i). We have shown numerically that there exists a
region in $(\alpha,\beta)$ space in which prey species persist in the system
(refer Fig. 12). Further, we discussed the existence of each and every
probable equilibrium along with their stability analysis, analytically and
numerically. According to Theorem 4(i), only the population extinct
equilibrium point exists and is surrounded by a global basin of attraction
when the prey’s natural death rate is higher than its birth rate. In this
circumstance, both species perish. We found a sufficient condition for the
stability of the only predator extinct equilibrium point in Theorem 4(ii). As
a result, there’s a basin of attraction around it, and depending on the size
of the initial population, predator species may become extinct.

Further, we have also explored all possible bifurcations including
transcritical bifurcation, saddle-node bifurcation, and Hopf bifurcation with
respect to different bifurcation parameters such as the level of fear, the
prey’s death rate caused by intra-prey competition, and the prey’s birth rate
as bifurcation parameters. We have numerically validated all of the derived
bifurcations and stability results through one/ two parameter bifurcation
diagram. We have found that our model’s Hopf bifurcation is subcritical in
contrast to traditional prey-predator models, which ignore the cost of fear
and have supercritical Hopf bifurcations in general. We acknowledge the
conclusion of Wang et al. [39] demonstrating numerically that a high amount of
fear can stabilise the prey-predator system by eliminating the occurrence of
periodic solutions. For both bifurcation parameters, the prey’s death rate
caused by the intra-prey competition and the prey’s birth rate, there exists a
region of bi-stability. The region occurs between the predator extinct
equilibrium point and an interior equilibrium point. It may cause the
extinction of predator species depending on the size of initial population. We
have obtained that there is a possibility for the coexistence of the species
for all values of the death rate caused by intraprey competition (refer the
Fig. 3b). In addition, we have also obtained that predators become extinct for
higher values of the prey birth rate ( refer the Fig. 4). We have also
investigated the impact of the level of fear and the birth rate of prey on the
equilibrium prey density and equilibrium predator densities. We have found
some conditions for which both prey and predator equilibrium densities
decrease as the fear level increases. Conditions have also been found under
which both prey and predator equilibrium densities increase as the prey’s
birth rate increases.

A diverse range of patterns has been found in biological systems. However
decades of research, detailed understanding of how these patterns arise still
away from us [56]. In this study, the dynamical behaviour of proposed spatial
model has been discussed analytically as well as numerically. We have also
obtained the condition for occurrence of Turing instability and found that
increasing value of parameter level of fear factor $k$, the system moves
toward spatial instability, whereas the non spatial model shows stable
dynamics (refer to Fig. 2b). The fear factor has an important impact on the
spatially inhomogeneous distribution of the two species [52]. We have observed
that although varying value of fear factor does not much change in the
spatiotemporal pattern but upper bound of density of predator population
decreases with increasing value of this parameter. In case of temporal prey-
predator model, one can observe that as the level of fear factor $k$
increases, the predator density gradually decreases and agrees with those in
[42, 57]. Thus fear factor shows the same impact on the density distribution
for spatial system as in temporal system. The predator density is distributed
in the form of mixture of spot and stripes and higher density of predator is
concentrated near the axes of domain of consideration. Diffusion coefficient
of predator population has also great impact on the spatial density
distribution of the predator population, as varying value of $D_{2}$, we have
found the irregular stripe and mixture of stripe and spots. Time evolution
plays very important role to understand the density distribution of species
over space and variation time leads the cluster of spot pattern.
Ecological Implications: Extinction and persistence are two important
phenomena of population ecology, as they allow population ecologists to
understand whether a species will survive or become extinct in a given
ecosystem. A general scenario is that when the natural death rate of a prey
species exceeds its birth rate, the prey species is wiped out. It means that
irrespective of whether any other species exist or not in the system, the prey
species may go extinct. In this ecosystem, if the growth rate of predator
species is entirely dependent on prey species, they will also go extinct. In
this case, whatever the initial population size of each species be, both the
species will disappear from the system. In a prey-predator interaction, when
the predator density is higher than a certain threshold value, it may drive
the prey to extinction not only by killing them directly but also by inducing
fear in them. They create a substantially high level of fear in the prey
population. Therefore, prey populations change their natural behaviour to keep
themselves safe from predators. In this process, they lose their time for
reproduction. Predators may also go extinct when their densities are higher
than a certain threshold value. In this case, two ecological scenarios may
occur: (i) intra predator competition increases; and (ii) because of the fear
of predators, the prey population makes their group defence tight and does not
allow predators to kill them. In both cases, predator population declines and
predators go extinct in the system. Normally, we observe that as prey number
increases, predator number also increases. However, when prey species show
group defence, predator species may become extinct in the system. Increased
prey numbers strengthen the group defence of prey and make it rare to predate
on themselves. Therefore, predators starve because of the unavailability of
food and eventually die out of the system. Another thing that we observe is
that as intra-prey competition increases, prey number decreases. Because of
the increased intra-prey competition, prey populations are wiped out of the
system, and as a result, predator species are also wiped out. However, in the
presence of group defence, we observe something different. Both populations
may coexist at their lower population densities. Moreover different
bifurcations occur in our proposed model which have the following sense in the
context of ecology: a saddle node bifurcation in our study can lead to
critical shifts in population dynamics, extinction risks, and changes in
ecosystem stability. Understanding these implications is vital for
conservation and ecosystem management efforts, as well as for predicting and
mitigating potential ecological crises. The occurrence of transcritical
bifurcation signifies a change in the ecosystem’s dominant species. Before
this bifurcation, either the predator or prey species may hold greater
abundance and dominance. Following the bifurcation, these roles can swap,
resulting in an altered ecosystem structure. For instance, if the prey had
been dominant previously, the predator could become more abundant and
influential within the ecosystem. Further, the system’s resilience, or its
ability to recover from disturbances, can be influenced by the Hopf
bifurcation. Understanding the conditions under which the system exhibits Hopf
bifurcation can help identify when it is most resilient or vulnerable to
external perturbations, such as habitat loss, climate change, or invasive
species.

In conclusion, the fear effect induced group defence, the prey’s birth rate,
and the prey’s death rate caused by intra-prey competition are probably what
determine the dynamics of the system. The inclusion of the Holling type IV
functional response, the fear factor, and the prey refuge results in rich
dynamics via occurence of different equilibrium points. The stability of these
equilibrium points may alter when local bifurcations take place. The dynamics
of our suggested model under various realistic time delays, various types of
group defence functions, as well as diverse prey refuges, would be fascinating
to observe in near future.

Acknowledgments
The research work of Jai Prakash Tripathi is supported by the Science and
Engineering Research Board (SERB), India [File No. MTR/2022/001028]. The
research work of Satish Kumar Tiwari is supported by RGIPT, Jais, Amethi,
India. The research of YK is partially supported by NSF-DMS (1716802 &
2052820); and The James S. McDonnell665 Foundation 21st Century Science
Initiative in Studying Complex Systems Scholar Award (DOI-27
10.37717/22002047).

## References

  * [1] E. H. Nelson, C. E. Matthews, and J. A. Rosenheim, “Predators reduce prey population growth by inducing changes in prey behavior,” Ecology, vol. 85, no. 7, pp. 1853-1858, 2004.
  * [2] M. Travers, M. Clinchy, L. Zanette, R. Boonstra, and T. D. Williams, “Indirect predator effects on clutch size and the cost of egg production,” Ecology Letters, vol. 13, no. 8, pp. 980-988, 2010.
  * [3] S. Eggers, M. Griesser, M. Nystrand, and J. Ekman, “Predation risk induces changes in nest-site selection and clutch size in the siberian jay,” Proceedings of the Royal Society B: Biological Sciences, vol. 273, no. 1587, pp. 701–706, 2006.
  * [4] L. Y. Zanette, A. F. White, M. C. Allen, and M. Clinchy, “Perceived predation risk reduces the number of offspring songbirds produce per year,” Science, vol. 334, no. 6061, pp. 1398–1401, 2011.
  * [5] S. L. Lima, “Nonlethal effects in the ecology of predator-prey interactions,” Bioscience, vol. 48, no. 1, pp. 25–34, 1998.
  * [6] E. L. Preisser, D. I. Bolnick, and M. F. Benard, “Scared to death? the effects of intimidation and con- sumption in predator–prey interactions,” Ecology, vol. 86, no. 2, pp. 501–509, 2005.
  * [7] B. L. Peckarsky, “Alternative predator avoidance syndromes of stream-dwelling mayfly larvae,” Ecology, vol. 77, no. 6, pp. 1888–1905, 1996.
  * [8] O. J. Schmitz, A. P. Beckerman, and K. M. O’Brien, “Behaviorally mediated trophic cascades: effects of predation risk on food web interactions,” Ecology, vol. 78, no. 5, pp. 1388–1399, 1997.
  * [9] L. B. Kats and L. M. Dill, “The scent of death: chemosensory assessment of predation risk by prey animals,” Ecoscience, vol. 5, no. 3, pp. 361–394, 1998.
  * [10] L. Y. Zanette and M. Clinchy, “Ecology of fear,” Current biology, vol. 29, no. 9, pp. R309–R313, 2019.
  * [11] M. Clinchy, M. J. Sheriff, and L. Y. Zanette, “Predator-induced stress and the ecology of fear,” Functional Ecology, vol. 27, no. 1, pp. 56–65, 2013.
  * [12] J. P. Tripathi, P. S. Mandal, A. Poonia, and V. P. Bajiya, “A widespread interaction between generalist and specialist enemies: The role of intraguild predation and allee effect,” Applied Mathematical Modelling, vol. 89, pp. 105–135, 2021.
  * [13] A. Sih, “Prey refuges and predator-prey stability,” Theoretical Population Biology, vol. 31, no. 1, pp. 1–12, 1987\.
  * [14] V. Bailey, A. Nicholson, and E. Williams, “Interaction between hosts and parasites when some host in-dividuals are more difficult to find than others,” Journal of Theoretical Biology, vol. 3, no. 1, pp. 1–18, 1962\.
  * [15] B. C. Bertram, “Living in groups: predators and prey,” Behavioural ecology: an evolutionary approach, pp. 221–248, 1978.
  * [16] S. Sarwardi, P. K. Mandal, and S. Ray, “Analysis of a competitive prey–predator system with a prey refuge,” Biosystems, vol. 110, no. 3, pp. 133–148, 2012.
  * [17] J. Maynard-Smith, Models in ecology. Cambridge university press, 1974.
  * [18] J. St Amant, “The mathematics of predator-prey interactions,” Ph.D. dissertation, MA Thesis, Univ. of Calif., Santa Barbara, Calif, 1970.
  * [19] W. W. Murdoch and A. Oaten, “Predation and population stability,” in Advances in ecological research. Elsevier, 1975, vol. 9, pp. 1–131.
  * [20] R. R. Vance, “Predation and resource partitioning in one predator–two prey model communities,” The American Naturalist, vol. 112, no. 987, pp. 797–813, 1978.
  * [21] P. Leslie and J. Gower, “The properties of a stochastic model for the predator-prey type of interaction between two species,” Biometrika, vol. 47, no. 3/4, pp. 219–234, 1960.
  * [22] M. P. Hassell and R. M. May, “Stability in insect host-parasite models,” The Journal of Animal Ecology, pp. 693–726, 1973.
  * [23] J. H. Connell, “Community interactions on marine rocky intertidal shores,” Annual review of ecology and systematics, vol. 3, no. 1, pp. 169–192, 1972.
  * [24] R. Upadhyay, “Chaotic dynamics in a three species aquatic population model with holling type ii functional response,” Nonlinear Analysis: Modelling and Control, vol. 13, no. 1, pp. 103–115, 2008.
  * [25] R. K. Upadhyay, R. D. Parshad, K. Antwi-Fordjour, E. Quansah, and S. Kumari, “Global dynamics of stochastic predator–prey model with mutual interference and prey defense,” Journal of Applied Mathematics and Computing, vol. 60, no. 1, pp. 169–190, 2019.
  * [26] S. Kim and K. Antwi-Fordjour, “Prey group defense to predator aggregated induced fear,” The European Physical Journal Plus, vol. 137, no. 6, pp. 1–17, 2022.
  * [27] S. Creel, P. Schuette, and D. Christianson, “Effects of predation risk on group size, vigilance, and foraging behavior in an african ungulate community,” Behavioral Ecology, vol. 25, no. 4, pp. 773–784, 2014.
  * [28] C. R. Boland, “An experimental test of predator detection rates using groups of free-living emus,” Ethology, vol. 109, no. 3, pp. 209–222, 2003.
  * [29] T. Caro, C. Graham, C. Stoner, and J. Vargas, “Adaptive significance of antipredator behaviour in artio-dactyls,” Animal behaviour, vol. 67, no. 2, pp. 205–228, 2004.
  * [30] S. K. Sasmal and Y. Takeuchi, “Dynamics of a predator-prey system with fear and group defense,” Journal of Mathematical Analysis and Applications, vol. 481, no. 1, p. 123471, 2020.
  * [31] D. Blank, “Vigilance, staring and escape running in antipredator behavior of goitered gazelle,” Behavioural processes, vol. 157, pp. 408–416, 2018.
  * [32] C. S. Holling, “The functional response of predators to prey density and its role in mimicry and population regulation,” The Memoirs of the Entomological Society of Canada, vol. 97, no. S45, pp. 5–60, 1965.
  * [33] Q. Cui, Q. Zhang, Z. Qiu, and Z. Hu, “Complex dynamics of a discrete-time predator-prey system with holling iv functional response,” Chaos, Solitons & Fractals, vol. 87, pp. 158–171, 2016.
  * [34] B. Boon and H. Laudelout, “Kinetics of nitrite oxidation by nitrobacter winogradskyi,” Biochemical Jour- nal, vol. 85, no. 3, p. 440, 1962.
  * [35] V. H. Edwards, “The influence of high substrate concentrations on microbial kinetics,” Biotechnology and bioengineering, vol. 12, no. 5, pp. 679–712, 1970.
  * [36] J. F. Andrews, “A mathematical model for the continuous culture of microorganisms utilizing inhibitory substrates,” Biotechnology and bioengineering, vol. 10, no. 6, pp. 707–723, 1968.
  * [37] V. Ajraldi, M. Pittavino, and E. Venturino, “Modeling herd behavior in population systems,” Nonlinear Analysis: Real World Applications, vol. 12, no. 4, pp. 2319–2338, 2011.
  * [38] H. I. Freedman and G. S. Wolkowicz, “Predator-prey systems with group defence: the paradox of enrichment revisited,” Bulletin of Mathematical Biology, vol. 48, no. 5-6, pp. 493–508, 1986.
  * [39] X. Wang, L. Zanette, and X. Zou, “Modelling the fear effect in predator–prey interactions,” Journal of mathematical biology, vol. 73, no. 5, pp. 1179–1204, 2016.
  * [40] Z. Xiao, Z. Li et al., “Stability analysis of a mutual interference predator-prey model with the fear effect,” Journal of Applied Science and Engineering, vol. 22, no. 2, pp. 205–211, 2019.
  * [41] T. K. Kar, “Stability analysis of a prey–predator model incorporating a prey refuge,” Communications in Nonlinear Science and Numerical Simulation, vol. 10, no. 6, pp. 681–691, 2005.
  * [42] H. Zhang, Y. Cai, S. Fu, and W. Wang, “Impact of the fear effect in a prey-predator model incorporating a prey refuge,” Applied Mathematics and Computation, vol. 356, pp. 328–337, 2019.
  * [43] B. Xie and N. Zhang, “Influence of fear effect on a holling type iii prey-predator system with the prey refuge,” AIMS Mathematics, vol. 7, no. 2, pp. 1811–1830, 2022.
  * [44] Y. Dong, D. Wu, C. Shen, and L. Ye, “Influence of fear effect and predator-taxis sensitivity on dynamical behavior of a predator–prey model,” Zeitschrift f$\ddot{u}$r angewandte Mathematik und Physik, vol. 73, no. 1, pp. 1–17, 2022.
  * [45] S. Wu, J. Wang, and J. Shi, “Dynamics and pattern formation of a diffusive predator–prey model with predator-taxis,” Mathematical Models and Methods in Applied Sciences, vol. 28, no. 11, pp. 2275–2312, 2018\.
  * [46] D. Jana, S. Batabyal, and M. Lakshmanan, “Self-diffusion-driven pattern formation in prey–predator sys- tem with complex habitat under fear effect,” The European Physical Journal Plus, vol. 135, no. 11, pp. 1–42, 2020.
  * [47] B. Chakraborty and N. Bairagi, “Complexity in a prey-predator model with prey refuge and diffusion,” Ecological Complexity, vol. 37, pp. 11–23, 2019.
  * [48] L. N. Guin, S. Pal, S. Chakravarty, and S. Djilali, “Pattern dynamics of a reaction-diffusion predator-prey system with both refuge and harvesting,” International Journal of Biomathematics, vol. 14, no. 01, p. 2050084, 2021.
  * [49] L. N. Guin, B. Mondal, and S. Chakravarty, “Existence of spatiotemporal patterns in the reaction–diffusion predator–prey model incorporating prey refuge,” International Journal of Biomathematics, vol. 9, no. 06, p. 1650085, 2016.
  * [50] J. Liu and Y. Kang, “Spatiotemporal dynamics of a diffusive predator–prey model with fear effect,” Non- linear Analysis: Modelling and Control, vol. 27, pp. 1–22, 2022.
  * [51] S. Chen, Z. Liu, and J. Shi, “Nonexistence of nonconstant positive steady states of a diffusive predatorprey model with fear effect, j,” diffusive predator–prey model with fear effect,” Nonlinear Model. Anal, vol. 1, pp. 47–56, 2019.
  * [52] R. Han, L. N. Guin, and B. Dai, “Cross-diffusion-driven pattern formation and selection in a modified leslie–gower predator–prey model with fear effect,” Journal of Biological Systems, vol. 28, no. 01, pp. 27–64, 2020.
  * [53] F. Chen, “On a nonlinear nonautonomous predator–prey model with diffusion and distributed delay,” Journal of Computational and Applied Mathematics, vol. 180, no. 1, pp. 33–49, 2005.
  * [54] L. Perko, Differential equations and dynamical systems. Springer Science & Business Media, 2013, vol. 7.
  * [55] C. Cosner, “Reaction–diffusion equations and ecological modeling,” in Tutorials in mathematical biosciences IV. Springer, 2008, pp. 77–115
  * [56] P. K. Maini and T. E. Woolley, “The turing model for biological pattern formation,” in The dynamics of biological systems. Springer, 2019, pp. 189–204.
  * [57] S. K. Sasmal, “Population dynamics with multiple allee effects induced by fear factors–a mathematical study on prey-predator interactions,” Applied Mathematical Modelling, vol. 64, pp. 1–14, 2018.

##  Appendix A Appendix

###  A.1 Proof of Theorem 1

Positivity

###### Proof.

Considering $x$ and $y$ as functions of $s$, we integrate first equation of
the model system (6) to get,

$\displaystyle\int_{0}^{t}\ \frac{dx}{x}$
$\displaystyle=\int_{0}^{t}\left(\left(\frac{r}{1+k\alpha
y}\right)-d_{1}-d_{2}x-\frac{\left(-m+1\right)\beta
y}{a+\left(-m+1\right)\alpha bx+\left(-m+1\right)^{2}x^{2}}\right)\ ds,$
$\displaystyle\ln{x(t)}-\ln{x(0)}$
$\displaystyle=\int_{0}^{t}\left(\left(\frac{r}{1+k\alpha
y}\right)-d_{1}-d_{2}x-\frac{\left(-m+1\right)\beta
y}{a+\left(-m+1\right)\alpha bx+\left(-m+1\right)^{2}x^{2}}\right)\ ds,$
$\displaystyle x(t)$
$\displaystyle=x(0)exp\left(\int_{0}^{t}\left(\left(\frac{r}{1+k\alpha
y}\right)-d_{1}-d_{2}x-\frac{\left(-m+1\right)\beta
y}{a+\left(-m+1\right)\alpha bx+\left(-m+1\right)^{2}x^{2}}\right)\
ds\right).$

Since $x(0)>0$ and exponential is always positive, hence
$x(t)>0~{}\forall~{}t.$ Similarly, using the second equation of the model
system (6), we get

$\displaystyle\int_{0}^{t}\ \frac{dy}{y}$
$\displaystyle=\int_{0}^{t}\left(-d-ey+\frac{\left(-m+1\right)\beta
cx}{a+\left(1-m\right)\alpha bx+\left(1-m\right)^{2}x^{2}}\right)\ ds,$
$\displaystyle\ln{y(t)}-\ln{y(0)}$
$\displaystyle=\int_{0}^{t}\left(-d-ey+\frac{\left(-m+1\right)\beta
cx}{a+\left(1-m\right)\alpha bx+\left(1-m\right)^{2}x^{2}}\right)\ ds,$
$\displaystyle y(t)$
$\displaystyle=y(0)exp\left(\int_{0}^{t}\left(-d-ey+\frac{\left(-m+1\right)\beta
cx}{a+\left(1-m\right)\alpha bx+\left(1-m\right)^{2}x^{2}}\right)\ ds\right).$

Since $y(0)>0$ and exponential is always positive. Hence
$y(t)>0~{}~{}\forall~{}t$. ∎

Boundedness

###### Proof.

  1. (i)

Using the first equation of the model system (6), we get

$\displaystyle\frac{dx}{dt}$ $\displaystyle\leq\frac{rx}{1+k\alpha
y}-d_{1}x-d_{2}x^{2}\leq rx-
d_{1}x-d_{2}x^{2}=x\left(\left(r-d_{1}\right)-d_{2}x\right).$

Now we use the comparison lemma ([53]) to obtain

$\displaystyle\limsup\limits_{t\to+\infty}x(t)\leq\left(\frac{r-d_{1}}{d_{2}}\right);~{}~{}provided~{}~{}r>d_{1}.$

  2. (ii)

Using the second equation of the model system (6), we get

$\displaystyle\frac{dy}{dt}$ $\displaystyle\leq-dy-
ey^{2}+\frac{c\beta\left(1-m\right)xy}{b\alpha\left(1-m\right)x}=-dy-
ey^{2}+\frac{c\beta y}{b\alpha}=y\left(\frac{c\beta-
db\alpha}{b\alpha}-ey\right).$

Now we use the comparison lemma ([53]) to obtain

$\displaystyle\limsup\limits_{t\to+\infty}y(t)\leq\left(\frac{c\beta-
bd\alpha}{b\alpha e}\right);~{}provided~{}c\beta>bd\alpha.$

∎

###  A.2 Proof of Theorem 4

###### Proof.

In order to analyze the stability of equilibrium points, we use the jacobian
matrix of model system (6) at any point $(x,y)$ which is given by

$J=\left(\begin{array}[]{cc}g_{x}&g_{y}\\\ \\\ h_{x}&h_{y}\end{array}\right),$

where

$\displaystyle g_{x}=$ $\displaystyle\frac{(-m+1)\beta xy\left((-m+1)\alpha
b+2(-m+1)^{2}x\right)}{\left(a+(-m+1)\alpha
bx+(-m+1)^{2}x^{2}\right)^{2}}-\frac{(-m+1)\beta y}{a+(-m+1)\alpha
bx+(-m+1)^{2}x^{2}}-2d_{2}x-d_{1}$ $\displaystyle+\frac{r}{\alpha ky+1},$
$\displaystyle g_{y}=$ $\displaystyle-\frac{(-m+1)\beta x}{a+(-m+1)\alpha
bx+(-m+1)^{2}x^{2}}-\frac{\alpha krx}{(\alpha ky+1)^{2}},$ $\displaystyle
h_{x}=$ $\displaystyle\frac{(m-1)\beta
cy\left(-a+(m-1)^{2}x^{2}\right)}{(a+(m-1)x(-\alpha b+(m-1)x))^{2}},$
$\displaystyle h_{y}=$ $\displaystyle\frac{(-m+1)\beta cx}{a+(-m+1)\alpha
bx+(-m+1)^{2}x^{2}}-d-2ey.$

  1. (i)

The jacobian matrix at $E_{0}(0,0)$ is given by

$J(E_{0})=(J)_{(E_{0})}=\left(\begin{array}[]{cc}-d_{1}+r&0\\\ \\\ 0&-d\\\
\end{array}\right).$

Here, $J(E_{0})$ is diagonal matrix whose eigenvalue are given by diagonal
entries. $\lambda_{1}=r-d_{1}$ and $\lambda_{2}=-d$ are eigenvalues of
$J(E_{0})$. So, $\lambda_{2}$ is always negative because $d>0$. Therefore,
$E_{0}(0,0)$ is saddle if $\lambda_{1}>0$ $(i.e.~{}-d_{1}+r>0)$ and
$E_{0}(0,0)$ is locally asymptotically stable if $\lambda_{1}<0$
$(i.e.~{}-d_{1}+r<0)$. The global asymptotic stability of $E_{0}(0,0)$ under
the condition $r-d_{1}<0$ follows from the Theorem 1 and the local stability
of $E_{0}(0,0)$.

  2. (ii)

The jacobian matrix at $E_{1}(\frac{r-d_{1}}{d_{2}},0)$ is given by

$J(E_{1})=(J)_{E_{1}}=\left(\begin{array}[]{cc}d_{1}-r&-\frac{kr\alpha\left(r-d_{1}\right)}{d_{2}}-\frac{(-m+1)\beta\left(r-d_{1}\right)}{\left(\frac{(-m+1)^{2}\left(r-d_{1}\right){}^{2}}{d_{2}^{2}}+\frac{b(-m+1)\alpha\left(r-d_{1}\right)}{d_{2}}+a\right)d_{2}}\\\
\\\ 0&\frac{\left(-m+1\right)\beta
c\left(r-d_{1}\right)d_{2}-d\left(ad_{2}^{2}+\left(-m+1\right)\alpha
b\left(r-d_{1}\right)d_{2}+\left(-m+1\right)^{2}\left(r-d_{1}\right)^{2}\right)}{\left(ad_{2}^{2}+\left(-m+1\right)\alpha
b\left(r-d_{1}\right)d_{2}+\left(-m+1\right)^{2}\left(r-d_{1}\right)^{2}\right)}\end{array}\right).$

Here, $J(E_{1})$ is an upper triangular matrix, so eigenvalue of $J(E_{1})$
are given by diagonal entries. Hence, $\lambda_{1}=d_{1}-r$ and
$\lambda_{2}=\frac{\left(-m+1\right)\beta
c\left(r-d_{1}\right)d_{2}-d\left(ad_{2}^{2}+\left(-m+1\right)\alpha
b\left(r-d_{1}\right)d_{2}+\left(-m+1\right)^{2}\left(r-d_{1}\right)^{2}\right)}{\left(ad_{2}^{2}+\left(-m+1\right)\alpha
b\left(r-d_{1}\right)d_{2}+\left(-m+1\right)^{2}\left(r-d_{1}\right)^{2}\right)}$
are eigenvalues of $J(E_{1})$. So when $r>d_{1}$, and
$d\left(ad_{2}^{2}+\left(-m+1\right)\alpha
b\left(r-d_{1}\right)d_{2}+\left(-m+1\right)^{2}\left(r-d_{1}\right)^{2}\right)>\left(-m+1\right)\beta
c\left(r-d_{1}\right)d_{2}$, the predator extinct equilibrium point
$E_{1}(\frac{r-d_{1}}{d_{2}},0)$ is asymptotically stable.

∎

###  A.3 Proof of Theorem 11

###### Proof.

At the interior equilibrium $E^{*}$, we take into account the Jacobian matrix
(8). The following is the characteristic equation at $E^{*}$ :

$\lambda^{2}-S(k)\lambda+T(k)=0,$ (15)

where $S=Tr(J(E^{*}))$ and $T=Det(J(E^{*})).$ Eq. (15) has two simple complex
roots

$\lambda_{1,2}=u(k)\pm\iota v(k).$ (16)

The characteristic equation (15) becomes

$\lambda^{2}+T(k)=0,$ (17)

at $k=k_{H}$ since $S(k_{H})=0$. Computing for the roots of Eq. (17) yields
$\lambda_{1,2}=\pm\iota\sqrt{T}$. As a result, we have two eigenvalues that
are entirely imaginary. We also confirm the transversality criterion. In Eq.
(16), for any $k$ in the neighborhood of $k_{H}$ , we have
$u(k)=Re\left[\lambda_{i}(k)\right]=\frac{1}{2}S(k)$ and
$v(k)=\sqrt{T(k)-\frac{[S(k)]^{2}}{4}}$. Thus,

$\displaystyle u(k)=$ $\displaystyle\frac{1}{2}\left(\frac{c(-m+1)\beta
x^{*}}{a+(-m+1)b\alpha
x^{*}+(x^{*})^{2}(-m+1)^{2}}+\frac{\beta(-m+1)x^{*}y^{*}\beta\left((-m+1)b\alpha+2(-m+1)^{2}x^{*}\right)}{\left(a+(-m+1)b\alpha
x^{*}+(x^{*})^{2}(-m+1)^{2}\right)^{2}}\right.$
$\displaystyle-\left.\frac{(-m+1)\beta y^{*}}{a+(-m+1)b\alpha
x^{*}+(-m+1)^{2}(x^{*})^{2}}-2d_{2}x^{*}-d-d_{1}-2ey^{*}+\frac{r}{\alpha
ky^{*}+1}\right).$

When $k=k_{H}$, it is evident that $u(k)=0$. Now we differentiate $u(k)$ with
respect to $k$ to get

$\displaystyle\dfrac{du(k)}{dk}\Big{|}_{k=k_{H}}=$
$\displaystyle\frac{1}{2}\left[\dfrac{d}{dk}\left(\frac{c(-m+1)\beta
x^{*}}{a+(-m+1)b\alpha
x^{*}+(x^{*})^{2}(-m+1)^{2}}+\frac{(-m+1)x^{*}y^{*}\beta\left((-m+1)b\alpha+2(-m+1)^{2}x^{*}\right)}{\left(a+(-m+1)b\alpha
x^{*}+(x^{*})^{2}(-m+1)^{2}\right)^{2}}\right.\right.$
$\displaystyle-\left.\left.\frac{(-m+1)\beta y^{*}}{a+(-m+1)b\alpha
x^{*}+(x^{*})^{2}(-m+1)^{2}}-2d_{2}x^{*}-d-d_{1}-2ey^{*}+\frac{r}{\alpha
ky^{*}+1}\right)\right]\Bigg{|}_{k=k_{H}}$

We point out that since $x^{*}$ and $y^{*}$ can be dependent on the critical
parameter K, it is impossible to determine their exact derivatives because we
lack their explicit expression. Thus, if $\dfrac{d}{dk}\left(S(k)\right)\neq
0$ at $k=k_{H}$, the transversality condition is satisfied. As a result, the
model system (6) experiences Hopf-bifurcation at $k=k_{H}$ at the coexisting
equilibrium point $E^{*}$. ∎

##  Appendix B Appendix

###  B.1 Coefficients of Eq. (7)

$\displaystyle A_{1}=$ $\displaystyle d_{2}(-e)(m-1)^{6}(e-\alpha dk),$
$\displaystyle A_{2}=$ $\displaystyle e(m-1)^{5}\left(\alpha d_{2}(-3\alpha
bdk+3be+\beta ck)-d_{1}(-1+m)(e-\alpha dk)+e(-1+m)r\right),$ $\displaystyle
A_{3}=$
$\displaystyle-e(m-1)^{4}\left(d_{2}\left(3\left(a+\alpha^{2}b^{2}\right)(e-\alpha
dk)+2\alpha^{2}b\beta ck\right)+\alpha d_{1}(m-1)(3\alpha bdk-3be-\beta
ck)+3\alpha be(m-1)r\right),$ $\displaystyle A_{4}=$
$\displaystyle(m-1)^{3}\left(\alpha d_{2}e\left(\beta
ck\left(2a+\alpha^{2}b^{2}\right)+b\left(6a+\alpha^{2}b^{2}\right)(e-\alpha
dk)\right)-d_{1}e(m-1)\left(3\left(a+\alpha^{2}b^{2}\right)(e-\alpha
dk)\right.\right.$ $\displaystyle\left.\left.+2\alpha^{2}b\beta
ck\right)+(m-1)\left(3e^{2}r\left(a+\alpha^{2}b^{2}\right)+\beta d(m-1)(\alpha
dk-e)\right)\right),$ $\displaystyle A_{5}=$
$\displaystyle-(m-1)^{2}\left((m-1)\left(\alpha
be^{2}r\left(6a+\alpha^{2}b^{2}\right)+2\alpha b\beta d(m-1)(\alpha
dk-e)+\beta^{2}c(m-1)(e-2\alpha dk)\right)+ad_{2}e\right.$
$\displaystyle\left.\left(3\left(a+\alpha^{2}b^{2}\right)(e-\alpha
dk)+2\alpha^{2}b\beta ck\right)-\alpha d_{1}e(m-1)\left(\beta
ck\left(2a+\alpha^{2}b^{2}\right)+b\left(6a+\alpha^{2}b^{2}\right)(e-\alpha
dk)\right)\right),$ $\displaystyle A_{6}=$
$\displaystyle(m-1)\left((m-1)\left(3a^{2}e^{2}r+3a\alpha^{2}b^{2}e^{2}r+2a\beta
d(m-1)(\alpha dk-e)+\alpha\beta(m-1)(\alpha bd-\beta c)\right.\right.$
$\displaystyle\left.\left.(\alpha bdk-be-\beta ck)\right)+ae\left(a\alpha
d_{2}(3b(e-\alpha dk)+\beta
ck)-d_{1}(m-1)\left(3\left(a+\alpha^{2}b^{2}\right)(e-\alpha
dk)+2\alpha^{2}b\beta ck\right)\right)\right),$ $\displaystyle A_{7}=$
$\displaystyle-a\left((m-1)\left(3a\alpha be^{2}r+\beta(m-1)(2\alpha bd(\alpha
dk-e)+\beta c(e-2\alpha dk))\right)+ae\left(ad_{2}(e-\alpha dk)+\alpha
d_{1}(m-1)\right.\right.$ $\displaystyle\left.\left.(3\alpha bdk-3be-\beta
ck)\right)\right),$ $\displaystyle A_{8}=$ $\displaystyle
a^{2}\left(ad_{1}e(\alpha dk-e)+ae^{2}r+\beta d(m-1)(\alpha dk-e)\right).$

###  B.2 Expressions of differentiation of functions ‘f’ and ‘g’:

$\displaystyle g_{x}=$ $\displaystyle\frac{(-m+1)\beta
xy\left((-m+1)b\alpha+2(-m+1)^{2}x\right)}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{2}}-\frac{(-m+1)\beta y}{a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}}-2d_{2}x-d_{1}$ $\displaystyle+\frac{r}{\alpha
ky+1},~{}g_{y}=-\frac{(-m+1)\beta x}{a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}}-\frac{\alpha krx}{(\alpha ky+1)^{2}},$ $\displaystyle
g_{xx}=$ $\displaystyle(-(-m+1))\beta
xy\left(\frac{2\left((-m+1)b\alpha+2(-m+1)^{2}x\right)^{2}}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{3}}-\frac{2(-m+1)^{2}}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{2}}\right)$ $\displaystyle-\frac{2(-1+m)\beta
y\left((-m+1)b\alpha+2(-m+1)^{2}x\right)}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{2}}-2d_{2},$ $\displaystyle g_{xy}=$
$\displaystyle\frac{(-m+1)\beta
x\left((-m+1)b\alpha+2(-m+1)^{2}x\right)}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{2}}-\frac{(-m+1)\beta}{a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}}-\frac{\alpha kr}{(\alpha ky+1)^{2}},$ $\displaystyle
g_{yx}=$ $\displaystyle\frac{(-m+1)\beta
x\left((-m+1)b\alpha+2(-m+1)^{2}x\right)}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{2}}-\frac{(-m+1)\beta}{a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}}-\frac{\alpha kr}{(\alpha ky+1)^{2}},$ $\displaystyle
g_{yy}=$ $\displaystyle\frac{2\alpha^{2}k^{2}rx}{(\alpha
ky+1)^{3}},~{}~{}h_{y}=\frac{(-m+1)\beta cx}{a+(-m+1)b\alpha
x+(-m+1)^{2}x^{2}}-d-2ey,~{}h_{yy}=-2e$ $\displaystyle h_{x}=$
$\displaystyle\frac{(-m+1)\beta cy}{a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}}-\frac{(-m+1)\beta
cxy\left((-m+1)b\alpha+2(-m+1)^{2}x\right)}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{2}}$

$\displaystyle h_{xx}=$ $\displaystyle(-m+1)\beta
cxy\left(\frac{2\left((-m+1)b\alpha+2(-m+1)^{2}x\right)^{2}}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{3}}-\frac{2(-m+1)^{2}}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{2}}\right)$ $\displaystyle-\frac{2(-m+1)\beta
cy\left((-m+1)b\alpha+2(-m+1)^{2}x\right)}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{2}},$ $\displaystyle h_{xy}=$
$\displaystyle\frac{c(-m+1)\beta}{a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}}-\frac{c(-m+1)\beta
x\left((-m+1)b\alpha+2(-m+1)^{2}x\right)}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{2}},$ $\displaystyle h_{yx}=$
$\displaystyle\frac{c(-m+1)\beta}{a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}}-\frac{c(-m+1)\beta
x\left((-m+1)b\alpha+2(-m+1)^{2}x\right)}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{2}}.$

### B.3

Expressions of $(\alpha_{ij})$ and $(\beta_{ij})$

$\displaystyle\alpha_{10}=$ $\displaystyle\frac{(-m+1)\beta
xy\left((-m+1)b\alpha+2(-m+1)^{2}x\right)}{\left(a+(-m+1)b\alpha
x+(-m+1)^{2}x^{2}\right)^{2}}-\frac{\beta(-m+1)y}{a+(-m+1)b\alpha
x+(-m+1)^{2}x^{2}}-2d_{2}x-d_{1}+\frac{r}{\alpha ky+1},$
$\displaystyle\alpha_{01}=$ $\displaystyle-\frac{\beta(-m+1)x}{a+(-m+1)\alpha
bx+(-m+1)^{2}x^{2}}-\frac{\alpha krx}{(\alpha ky+1)^{2}},$
$\displaystyle\alpha_{11}=$ $\displaystyle\frac{1}{2}\left(\frac{(-m+1)\beta
x\left(b(-m+1)\alpha+2(-m+1)^{2}x\right)}{\left(a+(-m+1)\alpha
bx+x^{2}(-m+1)^{2}\right)^{2}}-\frac{(-m+1)\beta}{a+(-m+1)\alpha
bx+x^{2}(-m+1)^{2}}-\frac{\alpha kr}{(\alpha ky+1)^{2}}\right),$
$\displaystyle\alpha_{12}=$ $\displaystyle\frac{\alpha^{2}k^{2}r}{3(\alpha
ky+1)^{3}},~{}\alpha_{02}=\frac{\alpha^{2}k^{2}rx}{(\alpha
ky+1)^{3}},~{}\alpha_{03}=-\frac{\alpha^{3}k^{3}rx}{(\alpha ky+1)^{4}},$
$\displaystyle\alpha_{21}=$ $\displaystyle\frac{1}{6}\left((-(-m+1))\beta
x\left(\frac{2\left((-m+1)b\alpha+2(-m+1)^{2}x\right)^{2}}{\left(a+(-m+1)\alpha
bx+x^{2}(-m+1)^{2}\right)^{3}}-\frac{2(-m+1)^{2}}{\left(a+(-m+1)\alpha
bx+x^{2}(-m+1)^{2}\right)^{2}}\right)\right.$
$\displaystyle\left.-\frac{2(m-1)\beta\left((-m+1)b+2(-m+1)^{2}x\right)}{\left(a+(-m+1)\alpha
bx+(-m+1)^{2}x^{2}\right)^{2}}\right),$ $\displaystyle\alpha_{20}=$
$\displaystyle\frac{1}{2}\left(\frac{2(m-1)^{2}\beta y\left(a\alpha
b-3a(m-1)x+(m-1)^{3}x^{3}\right)}{(a+(m-1)x((m-1)x-\alpha
b))^{3}}-2d_{2}\right),$ $\displaystyle\alpha_{30}=$
$\displaystyle-\frac{(m-1)^{3}\beta y\left(a^{2}+a\left(4b(m-1)\alpha
x-\alpha^{2}b^{2}-6x^{2}(m-1)^{2}\right)+x^{4}(m-1)^{4}\right)}{(a+(m-1)x(-\alpha
b+(m-1)x))^{4}},$ $\displaystyle\beta_{10}=$ $\displaystyle\frac{c(-m+1)\beta
y}{a+b(-m+1)\alpha x+x^{2}(-m+1)^{2}}-\frac{c(-m+1)\beta
xy\left(b(-m+1)\alpha+2(-m+1)^{2}x\right)}{\left(a+b(-m+1)\alpha
x+x^{2}(-m+1)^{2}\right)^{2}},$ $\displaystyle\beta_{01}=$
$\displaystyle\frac{c(-m+1)\beta x}{a+b(-m+1)\alpha
x+x^{2}(-m+1)^{2}}-d-2ey,~{}\beta_{03}=0,~{}\beta_{12}=0,~{}\beta_{02}=-e,$
$\displaystyle\beta_{20}=$ $\displaystyle\frac{1}{2}\left(c(-m+1)\beta
xy\left(\frac{2\left(b(-m+1)\alpha+2(-m+1)^{2}x\right)^{2}}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{3}}-\frac{2(-m+1)^{2}}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}x^{2}\right)^{2}}\right)\right.$
$\displaystyle\left.-\frac{2(-m+1)c\beta
y\left((-m+1)b\alpha+2x(-m+1)^{2}\right)}{\left(a+(-m+1)b\alpha
x+x^{2}(-m+1)^{2}\right)^{2}}\right),$ $\displaystyle\beta_{30}=$
$\displaystyle\frac{c(m-1)^{3}\beta y\left(a^{2}+a\left(4b(m-1)\alpha
x-\alpha^{2}b^{2}-6x^{2}(m-1)^{2}\right)+x^{4}(m-1)^{4}\right)}{(a+(m-1)x(-\alpha
b+(m-1)x))^{4}},$ $\displaystyle\beta_{21}=$
$\displaystyle\frac{1}{6}\left(c(-m+1)\beta