# Impact of Fear and Habitat Complexity in a Predator-Prey System with Two Different Shaped Functional Responses: A Comparative Study

Authors: Debgopal Sahoo; Guruprasad Samanta; Manuel De la Sen
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

Habitat complexity or the structural complexity of habitat reduces the available space for interacting species, and subsequently, the encounter rate between the prey and predator is decreased significantly. Different experimental shreds of evidence validate that the presence of the predator strongly affects the physiological behaviour of prey individuals and dramatically reduces their reproduction rate. In this study, we investigate the interplay between the level of fear and the degree of habitat complexity in a predator-prey model with two different shaped functional responses. We, therefore, develop the functional response using the timescale separation method, and the shape of the resulting functional response depends upon the monotonous property of catch rate,gN where N is the prey biomass. Whenever gN increases strictly, a saturating functional response occurs, but for nonmonotonic gN, a dome-shaped functional response arises. For saturating case, it has been revealed that both prey and predator biomass may oscillate for lower levels of fear and a lower degree of habitat complexity. To stabilize this oscillatory behaviour to a coexistence state, we have to adequately increase the level of fear or degree of habitat complexity. However, for dome-shaped case, more complicated dynamics are observed. In this case, coexistence steady state, if exists, may be locally asymptotically stable for a lower degree of habitat complexity, but for intermediate values, the system is capable of producing multiple coexistence steady states with a bistable phenomenon between predator-free steady state and a coexistence steady state. Moreover, if the level of fear is sufficiently low, the system may experience a supercritical or/and subcritical Hopf bifurcation. In the dynamics of parametric disturbance for the degree of habitat complexity parameter, dome-shaped functional response predicts that disturbance may trap the system into a nearest attractor (either a large amplitude stable limit cycle or predator-free steady state); this can be overcome only by a larger alteration, or sometimes it is impossible to overcome (hysteresis phenomena), whereas the saturating-shaped functional response predicts a system resilience. For both the functional responses, a higher degree of habitat complexity always increases the extinction possibility of the predator, and no level of fear can compensate this biodiversity loss.


---

## Body

## 1. Introduction

In ecology, food is one of the basic substances essential for maintaining life and growth of a species. For the purpose of collecting food, members of one species may destroy the members of the other, which establish the most common predator-prey relationship in ecology. From the pioneering invention of Lotka–Volterra model [1, 2], which is a concrete mathematical elucidation of a predator-prey relationship, various rich and prosperous ecological systems have been studied extensively to understand and interpret the dynamical consequences of interacting populations [3–8]. Most of these researches have emphasized qualitative analyses and dynamical behaviour together with stability properties, bifurcation analysis, creation of limit cycle, and permanence. These investigations reveal that the qualitative dynamics of a predator-prey model critically depends on the way of predator’s response towards the changing biomass of its prey [3]. Functional responses or numerical responses are widely used to describe the prey consumption rate as a function of prey density (biomass). In general, C. S. Holling (1959) classifies the functional responses in three types [9, 10]. In type I, consumption rate is linearly increasing with prey density, while in type II, predator is restricted by its capacity to process food, and the number of consumed prey saturates at high levels of prey density. Type III is characterized by a sigmoidal relation in which the consumption rate is positively density-dependent over some regions of prey density. Also, an extensive range of other functional responses, such as Crowley–Martin, Monod–Haldane, and Beddington–DeAngelis, are derived as well [11–14]. Moreover, it has been observed that similar shaped functional responses may not necessarily predict the similar model dynamics [15]. For instance, Aldebert et al. [16] explained the different model predictions like equilibrium versus oscillations, single attractor versus multiple attractors, and resilience versus hysteresis for two similar shaped functional responses: Holling type II and Ivlev’s functional responses. This phenomenon is known as structural sensitivity, which is one of the significant impediment to enhancing the predictive power of ecological models.In spite of being one of the most instinctive concepts in ecology, understanding the role of habitat complexity, which is critically significant for maintaining biodiversity, community persistence, and ecosystem stability, is still an expanding branch in modern research. Populations survive in natural habitats which are spatially or environmentally heterogeneous and could be split up into small patches with various attributes like different levels of resources and environmental capacity. In this article, we define the heterogeneity of habitat as habitat complexity, which refers to the morphological characteristics within a structure itself or the heterogeneity in the arrangement of objects [17]. Habitat complexity is observed in both terrestrial or aquatic ecological systems [18, 19]. For example, in marine habitat, complexity is increased in the presence of oyster and coral reefs, sea grass beds, mangroves, and salt marshes [20]. In lakes, habitat heterogeneity is most commonly present in the form of littoral zone vegetation or depth-gradient diversity [21]. MacArthur and MacArthur (1961) first performed the quantitative analysis of the impact of habitat complexity on species diversity, as the repercussions of which help to explore the local effects of habitat complexity related to increased richness and abundance in more complex habitats [22]. Lots of experimental studies reveal that habitat complexity may decrease the predation rates by reducing the encounter rates or attack rates between predator and prey [23, 24]. In a series of laboratory experiments, Manatunge et al. [25] observed the behaviour of planktivorous fish Pseudorasbora parva (Cyprinidae) and their zooplankton prey Daphnia pulex with artificial vegetation at densities of 0, 350, 700, 1400, 2100, and 2800stems.m−2 and recorded the foraging rates of the fish at different prey densities for all stem densities. They noticed that the foraging efficiency of P. parva decreased significantly with increasing habitat complexity. In a recent study, Bartholomew et al. (2000) proposed two nondimensional indices Ct/At and Sp/Pr that directly measure how the structural complexity interferes with the foraging ability of predator [26]. Ct/At estimates the amount of cover available within a habitat that interferes with predator’s capability to see, and Sp/Pr measures the extent to which the structure interferes with predator’s potentiality to move through the habitat in search of prey. Here, Ct is the total covered area within a habitat, At is the total area of the habitat, Sp is the average interstructural space size, and Pr is the size of the predator. The field experiment with the fish (Fundulus heteroclitus) as predator and amphipods as prey validate that prey survivorship will increase with increasing Ct/At and that survivorship will decrease with increasing Sp/Pr. In mathematical ecology, researchers are incorporating the effect of habitat complexity in a predator-prey system by modifying predator response functions, since this complexity can decrease the encounter coefficient of predators [17, 27, 28].Theoretical ecologists and evolutionary biologists have clarified that the presence of predator may affect the physiological behaviour of prey species to such an extent that it can be more effective in reducing prey biomass than direct predation [29, 30]. However, due to the lack of experimental evidences, this behavioural change of prey species has not been incorporated in a mathematical model. Recently, a field experiment on song sparrows (Melospiza melodia) validated the effect of fear on prey population [31]. In that experiment, it is observed that the number of offspring produced per year is reduced by around 40% due to perception of predation risk. This reduction occurred due to the antipredator activities of birds such as high levels of perceived predation risk which can cause prey to decamp from their habitat temporarily; prey species are worried and frightened to come to an open habitat, so they do not have any free atmosphere to do their daily works including coupling. Therefore, the costs of fear may additionally decrease prey fecundity and survival, and the overall reduction in prey numbers may exceed compared to direct predation alone. In another field experiment on Drosophila melanogaster as prey species and mantid as their predator species, Elliott et al. [32] observed that the presence of mantids decreased the reproductive performance of drosophila in their breeding as well as nonbreeding seasons. Therefore, all these experimental pieces of evidence confirm that the perception of predation risk is persuasive enough to control the community structure and should thus be given greater attention in ecological conservation and management. These experimental findings motivate the researchers to improve the predator-prey model by introducing the fear effects in prey population [33–36]. Wang et al. [37] considered a predator-prey model together with the cost of fear into the prey reproduction, and their analyses exhibit that a high level of fear stabilizes the system by ruling out the periodic oscillations. It has also been proved that the growing birth rate of prey or growing death rate of predator reduces the sensitivity of prey species in perceived predation risk. Wang and Zou [38] proposed a stage-structured predator-prey model with fear effect and adaptive avoidance of predators. Mathematical inspections analyzed the role of fear and the maturation delay between juvenile prey and adult prey in determining the long-term population dynamics. They proved that both high level of fear and strong adaptation of adult prey have a destabilizing effect, while large population of predators has a stabilizing effect on the predator-prey interactions. Most of the previous studies have considered that for high level of fear or larger predator biomass, fear function may affect the reproduction of prey to such an extent that the number of offspring decreases to zero. However, this is not the real scenario observed in various experimental evidences. Therefore, in this work, we assume that even if the level of fear or size of predator population increases to a larger quantity, the prey population will not stop reproduction. It is biologically significant since the prey population become conscious (or aware of) and habituated to the presence of predator.In the above, we have discussed various biological aspects that can make the dynamics of a predator-prey model more complex and rich. Arguably, it is observed that the complexity of a predator-prey model is controlled by the choice of functional response. In this study, we construct a functional response utilizing the quasi-steady-state assumptions in which the dynamical model is simplified based on the idea that the predation processes happen in a much faster timescale than the birth and death processes of predator. We shall show that the monotonic property of catch rate shapes the resulting functional response in different forms. For example, the functional response will be saturating (dome-shaped) if the catch rate is strictly increasing (initially increases and then decreases) with prey density. Further, we incorporate the fear effect in prey reproduction and modify the encounter rate of predator by introducing the habitat complexity. Now, the questions we arise in this paper are the following. (i) How does the choice of two different shaped functional responses affect the qualitative behaviour of a system? (ii) How much the structural sensitivity can influence the dynamics of an ecological model? (iii) How do the habitat complexity and cost of fear jointly affect the population dynamics?Inspection of the predator-prey interaction can be exploited using two major pathways. The most familiar pathway involves numerous methodical steps to characterize the nature of a specific predator-prey relationship by split up the interaction into various fundamental components. Although this approach helps to grasp the insights completely of a specific interaction, it is not fruitful to make that consequence in general statements. For that reason, the second pathway, that is, that of generalization, has been endeavoured. In this paper, we initially analyze the predator-prey model analytically by considering a general fear function and a prey-dependent catch rate and later examine whether this generalization is appropriate for a specific interaction. In recent times, a lot of species are facing the threat of extinction due to rapid climate change around the globe, pollution, or habitat loss, which may affect the ecological balance and biodiversity. So, coexistence in the ecosystem and maintaining biodiversity are the primary intentions of our investigation. In this paper, the study of multiple steady states receives special attention in the context of ecological conservation [39]. It explores the sensitivity of the system behaviour when that system faces some disruptions.We arrange this work in the following manner. In Section2, we upgrade the classical predator-prey model by incorporating the cost of fear (in prey reproduction) and habitat complexity that affects the encounter rate of predator. In this model, instead of using a particular functional response, a general response function has been employed whose graphical shape depends on the monotonic property of the catch rate. All the basic dynamical properties of the model, including positivity and boundedness of the solutions, permanence of the system, and local stability of the possible equilibrium points, are studied in this section. The entire discussion of Section 3 is divided into two sections. Section 3.1 is devoted to analyze the model dynamics with a particular fear function and a monotone increasing catch rate, whereas in Section 3.2, catch rate is taken as a dome-shaped form. In both sections, the feasibility of steady states, the stability analysis, bistability phenomenon (if exists), local bifurcations, and existence of limit cycles are discussed in detail. Our model analysis also includes the effect of fear and habitat complexity in forming the population demography. All of the theoretical findings are convinced through numerical simulations, which show some potential roles of fear effect and habitat complexity that may significantly control the dynamics of predator-prey interplay. In Section 4, we discuss and summarize the biological significance of our analytical findings.


## 2. The Model

Our model is constructed by ameliorating the following classical predator-prey system with logistic growth in prey population:(1)dNdt=rN−dN−aN2−fNP,dPdt=θfNP−mP,with initial conditions N0=N0>0 and P0=P0>0. Here, Nt and Pt are the respective population size of prey and predator at time t. Other parameters bear the following biological interpretations: r, d, and a are the reproduction rate, natural death rate, and intraspecific competition rate, respectively, of the prey population. The term fN represents a functional response that quantifies the variation of captured prey per unit time by one predator on average with changing densities; θ0<θ<1 is the conversion efficiency and θfNP measures the biomass production of P due to predation; m is the mortality rate of predator. The most conventional approach to formulate a functional response, fN, involves various time budgets of predator for searching prey and handling the victim prey. For an alternative approach, we assume that the whole predator community is divided into two different states: searching S and handling H. Therefore, P=S+H. Also, note that the transformation of searching individuals into handling individuals or vice versa occurs on a much faster timescale than the birth and death processes of predator. Hence, the dynamics of these alternations of states are given by(2)dSdt=−cgNS+γH,dHdt=cgNS−γH,with initial conditions S0=S0≥0 and H0=H0≥0. But both S0 and H0 are not simultaneously zero, since P0=S0+H0>0. Biologically, the term gN represents the successful catch rate per searching predator, whereas c represents the attack rate or encounter rate. It is biologically feasible to consider g0=0, that is, in the absence of prey, predator has nothing to catch. The time required for killing, consuming, and digesting prey (or in together handling prey) is symbolized as h=γ−1, where γ is the handling rate. After this handling time, handling individuals again turn back into searching individuals. Now, in this predator-prey system, we incorporate the effect of habitat complexity which is more likely to affect the attack coefficient c [23, 24]. Hence, we replace the attack coefficient c by c1−α, where α0<α<1 is a dimensionless parameter that quantifies the reduction of encounter rate due to habitat complexity and is interpreted as the degree or strength of habitat complexity. Hence, system (2) is modified as(3)dSdt=−c1−αgNS+γH,dHdt=c1−αgNS−γH.Implementing the timescale separation and holdingP=S+H=constant for this faster timescale, a quasi-stationary solution for the searching predator is obtained as(4)S∗=γPγ+c1−αgN=P1+ch1−αgN.We also presume that the predation depends only on available searching individuals, and then the total number of victim prey is given by(5)fNP=c1−αgNS∗=c1−αgNP1+ch1−αgN,which introduces the functional response, fN, as(6)fN=c1−αgN1+ch1−αgN.It should be mentioned that the monotonic property ofgN plays a crucial role in characterizing different shapes of the functional response. For strictly increasing gN, the derived functional response will be a saturating functional response. Otherwise, if gN is growing initially and after some certain prey density, due to group defense, gN will decrease and in that case, a dome-shaped functional response will occur. From the practical perspective, derivation of the functional response in this manner and not to include it directly into the model has significant advantages. For many experiments, one may have to retrieve the catch rate from some artificial laboratory experiments and for that case, it is sufficient to initiate the experiment with an entirely searching (not satiated) predator population into a prey population of different sizes and extract the catch rate depending on the prey population.Various empirical research studies on the predator-prey system reveal that prey’s reproduction rate significantly decreases due to scare of being caught by the predator. Incorporating this biological effect, we modify the model by multiplying the reproduction term by a functionϕk,η,P that accounts for the cost of fear, where k represents the level of fear and η∈0,1 describes the maximum cost of fear. It is biologically relevant to consider the following conditions on ϕk,η,P: ϕ0,η,P=ϕk,η,0=1 means that in the absence of fear or in the absence of predator, the fear function does not affect the growth of prey population; ∂ϕ/∂k<0, ∂ϕ/∂P<0 signify that increase of fear level or enlarging the predator species impacts the negative effect on the reproduction of prey population; limk⟶∞ϕk,η,P=limP⟶∞ϕk,η,P=η implies that even if the level of fear or size of predator population increases infinitely large, the prey population will stress under certain level of fear due to physiological impact when the prey populations are accustomed to fear from predator species. In the experiment of Zanette et al. [31], although the reproduction of offspring songbirds reduces 40% per year due to fear of predator, they have not completely stopped the reproduction and that motivates to consider limk⟶∞ϕk,η,P=limP⟶∞ϕk,η,P=η instead of limk⟶∞ϕk,η,P=limP⟶∞ϕk,η,P=0. It is biologically significant since the prey population is conscious and habituated to the presence of predator. Combining all of these aspects and incorporating these in (1), our resulting model is given by the following system of nonlinear ordinary differential equations:(7)dNdt=rϕk,η,PN−dN−aN2−c1−αgNP1+ch1−αgN,N0=N0,dPdt=θc1−αgNP1+ch1−αgN−mP,P0=P0.


### 2.1. Positivity and Boundedness

To examine the well-posedness of system (7), we prove that every positive solution remains positive and will not grow abruptly for long time interval. Biologically, positivity and uniform boundedness assert that interacting species in our model are ecologically well behaved, and abundance of each species is restricted due to limited resource.Theorem 1.
Each solution of system (7) is positively invariant and uniformly bounded for any time t≥0.Proof.
It can be easily proved that the RHS of system (7) is continuous and locally Lipschitzian in the domain ℝ+2 and hence the solution Nt,Pt of the system exists uniquely on 0,ξ, where 0<ξ≤+∞. Further, it is assumed that limN⟶0+gN/N=g′0≠0. Then, from system (7) with initial conditions N0>0 and P0>0, we have(8)Nt=N0exp∫0trϕk,η,Ps−d−aNs−c1−αgNs/NsPs1+ch1−αgNsds>0,Pt=P0exp∫0tθc1−αgNs1+ch1−αgNs−mds>0.
This proves that the system is positively invariant for any timet>0.
To prove the boundedness of the solution, we consider the functionW=N+1/θP. Then,(9)dWdt=dNdt+1θdPdt=rϕk,η,PN−dN−aN2−mθP≤rN−dN−aN2−mθP.
For eachWt>0, we get(10)dWdt+λW≤rN−d−λN−aN2−m−λθP≤rN−aN2whereλ<minm,d≤−aN−r2a2+r24a≤r24a=Ωsay.
Using the theory of differential inequality [40] for Wt, we have 0<Wt≤Ω/λ1−e−λt+W0e−λt. Taking limit when t⟶∞, we obtain 0<Wt≤Ω/λ. Hence, every positive solution in ℝ+2 confined in the region B=N,P∈ℝ+2:0<N+P/θ≤Ω/λ. This shows that the solution of the system is bounded.


### 2.2. Uniform Persistence

From mathematical point of view, persistence of a system assures that the solutions of that system are always away from zero. Ecologically, it means the sustainability of all the interacting species for a long time interval. Here, we prove the permanence of system (7) directly by utilizing average Lyapunov function [41].Theorem 2.
System (7) is uniformly persistent if (i) r>d and (ii) gr−d/a>m/c1−αθ−mh with θ>mh.Proof.
ForN,P∈ℝ+2, we consider the following average Lyapunov function:(11)ℒN,P=Nσ1Pσ2,where σ1 and σ2 are some positive constants. Now, differentiating with respect to t, we have(12)ΔN,P=ℒ˙ℒ=σ1N˙N+σ2P˙P,=σ1rϕk,η,P−d−aN−c1−αgN/NP1+ch1−αgN+σ2θc1−αgN1+ch1−αgN−m.
The system will be uniformly persistent ifΔN,P>0 at the boundary equilibrium points 0,0 and r−d/a,0 for some σ1>0 and σ2>0. Now, evaluating Δ at the boundary equilibrium points, we obtain(13)Δ0,0=σ1r−d+σ2−m,Δr−da,0=σ2θc1−αgr−d/a1+ch1−αgr−d/a−m.
Therefore, for some positiveσii=1,2, Δ0,0>0 if r>d and Δr−d/a,0>0 if gr−d/a>m/c1−αθ−mh with θ>mh. This completes the proof.Note that the cost of fear in prey reproduction does not affect the persistent of the model system under consideration.


### 2.3. Equilibria and Its Stability

Ecological equilibrium is a dynamical state at which a community of organisms ceases to grow. In predator-prey system, equilibrium occurs at the intersection of their respective nullclines in the nonnegative quadrant. Hence, all the steady states satisfydN/dt=0 and dP/dt=0 simultaneously. Two boundary steady states of this model are given by(14)E0=0,0,E1=r−da,0,providedr>d.Note that higher mortality rate is harmful to any species and may lead that species to extinction, which is biologically intuitive, so throughout this paper, it is assumed thatr>d. Due to biodiversity loss, these steady states are harmful to the effective functioning of the ecosystem. So, coexistence in the ecosystem and maintaining biodiversity are the primary intention of our investigation. Existence of nontrivial steady state E2=N∗,P∗ depends on the form of gN and ϕk,η,P, and this takes the form of(15)gN∗=mc1−αθ−mh,(16)ϕk,η,P∗=dr+aN∗r+mP∗rθN∗.By definition,gN is a catch rate and hence, gN∗>0. Then, for existence of positive N∗, (15) yields(17)θ>mh.IfgN is bounded, then another additional requirement is(18)supgN=M>mc1−αθ−mh.From the biological viewpoint, (17) signifies that the mortality rate of predator needs to be smaller than the product of conversion efficiency θ and handling rate 1/h, which are both related to predation capabilities. However, we presume that searching and handling prey happen on a much smaller timescale than birth and death processes; hence, this inequality is feasible and likely holds. It is worth mentioning that increase of searching rate c fails to compensate for lower handling rate concerning the existence of the interior steady state. Since Z=ϕk,η,P is decreasing in P with ϕk,η,0=1 and Z=d/r+aN∗/r+m/rθN∗P is a straight line with Z intercepts d/r+aN∗/r; they intersect at exactly one point for each N∗, provided(19)N∗<r−da.So, the number of coexistence steady states depends on the form ofgN only: if gN is strictly increasing, then precisely one E2N∗,P∗ exists; otherwise, for dome-shaped gN at-most two E2, namely, E21N1∗,P1∗ and E22N2∗,P2∗ occur. For multiple cases, if we assume N1∗<N2∗, then gN1∗ must lie in the increasing branch, and gN2∗ must lie in the decreasing branch of gN (see Figure 1). Also, it should be noted that the number of coexistence steady states is not affected by the level of fear and from (16), we obtain(20)dP∗dk=−∂ϕ/∂k∂ϕ/∂P∗−m/rθN∗<0,as,∂ϕ∂k<0and∂ϕ∂P∗<0,which means the predator biomass of coexistence equilibrium point decreases due to increased level of fear.Figure 1
Geometrical explanation for existence of coexistence steady state(s). Here,N∗ is the intersection of the curve F1N=gN and the horizontal line F1N=m/c1−αθ−mh. For this particular N∗, P∗ is the meet point of the curve F2P=ϕk,η,P and the line F2P=d/r+aN∗/r+m/rθN∗P. (a, b) The horizontal line F1N=m/c1−αθ−mh intersects the increasing gN at exactly one point provided conditions (17) and (18) are fulfilled, and the corresponding P∗ will be possible if (19) holds. (c, d) Multiple coexistences may occur when gN is a nonmonotonic function. Particular forms of gN and ϕk,η,P are mentioned in the figure. Parameters values are taken as r=0.6,η=0.4,k=0.1,d=0.2,α=0.6,a=0.09,c=0.50.4,h=0.007,θ=0.4,m=0.010.08. Values in the first bracket corresponding to c and m are used for plotting the last two figures.
(a)(b)(c)(d)For the purpose of inspecting local stability around the steady states, we consider the following Jacobian matrix:(21)J=rϕk,η,P−d−2aN−c1−αg′NP1+ch1−αgN2rN∂ϕ∂P−c1−αgN1+ch1−αgNθc1−αg′NP1+ch1−αgN2θc1−αgN1+ch1−αgN−m.Eigenvalues at the trivial equilibrium pointE0 are evaluated as r−d>0 and −m<0. Hence, E0 is always a saddle point. At the semitrivial equilibrium point E1, Jacobian matrix has the eigenvalues λ1=−r−d and λ2=c1−αgr−d/aθ−mh−m/1+ch1−αgr−d/a. If θ<mh, no coexistence steady states exist and the semitrivial solution is locally asymptotically stable (LAS), but this inequality hardly holds. Hence, for stability of E1, we must have(22)gr−da<mc1−αθ−mh.Therefore, we have the following theorems.Theorem 3.
The trivial equilibrium pointE00,0 is always unstable saddle.Theorem 4.
The axial equilibrium pointE1=r−d/a,0 is locally asymptotically stable if gr−d/a<m/c1−αθ−mh.Further, we assume that the coexistence is possible and soN∗<r−d/a. If gN is strictly increasing in N, then exactly one coexistence steady state E2N∗,P∗ exists and gN∗<gr−d/a⇒gr−d/a>m/c1−αθ−mh, which establishes that existence of E2 forces E1 to be an unstable steady state. However, for a dome-shaped predation rate, multiple coexistence steady states, namely, E21N1∗,P1∗ and E22N2∗,P2∗, with N1∗<N2∗ may exist. In that multiple case, N2∗<r−d/a⇒gr−d/a<m/c1−αθ−mh and hence E1 exhibits stable behaviour; that is, bistability may occur between coexistence and semitrivial or axial equilibrium point.Since the explicit form ofE2N∗,P∗ cannot be assessed from (15) and (16), we address the local stability of E2 using the sign of trace and determinant of the Jacobian matrix evaluated at E2. After some algebraic simplifications using (15) and (16), we obtain(23)TrJE2=−aN∗+mθ−c1−αθ−mh2g′N∗N∗P∗θ2N∗,(24)DetJE2=mθ−rN∗∂ϕ∂P|P=P∗c1−αθ−mh2g′N∗P∗θ.Since∂ϕ/∂P<0, for positive determinant, g′N∗>0 must holds. For strictly increasing catch rate, this is always true and hence, for this case, stability of coexistence steady state (exactly one E2N∗,P∗ exists) depends on the sign of TrJE2. Coexistence steady state will be locally asymptotically stable (LAS) if TrJE2<0; that is, aθ2N∗2+c1−αθ−mh2g′N∗N∗P∗>mθP∗. For the either case, that is, for unstable coexistence state, there must be at least one limit cycle that occured [42]. Next, for dome-shaped gN, multiple coexistence steady states E21N1∗,P1∗ and E22N2∗,P2∗ with N1∗<N2∗ may occur, and the corresponding trace and determinant can be easily calculated from (23) and (24). In that case, gN2∗ always lies in the decreasing branch of gN, and then g′N2∗<0, and this implies that the determinant at E22 is less than zero. Therefore, E22 is always a saddle point. Stability behaviour of E21 is the same as E2 in monotone increasing case. Also, for this multiple case, N2∗<r−d/a⇒gr−d/a<m/c1−αθ−mh, and hence, axial steady state is locally asymptotically stable. Therefore, for dome-shaped functional response, the system is capable of producing a bistable phenomenon between E21 and E1.Theorem 5.
(i)
IfgN is strictly increasing, precisely one coexistence steady state, E2N∗,P∗ exists and will be locally asymptotically stable provided TrJE2<0.(ii)
For dome-shapedgN, at most two coexistence steady states E21N1∗,P1∗ and E22N2∗,P2∗ occur. E22 is always unstable, and E21 is locally asymptotically stable if TrJE21<0.We summarize all of these parametric conditions related to the existence of steady states and their stability in Table1.Table 1
Existence of steady states and their stability.
Property ofgNCond. 1Cond. 2(For bounded gN)Consequence of cond. 1 and cond. 2Cond. 3Consequence of cond. 3No. ofE2StabilitygN is strictly increasingθ<mh—NoN∗ exists——0E0: unstableE1: LASθ>mhm/c1−αθ−mh>MNoN∗ exists——0E0: unstableE1: LASm/c1−αθ−mh<MExactly oneN∗ existsN∗>r−d/aNoP∗ exists0E0: unstableE1: LASN∗<r−d/aExactly oneP∗ exists1E0: unstableE1: unstableE2: LAS/limit cyclegN is dome-shapedθ<mh—NoN∗ exists——0E0: unstableE1: LASθ>mhm/c1−αθ−mh>MNoN∗ exists--0E0: unstableE1: LASm/c1−αθ−mh<MAt least oneN∗ (say N1∗) or at-most two N∗ exist (say N1∗, N2∗ with N1∗<N2∗)r−d/a<N1∗<N2∗NoP1∗ and P2∗ exist0E0: unstableE1: LASN1∗<r−d/a<N2∗OnlyP1∗ exists, P2∗ does not exist1E0: unstableE1: unstableE21: LAS/Limit cycleN1∗<N2∗<r−d/aBothP1∗ and P2∗ exist2E0: unstableE1: LASE21: LAS/limit cycleE22: unstableRemark 1.
Forθ=mh, the second equation of model system (7) becomes(25)dPdt=−m1+ch1−αgNP.
Then,dP/dt=0 only when P=0, and therefore no coexistence equilibrium point is possible for θ=mh. Further, we have dP/dt<0 for nonzero P as 0<α<1 and gN≥0. Therefore, the predator biomass decreases with time and the predator-free equilibrium point E1 will be locally asymptotically stable. From the eigenvalue analysis, we can observe that both eigenvalues of the Jacobian matrix evaluated at E1 are negative for θ=mh, which implies that E1 is locally asymptotically stable.Remark 2.
IfsupgN=M=m/c1−αθ−mh and gN is strictly increasing in N, no coexistence equilibrium is possible. Also, gr−d/a<M=m/c1−αθ−mh implies both eigenvalues of the Jacobian matrix evaluated at E1 are negative. Hence, E1 is locally asymptotically stable.
Now, ifgN is dome-shaped and supgN=M=m/c1−αθ−mh, the model system experiences a saddle-node bifurcation at these parametric values. Detailed analysis of this bifurcation has been discussed later by taking a particular form of gN as gN=N/1+N2.


## 2.1. Positivity and Boundedness

To examine the well-posedness of system (7), we prove that every positive solution remains positive and will not grow abruptly for long time interval. Biologically, positivity and uniform boundedness assert that interacting species in our model are ecologically well behaved, and abundance of each species is restricted due to limited resource.Theorem 1.
Each solution of system (7) is positively invariant and uniformly bounded for any time t≥0.Proof.
It can be easily proved that the RHS of system (7) is continuous and locally Lipschitzian in the domain ℝ+2 and hence the solution Nt,Pt of the system exists uniquely on 0,ξ, where 0<ξ≤+∞. Further, it is assumed that limN⟶0+gN/N=g′0≠0. Then, from system (7) with initial conditions N0>0 and P0>0, we have(8)Nt=N0exp∫0trϕk,η,Ps−d−aNs−c1−αgNs/NsPs1+ch1−αgNsds>0,Pt=P0exp∫0tθc1−αgNs1+ch1−αgNs−mds>0.
This proves that the system is positively invariant for any timet>0.
To prove the boundedness of the solution, we consider the functionW=N+1/θP. Then,(9)dWdt=dNdt+1θdPdt=rϕk,η,PN−dN−aN2−mθP≤rN−dN−aN2−mθP.
For eachWt>0, we get(10)dWdt+λW≤rN−d−λN−aN2−m−λθP≤rN−aN2whereλ<minm,d≤−aN−r2a2+r24a≤r24a=Ωsay.
Using the theory of differential inequality [40] for Wt, we have 0<Wt≤Ω/λ1−e−λt+W0e−λt. Taking limit when t⟶∞, we obtain 0<Wt≤Ω/λ. Hence, every positive solution in ℝ+2 confined in the region B=N,P∈ℝ+2:0<N+P/θ≤Ω/λ. This shows that the solution of the system is bounded.


## 2.2. Uniform Persistence

From mathematical point of view, persistence of a system assures that the solutions of that system are always away from zero. Ecologically, it means the sustainability of all the interacting species for a long time interval. Here, we prove the permanence of system (7) directly by utilizing average Lyapunov function [41].Theorem 2.
System (7) is uniformly persistent if (i) r>d and (ii) gr−d/a>m/c1−αθ−mh with θ>mh.Proof.
ForN,P∈ℝ+2, we consider the following average Lyapunov function:(11)ℒN,P=Nσ1Pσ2,where σ1 and σ2 are some positive constants. Now, differentiating with respect to t, we have(12)ΔN,P=ℒ˙ℒ=σ1N˙N+σ2P˙P,=σ1rϕk,η,P−d−aN−c1−αgN/NP1+ch1−αgN+σ2θc1−αgN1+ch1−αgN−m.
The system will be uniformly persistent ifΔN,P>0 at the boundary equilibrium points 0,0 and r−d/a,0 for some σ1>0 and σ2>0. Now, evaluating Δ at the boundary equilibrium points, we obtain(13)Δ0,0=σ1r−d+σ2−m,Δr−da,0=σ2θc1−αgr−d/a1+ch1−αgr−d/a−m.
Therefore, for some positiveσii=1,2, Δ0,0>0 if r>d and Δr−d/a,0>0 if gr−d/a>m/c1−αθ−mh with θ>mh. This completes the proof.Note that the cost of fear in prey reproduction does not affect the persistent of the model system under consideration.


## 2.3. Equilibria and Its Stability

Ecological equilibrium is a dynamical state at which a community of organisms ceases to grow. In predator-prey system, equilibrium occurs at the intersection of their respective nullclines in the nonnegative quadrant. Hence, all the steady states satisfydN/dt=0 and dP/dt=0 simultaneously. Two boundary steady states of this model are given by(14)E0=0,0,E1=r−da,0,providedr>d.Note that higher mortality rate is harmful to any species and may lead that species to extinction, which is biologically intuitive, so throughout this paper, it is assumed thatr>d. Due to biodiversity loss, these steady states are harmful to the effective functioning of the ecosystem. So, coexistence in the ecosystem and maintaining biodiversity are the primary intention of our investigation. Existence of nontrivial steady state E2=N∗,P∗ depends on the form of gN and ϕk,η,P, and this takes the form of(15)gN∗=mc1−αθ−mh,(16)ϕk,η,P∗=dr+aN∗r+mP∗rθN∗.By definition,gN is a catch rate and hence, gN∗>0. Then, for existence of positive N∗, (15) yields(17)θ>mh.IfgN is bounded, then another additional requirement is(18)supgN=M>mc1−αθ−mh.From the biological viewpoint, (17) signifies that the mortality rate of predator needs to be smaller than the product of conversion efficiency θ and handling rate 1/h, which are both related to predation capabilities. However, we presume that searching and handling prey happen on a much smaller timescale than birth and death processes; hence, this inequality is feasible and likely holds. It is worth mentioning that increase of searching rate c fails to compensate for lower handling rate concerning the existence of the interior steady state. Since Z=ϕk,η,P is decreasing in P with ϕk,η,0=1 and Z=d/r+aN∗/r+m/rθN∗P is a straight line with Z intercepts d/r+aN∗/r; they intersect at exactly one point for each N∗, provided(19)N∗<r−da.So, the number of coexistence steady states depends on the form ofgN only: if gN is strictly increasing, then precisely one E2N∗,P∗ exists; otherwise, for dome-shaped gN at-most two E2, namely, E21N1∗,P1∗ and E22N2∗,P2∗ occur. For multiple cases, if we assume N1∗<N2∗, then gN1∗ must lie in the increasing branch, and gN2∗ must lie in the decreasing branch of gN (see Figure 1). Also, it should be noted that the number of coexistence steady states is not affected by the level of fear and from (16), we obtain(20)dP∗dk=−∂ϕ/∂k∂ϕ/∂P∗−m/rθN∗<0,as,∂ϕ∂k<0and∂ϕ∂P∗<0,which means the predator biomass of coexistence equilibrium point decreases due to increased level of fear.Figure 1
Geometrical explanation for existence of coexistence steady state(s). Here,N∗ is the intersection of the curve F1N=gN and the horizontal line F1N=m/c1−αθ−mh. For this particular N∗, P∗ is the meet point of the curve F2P=ϕk,η,P and the line F2P=d/r+aN∗/r+m/rθN∗P. (a, b) The horizontal line F1N=m/c1−αθ−mh intersects the increasing gN at exactly one point provided conditions (17) and (18) are fulfilled, and the corresponding P∗ will be possible if (19) holds. (c, d) Multiple coexistences may occur when gN is a nonmonotonic function. Particular forms of gN and ϕk,η,P are mentioned in the figure. Parameters values are taken as r=0.6,η=0.4,k=0.1,d=0.2,α=0.6,a=0.09,c=0.50.4,h=0.007,θ=0.4,m=0.010.08. Values in the first bracket corresponding to c and m are used for plotting the last two figures.
(a)(b)(c)(d)For the purpose of inspecting local stability around the steady states, we consider the following Jacobian matrix:(21)J=rϕk,η,P−d−2aN−c1−αg′NP1+ch1−αgN2rN∂ϕ∂P−c1−αgN1+ch1−αgNθc1−αg′NP1+ch1−αgN2θc1−αgN1+ch1−αgN−m.Eigenvalues at the trivial equilibrium pointE0 are evaluated as r−d>0 and −m<0. Hence, E0 is always a saddle point. At the semitrivial equilibrium point E1, Jacobian matrix has the eigenvalues λ1=−r−d and λ2=c1−αgr−d/aθ−mh−m/1+ch1−αgr−d/a. If θ<mh, no coexistence steady states exist and the semitrivial solution is locally asymptotically stable (LAS), but this inequality hardly holds. Hence, for stability of E1, we must have(22)gr−da<mc1−αθ−mh.Therefore, we have the following theorems.Theorem 3.
The trivial equilibrium pointE00,0 is always unstable saddle.Theorem 4.
The axial equilibrium pointE1=r−d/a,0 is locally asymptotically stable if gr−d/a<m/c1−αθ−mh.Further, we assume that the coexistence is possible and soN∗<r−d/a. If gN is strictly increasing in N, then exactly one coexistence steady state E2N∗,P∗ exists and gN∗<gr−d/a⇒gr−d/a>m/c1−αθ−mh, which establishes that existence of E2 forces E1 to be an unstable steady state. However, for a dome-shaped predation rate, multiple coexistence steady states, namely, E21N1∗,P1∗ and E22N2∗,P2∗, with N1∗<N2∗ may exist. In that multiple case, N2∗<r−d/a⇒gr−d/a<m/c1−αθ−mh and hence E1 exhibits stable behaviour; that is, bistability may occur between coexistence and semitrivial or axial equilibrium point.Since the explicit form ofE2N∗,P∗ cannot be assessed from (15) and (16), we address the local stability of E2 using the sign of trace and determinant of the Jacobian matrix evaluated at E2. After some algebraic simplifications using (15) and (16), we obtain(23)TrJE2=−aN∗+mθ−c1−αθ−mh2g′N∗N∗P∗θ2N∗,(24)DetJE2=mθ−rN∗∂ϕ∂P|P=P∗c1−αθ−mh2g′N∗P∗θ.Since∂ϕ/∂P<0, for positive determinant, g′N∗>0 must holds. For strictly increasing catch rate, this is always true and hence, for this case, stability of coexistence steady state (exactly one E2N∗,P∗ exists) depends on the sign of TrJE2. Coexistence steady state will be locally asymptotically stable (LAS) if TrJE2<0; that is, aθ2N∗2+c1−αθ−mh2g′N∗N∗P∗>mθP∗. For the either case, that is, for unstable coexistence state, there must be at least one limit cycle that occured [42]. Next, for dome-shaped gN, multiple coexistence steady states E21N1∗,P1∗ and E22N2∗,P2∗ with N1∗<N2∗ may occur, and the corresponding trace and determinant can be easily calculated from (23) and (24). In that case, gN2∗ always lies in the decreasing branch of gN, and then g′N2∗<0, and this implies that the determinant at E22 is less than zero. Therefore, E22 is always a saddle point. Stability behaviour of E21 is the same as E2 in monotone increasing case. Also, for this multiple case, N2∗<r−d/a⇒gr−d/a<m/c1−αθ−mh, and hence, axial steady state is locally asymptotically stable. Therefore, for dome-shaped functional response, the system is capable of producing a bistable phenomenon between E21 and E1.Theorem 5.
(i)
IfgN is strictly increasing, precisely one coexistence steady state, E2N∗,P∗ exists and will be locally asymptotically stable provided TrJE2<0.(ii)
For dome-shapedgN, at most two coexistence steady states E21N1∗,P1∗ and E22N2∗,P2∗ occur. E22 is always unstable, and E21 is locally asymptotically stable if TrJE21<0.We summarize all of these parametric conditions related to the existence of steady states and their stability in Table1.Table 1
Existence of steady states and their stability.
Property ofgNCond. 1Cond. 2(For bounded gN)Consequence of cond. 1 and cond. 2Cond. 3Consequence of cond. 3No. ofE2StabilitygN is strictly increasingθ<mh—NoN∗ exists——0E0: unstableE1: LASθ>mhm/c1−αθ−mh>MNoN∗ exists——0E0: unstableE1: LASm/c1−αθ−mh<MExactly oneN∗ existsN∗>r−d/aNoP∗ exists0E0: unstableE1: LASN∗<r−d/aExactly oneP∗ exists1E0: unstableE1: unstableE2: LAS/limit cyclegN is dome-shapedθ<mh—NoN∗ exists——0E0: unstableE1: LASθ>mhm/c1−αθ−mh>MNoN∗ exists--0E0: unstableE1: LASm/c1−αθ−mh<MAt least oneN∗ (say N1∗) or at-most two N∗ exist (say N1∗, N2∗ with N1∗<N2∗)r−d/a<N1∗<N2∗NoP1∗ and P2∗ exist0E0: unstableE1: LASN1∗<r−d/a<N2∗OnlyP1∗ exists, P2∗ does not exist1E0: unstableE1: unstableE21: LAS/Limit cycleN1∗<N2∗<r−d/aBothP1∗ and P2∗ exist2E0: unstableE1: LASE21: LAS/limit cycleE22: unstableRemark 1.
Forθ=mh, the second equation of model system (7) becomes(25)dPdt=−m1+ch1−αgNP.
Then,dP/dt=0 only when P=0, and therefore no coexistence equilibrium point is possible for θ=mh. Further, we have dP/dt<0 for nonzero P as 0<α<1 and gN≥0. Therefore, the predator biomass decreases with time and the predator-free equilibrium point E1 will be locally asymptotically stable. From the eigenvalue analysis, we can observe that both eigenvalues of the Jacobian matrix evaluated at E1 are negative for θ=mh, which implies that E1 is locally asymptotically stable.Remark 2.
IfsupgN=M=m/c1−αθ−mh and gN is strictly increasing in N, no coexistence equilibrium is possible. Also, gr−d/a<M=m/c1−αθ−mh implies both eigenvalues of the Jacobian matrix evaluated at E1 are negative. Hence, E1 is locally asymptotically stable.
Now, ifgN is dome-shaped and supgN=M=m/c1−αθ−mh, the model system experiences a saddle-node bifurcation at these parametric values. Detailed analysis of this bifurcation has been discussed later by taking a particular form of gN as gN=N/1+N2.


## 3. Model with a Given Fear Function and Catch Rate

This section executes some extensive analysis of model (7) with some particular form of fear function and catch rate. Previously, we have already mentioned that the reproduction term of prey is modified by a fear function ϕk,η,P, where k is the level of fear and η∈0,1 is the maximum cost of fear. Since prey population is conscious and habituated with the presence of predator, we take limk⟶∞ϕk,η,P=limP⟶∞ϕk,η,P=η instead of limk⟶∞ϕk,η,P=limP⟶∞ϕk,η,P=0. Here, we assume the following particular form of ϕ is(26)ϕk,η,P=η+1−η1+kP,which agrees with all the conditions of ϕ as indicated earlier. Diverse biological phenomena, such as saturation and group defense, may arise from the monotonic property of gN, which characterize the different shapes of the functional response. Here, we investigate model (7) for two different forms of gN: one is increasing-shaped, and the other one is dome-shaped.


### 3.1. Strictly Increasing Catch Rate

AssuminggN=N for which the resulting functional response, fN, becomes a saturating functional response and taking ϕk,η,P as (26), we particularise model (7) as(27)dNdt=rη+1−η1+kPN−dN−aN2−c1−αNP1+c1−αhN,dPdt=θc1−αNP1+c1−αhN−mP,with initial population size N0=N0 and P0=P0. It has been already proved that the model (27) is positively invariant and uniformly bounded for any time t≥0. Two boundary equilibrium points that always exist are given by (i) the trivial equilibrium E0=0,0 and (ii) the predator-free equilibrium E1=r−d/a,0. Since gN=N is strictly increasing and unbounded, exactly one coexistence equilibrium, E2N∗,P∗, exists provided θ>mh, and N∗<r−d/a. From (15) and (16), we obtain N∗=m/c1−αθ−mh and P∗ is the positive root of the following quadratic equation:(28)mkθP2+kdN∗+aN∗2+mθ−rkηN∗P+dN∗+aN∗2−rN∗=0.The conditionN∗<r−d/a assures that the above equation has exactly one positive root. On the other hand, if N∗>r−d/a, kdN∗+aN∗2+m/θ−rkηN∗=kdN∗+aN∗2−rN∗+m/θ+rk1−ηN∗>0; hence, no positive P∗ exists. From the existence conditions of E2, we can define an existence region A1 in α−θ parametric plane:(29)A1=0<α<1,0<θ<1:α<α˜,θ>θ˜,where α˜=1−am/cr−dθ−mh and θ˜=mh+am/cr−d.The trivial equilibrium pointE0 is always an unstable saddle. Predator-free equilibrium point E1 will be locally asymptotically stable if r−d/a<m/c1−αθ−mh, which can be deduced directly from (22). Interestingly, stability of E1 ensures the nonexistence of E2. Thus, the stability region of E1 in α−θ parametric plane is given by(30)A2=0<α<1,0<θ<1:α>α˜,θ>θ˜.In the previous section, it has been established that stability of coexistence steady state for strictly increasing catch rate exclusively depends on the sign ofTrJE2, where JE2 is the Jacobian matrix for E2. Substituting gN∗=N∗ in (23) and using the equation of nontrivial predator nullcline, we obtain(31)TrJE2=−aN∗+c21−α2hP∗N∗1+ch1−αN∗2.Hence,E2 is locally asymptotically stable if TrJE2<0. When k=0, one can easily derive that system (27) will be locally asymptotically stable around E2 in the parametric region:(32)A3=0<α<1,0<θ<1:α′<α<α˜,θ>θ′,c>ahr−d,where α′=1−aθ+mh/chr−dθ−mh and θ′=mh1+chr−d/a/chr−d/a−1. In presence of fear, that is, when k>0, we can find α∗ such that α∗<α′ and E2 is locally asymptotically stable for α∗<α′<α<α˜. Hence, fear in prey promotes the possibility of coexistence in environment.


#### 3.1.1. Local Bifurcations: System (27)

Here, we closely observe the change of system behaviour through local bifurcations with respect to the parametersα and k in which various dynamical consequences, such as stability exchange, occurrence of new steady state, and appearance of limit cycle, may happen. Such local bifurcations can also be analyzed from the positions of nontrivial nullclines.(a) Transcritical Bifurcation. In system (27), coexistence steady state E2 exchanges its stability with axial equilibrium point E1 with respect to the bifurcating parameter α; that is, a transcritical bifurcation between E1 and E2 occurs. Earlier, it has been already proved that when α∗<α<α˜, the coexistence state E2 is locally asymptotically stable, but the axial equilibrium E1 is unstable. Whenever α=α˜, two equilibria coincide and exchange their stability. For α>α˜, E2 becomes unstable (in this case, E2 is biologically infeasible) and E1 becomes stable. Here, we prove this bifurcation scenario using Sotomayor’s theorem [43]. For α=α˜, the system exhibits only one axial equilibrium point E1. At E1, the Jacobian matrix is evaluated as(33)JE1;α=α˜=−r−d−rk1−ηr−da−mθ00.One of the eigenvalues ofJE1;α=α˜ is zero. The eigenvectors of JE1;α=α˜ and JE1;α=α˜T corresponding to the zero eigenvalue are, respectively, given by V=v11 and W=01, where v1=−rk1−η/a−m/θr−d. Using the expression for V and W and adopting the similar notations from Perko (2013) [43], one can easily calculate the transversality conditions for transcritical bifurcation as follows:(34)WTFαE1,α=α˜=0,WTDFαE1,α=α˜V=−cr−dθ−mh2aθ≠0,WTD2FE1,α=α˜V,V=−2amθ−mhθr−drk1−ηa+mθr−d≠0.Thus, from Sotomayor’s theorem [43], it can be concluded that the model system undergoes a transcritical bifurcation as parameter α passes through the threshold value α˜.(b) Hopf Bifurcation. We have already observed that coexistence steady state E2 (if exists) may lose their stability when the trace of the corresponding Jacobian matrix JE2 becomes positive. We show that this stability loss occurs via Hopf bifurcation. Taking α as a variable parameter, we assume the characteristic equation of JE2 is λ2−Tαλ+Dα=0, where Tα=TrJE2 and Dα=detJE2. At α=α∗, Tα∗=0 and Dα∗>0. Then, the characteristic equation has purely imaginary roots, which are λ1=iDα∗ and λ2=−iDα∗. Also, it can be shown that the transversality condition for Hopf bifurcation is fulfilled; that is, d/dαReλ|α∗≠0. Therefore, system (27) undergoes a Hopf bifurcation when α=α∗.For quantitative study of the system and to explore the effect of fear and habitat complexity on population dynamics, we execute some numerical study. In this regard, we therefore considerk and α as variable parameters. Other parametric values are hypothetically fixed at r1=2.8,η=0.4,d=0.9,a=0.008,c=0.15,h=0.08,θ=0.3,m=0.6. Figures 2(a)–2(f) are the bifurcation diagrams of system (27) for the bifurcation parameter α. In Figures 2(a) and 2(b), we observe that in absence of fear, both prey and predator oscillate for lower degree of habitat complexity. But as α crosses the lower threshold value αh=α∗=0.515, a stable behaviour is observed around the coexistence steady state and a Hopf bifurcation occurs at α∗=0.515. Numerically, it can be shown that the first Lyapunov coefficient is negative, and hence the Hopf bifurcation is supercritical and the bifurcating limit cycle is stable. For the intermediate values of α, precisely, if 0.515<α<0.9332, the coexistence steady state remains stable; in fact, it is globally asymptotically stable. A transcritical bifurcation occurs at the higher threshold value αtc=α˜=0.9332, and beyond this value, the predator population goes extinct, and the prey population settles to its carrying capacity r−d/a. Now, if we slightly increase the value of fear level to k=0.1, the Hopf bifurcation point drops to αh=α∗=0.165 (Figures 2(c) and 2(d)). Moreover, for higher level of fear, Hopf bifurcation disappears; that is, no oscillatory behaviour occurs (Figures 2(e) and 2(f)). Therefore, we can conclude that the higher level of fear promotes the coexistence in environment.Figure 2
Bifurcation diagram of system (27) for the bifurcation parameter α. (a, b) The change of population biomass for k=0, and in that case, the system undergoes a supercritical Hopf bifurcation at αh=0.515 and a transcritical bifurcation at αtc=0.9332. (c, d) The change of population biomass for k=0.1, in which case the values of αh drops significantly. Moreover, for higher level of fear, that Hopf bifurcation disappears; that is, no oscillatory behaviour is observed in (e, f). In other words, higher level of fear promotes the coexistence. Blue-coloured curve and red-coloured curve indicate the prey/predator biomass of E1=N1,P1=r−d/a,0 and E2=N∗,P∗, respectively. Solid line refers to the stable behaviour, and dashed line refers to the unstable behaviour of the system around any equilibrium point. The dashed-dot curve indicates the max and min of prey/predator population. The set of parameter values taken here is r1=2.8,η=0.4,d=0.9,a=0.008,c=0.15,h=0.08,θ=0.3,m=0.6. (a) Change of prey biomass for the bifurcation parameter α and fixed k = 0. (b) Change of predator biomass for the bifurcation parameter α and fixed k = 0. (c) Change of prey biomass for the bifurcation parameter α and fixed k = 0.1. (d) Change of predator biomass for the bifurcation parameter α and fixed k = 0.1. (e) Change of prey biomass for the bifurcation parameter α and fixed k = 1. (f) Change of predator biomass for the bifurcation parameter α and fixed k = 1.
(a)(b)(c)(d)(e)(f)In Figure3, we have plotted the stability regions of different steady states emerging from system (27) in α−k parametric plane. Here, E2 is unstable and exhibits a stable limit cycle below the Hopf bifurcation curve (red-coloured curve), and it is asymptotically stable in between Hopf bifurcation curve and a transcritical bifurcation curve (blue-coloured curve). Both equilibrium points, E0 and E1, are saddle points below the transcritical bifurcation curve, α=α˜=1−am/cr−dθ−mh, which is independent of fear parameter, k. Hence, if the degree of habitat complexity increases to α˜, whatever the value of k is, predator species always goes extinct, and the prey species survives to its environmental carrying capacity, r−d/a. Also, note that for α>α˜, E2 is biologically infeasible and E0 is a unstable saddle, and that suggests that E1 is globally asymptotically stable. From each region, we take a typical values of α,k and plot the phase trajectories to show the representative behaviour of the solutions of system (27) in that region (Figure 4). Phase trajectories are drawn for different initial biomass 100,20, 130,10, 80,60, and 150,45 to validate that the observed results are global.Figure 3
Stability regions of different equilibrium points of the model system (27) in α−k parametric plane. For lower values of α and k, both prey and predator populations oscillate and settle to a stable limit cycle. As the degree of habitat complexity crosses the lower threshold curve “Hopf bifurcation curve,” the system experiences a supercritical Hopf bifurcation. The coexistence steady state, E2, remains stable; in fact, globally asymptotically stability for the values of α,k belongs to the light blue shaded region. A transcritical bifurcation occurs as the degree of habitat complexity crosses the critical threshold α˜=αtc=0.9332. Other set of parameter values is taken as r1=2.8,η=0.4,d=0.9,a=0.008,c=0.15,h=0.08,θ=0.3,m=0.6..Figure 4
Phase portrait for different representative values ofα and k taken from each region of Figure 3. In each figure, all the steady states, specified by small red-coloured filled circles, are the intersection of the prey nullclines (red-coloured curves) and the predator nullclines (green-coloured curves). Trajectories (blue-coloured curves) are drawn for different initial conditions and all the initial conditions are marked by small blue-coloured circles. The black-coloured orbit represents a stable limit cycle that occurs due to supercritical Hopf bifurcation. Other set of parameter values is r1=2.8,η=0.4,d=0.9,a=0.008,c=0.15,h=0.08,θ=0.3,m=0.6.. (a) E0 (0,0) and E1 (237.5,0) are unstable; limit cycle occurs around E2 (17.7,8.6). (b) E0 (0,0) and E1 (237.5,0) are unstable; E2(39.7,6.2) is globally asymptotically stable. (c) E0 (0,0) is unstable; E1 (237.5,0) is globally asymptotically stable.
(a)(b)(c)


#### 3.1.2. Impact of Fear and Habitat Complexity in Population Biomass

The prey biomass of the coexistence steady stateE2N∗,P∗ is independent of fear parameter k, and so change of fear level does not affect the prey biomass. On the other hand, predator biomass in coexistence steady state that satisfies (28) explicitly depends on k. Also,(35)dP∗dk=−rθ1−ηN∗P∗m1+kP∗2+rkθ1−ηN∗<0,which indicates that predator biomass decreases due to increased level of fear. Biologically, this means that due to the scare of being caught by predator, prey species forage less and that reduces the predation process and consequently predator biomass decreases because of starvation.Now, we differentiateN∗ and P∗ with respect to α to obtain the proper behavioural fluctuations of species biomass in respect of habitat complexity. Derivative of N∗ with respect to α is(36)dN∗dα=mc1−α2θ−mh>0.Hence, the prey population increases as the degree of habitat complexity increases. Now, the differentiation ofP∗ with respect to α gives(37)dP∗dα=mP∗/rθN∗2−a/rm/rθN∗+1−ηk/1+kP∗2dN∗dα.Therefore, variation of predator biomass depends on the sign ofmP∗/rθN∗2−a/r. In Figures 2(a)–2(f), we have plotted the change of species biomass with respect to α and observed that prey population increases as α>αh increases and attains to its carrying capacity at α=αtc. On the other hand, predator biomass initially increases for α>αh, but after a certain value of α, predator biomass decreases steeply and disappears.We can now identify whether the dynamical characteristics emerging from this particular type of system are applicable for model (7) with a general strictly increasing catch rate, gN, that shapes the functional response saturating form. Suppose gN is unbounded increasing function, then the coexistence steady state occurs under the parametric conditions θ>mh and N∗<r−d/a, but for bounded and increasing gN with supgN=M, we must have another additional requirement m/c1−αθ−mh<M for existence of coexistence steady state. Therefore, one can observe different dynamical behaviour, in particular, existence versus nonexistence of interior equilibrium point (coexistence state) for two similar shaped functional responses. However, if the coexistence steady state occurs, all other qualitative behaviour, such as stability, bifurcations, occurrence of limit cycles, impact of fear, and habitat complexity, will be similar to this particular type of model.


### 3.2. Dome-Shaped Catch Rate

Here, we investigate the stability analysis, and the impact of fear and habitat complexity of model (7) in the presence of a dome-shaped catch rate. Particularly, gN is taken as gN=N/1+N2. Then,(38)limN⟶0g′N=1,limN⟶∞g′N=0,which signify that catch rate at lower prey biomass is proportional to the available prey biomass but at higher prey biomass, group defense leads to lower predation success. The maximum catch rate occurs at Nmax=1, then gNmax=1/2; we call it M. Incorporating the above catch rate and the particular fear function, ϕk,η,P, as (26), in the general model (7), we get(39)dNdt=rη+1−η1+kPN−dN−aN2−c1−αNP1+c1−αhN+N2,dPdt=θc1−αNP1+c1−αhN+N2−mP,with initial population size N0=N0 and P0=P0. By Theorem 1, we can say that model (39) is positively invariant and uniformly bounded for any time t≥0. Two boundary steady states (i) E0=0,0 and (ii) E1r−d/a,0 always exist. Since gN=N/1+N2 is dome-shaped, multiple coexistence steady states (in particular, at most two) may occur. The coexistence steady state E2N∗,P∗ satisfies the following two equations:(40)gN∗=N∗1+N∗2=mc1−αθ−mh,(41)1+kηP∗1+kP∗=dr+aN∗r+mP∗rθN∗.Two positiveN∗, say N1∗ and N2∗, exist if m/c1−αθ−mh<gNmax=1/2; that is, α<1−2m/cθ−mh=α˜ (say) and θ>mh+2/c. For each positive N∗, positive P∗ exists if N∗<r−d/a. Therefore, exactly one coexistence steady state E2N∗,P∗ exists if α<1−m/cθ−mhr−d/a+a/r−d=α^ (say) and θ>mh+m/cr−d/a+a/r−d. For α>α^, two possibilities occur: either no E2 exists or exactly two E2E21N1∗,P1∗,E22N2∗,P2∗ exist. Whenever r−d/a<Nmax, no E2 exists and for r−d/a>Nmax, precisely two coexistence state occur if α^<α<α˜. We summarize all of these existence conditions and the number of coexistence states in Table 2.Table 2
Different conditions for existence of interior steady states. Here,α˜=1−2m/cθ−mh and α^=1−m/cθ−mhr−d/a+a/r−d.
Cond. 1Cond. 2Cond. 3No. of coexistence statesr−d/a>Nmaxα˜<α<1—00<α<α˜0<α<α^1α^<α<α˜2r−d/a<Nmaxα˜<α<1—00<α<α˜0<α<α^1α^<α<α˜0The trivial equilibriumE0 is always unstable saddle. The axial equilibrium point E1 is locally asymptotically stable if condition (22) holds, which implies α>α^. Note that whenever exactly one coexistence state occurs, E1 becomes unstable. Now, we investigate the local stability of the coexistence steady state. We can easily get that(42)g′Ni∗=1−Ni∗21+Ni∗22,which is a crucial expression for stability of the nontrivial equilibrium. In particular, a necessary condition for stability is that g′Ni∗>0, which is always true if exactly one coexistence steady state exists. If two coexistence steady states, say E21∗N1∗,P1∗ and E22∗N2∗,P2∗ with N1∗<1<N2∗, exist, then g′N2∗<0; that is, E22∗N2∗,P2∗ must be saddle point. Other E21∗N1∗,P1∗ will be asymptotically stable if TrJE21=−aN1∗+mθ−c1−αθ−mh2g′N1∗P1∗/θ2N1∗<0; that is, α<1−m−aθN1∗θ/cθ−mh2g′N1∗P1∗. Hence, bistability occurs between E21 and E1. On the other hand, if only one coexistence steady state, namely, E2N∗,P∗, exists, then it will be asymptotically stable if TrJE2<0. In Figure 5(a), precisely one coexistence steady state E20.13,1.11 exists and as TrJE2<0, that E2 is locally asymptotically stable. For the initial size 0.1,1.0 or 0.3,1.3, population biomass approaches the coexistence state E20.13,1.11 for long time interval. If we slightly change the initial size to 0.5,1.5 or 2,2, population biomass exhibits some limit cycle behaviour. This happens due to the occurrence of subcritical Hopf bifurcation, which will be discussed in the next section. In Figure 5(c), multiple coexistence states E210.38,1.66 and E222.62,0.53 occur, where E21 is locally asymptotically stable and E22 is unstable saddle. Also, the predator-free steady state E14.44,0 is locally asymptotically stable. For initial size 1,3 or 2,4, population biomass approaches the coexistence steady state E210.38,1.66, and for 1,0.2 or 5,0.3, predator population goes extinct; that is, Nt,Pt tends to E14.44,0 for long time interval. Hence, in this case, bistability between E1 and E21 occurs. In Figure 5(e), no interior steady state exists; hence, the predator-free steady state E14.44,0 is globally asymptotically stable. Different trajectories for initial population biomass 0.1,3.5, 0.5,1.5, 3,4, and 4,2 are drawn in Figure 5. Parameter values are described in Figure 5.Figure 5
Phase portrait for different representative values ofα and k taken from each region of Figure 8. In each figure, all the steady states, specified by small red-coloured filled circles, are the intersection of the prey nullclines (red-coloured curves) and the predator nullclines (green-coloured curves). Trajectories (blue-coloured curves) are drawn for different initial conditions, and all the initial conditions are marked by small blue-coloured circles. The black dotted orbit represents an unstable limit cycle that occurs due to subcritical Hopf bifurcation. Other parameter values are the same as in Figure 6. (a) E0(0,0) and E1(4.44,0) are unstable; E2(0.13,1.11) is locally asymptotically stable. (b) E0(0,0) and E1(4.44,0) are unstable; limit cycle occurs around E2(0.17,2.19). (c) E0(0,0) is unstable; E1(4.44,0) and E210.38,1.66 are locally asymptotically stable; E222.62,0.53 is unstable. (d) E0(0,0) is unstable; E1(4.44,0) is locally asymptotically stable; limit cycle occurs around E210.27,2.41; E223.73,0.66 is unstable. (e) E0(0,0) is unstable; E1(4.44,0) is globally asymptotically stable. No coexistence is possible.
(a)(b)(c)(d)(e)


#### 3.2.1. Local Bifurcations: System (39)

Now, we observe that system (39) undergoes various local bifurcations as the parameters k and α cross through critical thresholds.(a) Transcritical Bifurcation. We prove that the axial equilibrium point E1 exchanges its stability with coexistence steady state E2 as the bifurcating parameter α crosses the threshold value α=α^; that is, a transcritical bifurcation occurs between E1 and E2. For α=α^, the Jacobian matrix of system (39) at E1 is evaluated as follows:(43)JE1;α=α^=−r−d−rk1−ηr−da−mθ00.Clearly one of the eigenvalues is zero. The eigenvectors ofJE1;α=α^ and JE1;α=α^T corresponding to the zero eigenvalue are given by V=v11 and W=01 respectively, where v1=−rk1−η/a−m/θr−d. Using the expressions for V and W, we compute the transversality conditions as follows:(44)WTFαE1,α=α^=0,WTDFαE1,α=α^V=−acθ−mh2r−dθa2+r−d2≠0,WTD2FE1,α=α^V,V=−2mθ−mha2−r−d2θa2+r−d2rk1−ηa+mθr−d≠0,provided r−d/a≠1. All the transversality conditions for transcritical bifurcation are satisfied, and hence by Sotomayor’s theorem, we can conclude that system (39) exhibits transcritical bifurcation for α=α^. In numerical demonstration, we call this α^ as αtc. We conclude this result in the following theorem.Theorem 6.
The predator-free equilibrium pointE1 exhibits a transcritical bifurcation when the bifurcation parameter α crosses the threshold value α^=1−m/cθ−mhr−d/a+a/r−d.(b) Saddle-Node Bifurcation. Let r−d/a>1 and choose α as a bifurcation parameter. Then, the equation obtained from (40), that is, N2−c1−αθ−mh/mN+1=0, has a double root N∗ if α=1−2m/cθ−mh=α˜. Geometrically, this means that for the parametric value α=α˜, the curve y=gN=N/1+N2 touches the straight line y=m/c1−αθ−mh at N=1=N∗, and essentially gN will be the maximum at this point. Hence, g′N∗=0. Now, as r−d/a>1, then corresponding this N∗=1, exactly one P∗ exists. Therefore, for α=α˜, exactly one coexistence steady state 1,P∗ exists with multiplicity 2. The Jacobian matrix of system (39) at E21,P∗ is given by(45)JE2;α=α˜=rη+1−η1+kP∗−d−2a−rk1−η1+kP∗2−mθ00.Hence, one eigenvalue ofJE2;α=α˜ is zero. The eigenvectors of JE2;α=α˜ and JE2;α=α˜T corresponding to the zero eigenvalue are computed as V=v11 and W=01, respectively, where v1=rk1−η/1+kP∗+m/θ/rη+1−η/1+kP∗−d−2a, provided the denominator is not zero. For zero denominator, Bogdanov–Takens bifurcation occurs. Two transversality conditions for saddle-node bifurcation are as follows:(46)WTFαE2,α=α˜=−cθ−mh2P∗2θ≠0,WTD2FE2,α=α˜V,V=−mθ−mhθv12≠0.Therefore, by Sotomayor’s theorem, we can conclude that system (39) undergoes a saddle-node bifurcation for α=α˜. In numerical investigation, we call this α˜ as αsn.Theorem 7.
System (39) undergoes a saddle-node bifurcation as the bifurcation parameter α crosses threshold value α˜=1−2m/cθ−mh.(c) Hopf Bifurcation. We know that the coexistence steady state E2 (if exists) lost its stability when the trace of the corresponding Jacobian matrix JE2 becomes positive. Also, note that for multiple coexistence states, only E21 may lose its stability as E22 is always unstable saddle. We show that this instability occurs via Hopf bifurcation. Let the characteristic equation of JE2 for a variable parameter α is given by λ2−Tαλ+Dα=0, where Tα=TrJE2 and Dα=detJE2. For the suitable choice of other parameter values, we may find a critical value of α, say α∗ from the expression of Tα, such that Tα∗=0 and Dα∗>0. Then, the characteristic equation has purely imaginary eigenvalues, which are λ1=iDα∗ and λ2=−iDα∗. The transversality condition for Hopf bifurcation d/dαReλ|α∗≠0 is also satisfied. Therefore, system (39) undergoes a Hopf bifurcation, when α=α∗.For better visualization of the bifurcation analysis, we execute some numerical study and demonstrate how the fear levelk and the strength of habitat complexity α control the system dynamics. In this regard, we consider the following set of parametric values hypothetically: r=0.6,η=0.4,d=0.2,a=0.09,c=0.5,h=0.007,θ=0.4,m=0.01. Figures 6(a)–6(f) illustrate the qualitative change in system dynamics for the variable parameter α in presence of different levels of fear. In Figures 6(a) and 6(b), we have plotted the change of population biomass and the corresponding bifurcation point for the bifurcation parameter α and the lower fear level, fixed at k=0.1. We can observe that for sufficiently lower strength of habitat complexity, coexistence steady state is locally asymptotically stable. Further, if we increase the value of α to α∗=αh=0.6293, system (39) exhibits limit cycle behaviour and undergoes a Hopf bifurcation around the coexistence state. Also, we calculate the first Lyapunov coefficient numerically which is positive in sign. Hence, the Hopf bifurcation is subcritical, and the bifurcating limit cycle is unstable. Potentially this subcritical case is much more dramatic and dangerous in system biology. Once the strength of habitat complexity crosses the threshold value αh, the trajectories of this system must bounce to the nearest attractor (in this case, it must be either a limit cycle with high amplitude or predator-free steady state) which is far away from the coexistence steady state. Whenever this attractor is a limit cycle, it is impossible to get back the stable coexistence steady state by reversing the value of α (see Figure 7(a)). For either case, that is, when the attractor is the predator-free steady state, we have to reverse the α value adequately to get back the coexistence state (see Figure 7(b)). In other words, population biomass is either trapped in a high amplitude limit cycle through the subcritical Hopf bifurcation or exhibits a looping behaviour. This phenomenon is known as hysteresis. Now at α^=αtc=0.7665, the system undergoes a transcritical bifurcation and the predator-free equilibrium becomes locally asymptotically stable. A saddle-node bifurcation also occurs at a higher degree of habitat complexity, α˜=αsn=0.9, through which the coexistence steady states disappear. Multiple coexistences (E21 and E22) and bistability phenomenon occur in between α^ and α˜, which is specified by the grey-shaded region. Figure 6(a) demonstrates that prey equilibrium biomass increases continuously due to increasing value of habitat complexity, whereas Figure 6(b) shows that predator equilibrium biomass initially increases, but for higher value of habitat complexity, predator species goes extinct.Figure 6
Bifurcation diagram of codimension-1 for the bifurcation parameterα. (a, b) The change of population biomass for the lower fear level, k=0.1, in which case a subcritical Hopf bifurcation at αh=0.6293, a saddle-node bifurcation at αsn=0.9, and a transcritical bifurcation at αtc=0.7665 occur. (c, d) The change of population biomass for the fear level k=0.8, in which case another Hopf bifurcation at αh2=0.89 occurs and that Hopf bifurcation is supercritical. But for larger fear level k=1.5, both Hopf bifurcations disappear in (e, f). Multiple coexistence occurs in the grey-shaded region. Blue-coloured curve and red-coloured curve indicate the prey/predator biomass of E1=N1,P1=r−d/a,0 and E2=N∗,P∗ (or E21N1∗,P1∗ and E22N2∗,P2∗), respectively. Solid lines and dashed lines refer to the stable and unstable behaviour of the system respectively. Dashed-dot curves indicate the max and min of prey/predator population. Other parameter values are fixed at r=0.6,η=0.4,d=0.2,a=0.09,c=0.5,h=0.007,θ=0.4,m=0.01. (a) Change of prey biomass for the bifurcation parameter α and fixed k = 0.1. (b) Change of predator biomass for the bifurcation parameter α and fixed k = 0.1. (c) Change of prey biomass for the bifurcation parameter α and fixed k = 0.8. (d) Change of predator biomass for the bifurcation parameter α and fixed k = 0.8. (e) Change of prey biomass for the bifurcation parameter α and fixed k = 1.5. (f) Change of predator biomass for the parameter α and fixed k = 1.5.
(a)(b)(c)(d)(e)(f)Figure 7
Evidence of hysteresis phenomenon. Using the black arrows, we show how the predator biomass is either trapped in a high amplitude limit cycle through the subcritical Hopf bifurcation (a) or exhibits a hysteresis loop (b). In (a), we takek=0.8, and in (b), k=1.5. Other parameter values are taken from Figure 6. Grey-shaded regions refer to the existence of multiple coexistence equilibrium points. Diagrams for prey populations are qualitatively similar.
(a)(b)To observe the impact of fear in population dynamics, we analyze the model numerically for slightly increasing value ofk. In Figures 6(c) and 6(d), we have plotted the bifurcation diagram and change of population biomass for the bifurcation parameter α and a fixed fear level at k=0.8. In this case, the subcritical Hopf bifurcation occurs at αh1=0.84, which is larger than the previous one αh=0.6293. Also, another Hopf bifurcation occurs at αh2=0.89 through which the coexistence steady state E21 becomes locally asymptotically stable. Numerically, it can be shown that the first Lyapunov coefficient is negative, and hence this Hopf bifurcation is supercritical. For higher fear level, k=1.5, both Hopf bifurcations disappear (see Figures 6(e) and 6(f)) and only transcritical and saddle-node bifurcations remained unchanged. Hence, a higher level of fear promotes the coexistence in environment provided the habitat complexity is sufficiently low. Note that α^=αtc=1−m/cθ−mhr−d/a+a/r−d and α˜=αsn=1−2m/cθ−mh are free from fear parameter, and, therefore, bistability phenomenon must occur between one coexistence steady state and the predator-free steady state. Also, for sufficiently higher degree of habitat complexity, coexistence in environment fails to exist and the system settles into the predator-free steady state.(d) Bogdanov–Takens Bifurcation. For a suitable set of parameter values, system (39) exhibits a saddle-node bifurcation, as well as Hopf bifurcation. For example, we have already proved that for r−d/a>1, two coexistence steady states E21 and E22 appear through saddle-node bifurcation below the threshold value of α˜=1−2m/cθ−mh. Here, E22 is always unstable saddle, and the stability of E21 depends on the sign of TrJE21. Therefore, we can find another critical value of α=α∗ at which TrJE21=0 and E21 becomes unstable through Hopf bifurcation. Geometrically, Bogdanov–Takens (BT) bifurcation is a point in two-parametric bifurcation plane, where the saddle-node bifurcation curve meets the Hopf bifurcation curve, and at that steady state, the corresponding Jacobian matrix has a zero eigenvalue with algebraic multiplicity two.(e) Bautin Bifurcation. The coexistence equilibrium point E21, if it exists, may change its stability through a supercritical or a subcritical Hopf bifurcation. For a supercritical case, a stable limit cycle occurs through Hopf bifurcation and the first Lyapunov coefficient l1>0, whereas for subcritical case, an unstable limit cycle appears and the first Lyapunov coefficient, l1<0. Now, if l1=0, system (39) exhibits a Bautin bifurcation or generalized Hopf (GH) bifurcation in two-parametric bifurcation plane. This bifurcation point separates the subcritical and supercritical Hopf bifurcations, and for neighbouring parametric values, the system exhibits one stable and one unstable limit cycles, which collide and vanish through a saddle-node bifurcation of periodic orbits.In Figure8, we have plotted the stability regions of various steady states emerging from model (39) in α−k parametric plane fixing all other parameter values at r=0.6,η=0.4,d=0.2,a=0.09,c=0.5,h=0.007,θ=0.4,m=0.01 and inspect how the system dynamics depends on the fear level, as well as the habitat complexity. The figure includes a transcritical bifurcation curve (magenta-coloured), a saddle-node bifurcation curve (blue-coloured), and a Hopf bifurcation curve (black-coloured). Both transcritical curve, α=α^=1−m/cθ−mhr−d/a+a/r−d, and saddle-node curve, α=α˜=1−2m/cθ−mh, are free from fear parameter k. Hence, the system always experiences these two local bifurcations for bifurcating parameter α without any restriction on k. One can observe that these bifurcation curves divide the whole parametric region into five small subregions R1−R5. Particularly in R1, unique coexistence steady state occurs which is locally asymptotically stable. That coexistence state eventually turns into an unstable state via Hopf bifurcation, if the parametric values α,k enter R2 region. Numerically, one can prove that this Hopf bifurcation is subcritical. Now, if the degree of habitat complexity crosses the transcritical bifurcation curve and lies between α^ and α˜, multiple coexistence steady state occurs, of which one is unstable saddle and the other one either is locally asymptotically stable (in R3) or exhibits stable/unstable limit cycle (in R4). Further, in these parametric regions, the predator-free steady state is also locally asymptotically stable. Hence, bistability phenomenon occurs in both R3 and R4. Further, for higher value of α, coexistence states disappear through a saddle-node bifurcation and all the trajectories settle into the predator-free steady state (in R5). At αGH,kGH=0.87,0.9, two Hopf bifurcations—one is subcritical and the other one is supercritical—collide and disappear through saddle-node bifurcation of periodic orbits, for which a generalized Hopf bifurcation (GH) or Bautin bifurcation occurs. Also, a Bogdanov–Takens (BT) bifurcation occurs at αBT,kBT=0.9,0.44, where the Hopf bifurcation curve touches the saddle-node bifurcation curve. Five phase portraits for five representative values of α and k, taken from each region R1−R5, are plotted in Figures 5(a)–5(e) and display the typical behaviour of the solution trajectories of system (39) in that region.Figure 8
Stability region of different equilibrium points and two-parametric bifurcation diagram of the model system (39) in α−k parametric plane. One transcritical bifurcation curve (magenta-coloured line), one saddle-node bifurcation curve (blue-coloured line), and one Hopf bifurcation curve (black-coloured curve) split up the whole region into five small subregions R1−R5. In R1, exactly one coexistence steady state E2 occurs which is locally asymptotically stable. In R2, that E2 becomes unstable through subcritical Hopf bifurcation. Multiple coexistence (E21 and E22) occurs whenever α and k together belong to R3 and R4. In R3, E21 is locally asymptotically stable and E22 is unstable saddle. Further, E21 becomes unstable in R4 through Hopf bifurcation. Also, in these two regions, predator-free equilibrium point E1 is locally asymptotically stable. Hence, bistablity phenomenon occurs in R3 and R4. In R5, no coexistence occurs and E1 is globally asymptotically stable. A two-parametric bifurcation Bogdanov–Takens (BT) and generalized Hopf (GH) appear in this diagram. Other parameter values are the same as in Figure 6.


#### 3.2.2. Impact of Fear and Habitat Complexity in Population Biomass

The prey biomass of the coexistence steady stateE2N∗,P∗ is independent of fear parameter k, and so we only discuss here the effects of fear upon the predator biomass. Differentiating (41) with respect to k, we obtain(47)dP∗dk=−rθ1−ηN∗P∗m1+kP∗2+rkθ1−ηN∗<0,which indicates that predator biomass decreases due to increased level of fear; that means due to the scare of being caught by predator, the interaction among predator and prey reduces significantly and consequently predator biomass decreases because of starvation.Now, to obtain the variation of species biomass, we differentiateN∗ and P∗ with respect to α. Here,(48)dN∗dα=mc1−α2θ−mh1−N∗21+N∗22.Hence, ifN∗<1, the prey population increases as the degree of habitat complexity increases, but for N∗>1, prey population decreases. In this particular model, if N∗>1, gN∗ lies in the decreasing branch of gN and hence the corresponding coexistence steady state is unstable saddle, therefore, from stability view point, prey population increases for growing habitat complexity. Now,(49)dP∗dα=mP∗/rθN∗2−a/rm/rθN∗+1−ηk/1+kP∗2dN∗dα.Therefore, variation of predator biomass depends on the sign ofmP∗/rθN∗2−a/r. From Figures 6(a)–6(f), one can observe the variation of species biomass with respect to α. Prey population increases as α increases and attains to its carrying capacity at α=αtc. On the other hand, predator biomass initially increases, but after a certain value of α, predator biomass goes extinct.After some detailed analysis of this particular type of model, we can now recognize the principal characteristics which can be generalized in model (7), where gN is a general dome-shaped catch rate. Let, gN attain its maximum at Nmax and gNmax=M. The number of coexistence steady states depends on the different restrictions of α, which are described in Table 2, with α˜=1−m/cMθ−mh and α^=1−m/cθ−mhgr−d/a. Whenever exactly one coexistence steady state E2 occurs, stability of that E2 depends on the sign of TrJE2. The system experiences a Hopf bifurcation around E2, if TrJE2 becomes zero for some threshold value of α=α∗ or k=k∗. For multiple cases, one coexistence state is either stable or becomes unstable via Hopf bifurcation and the other one is unstable saddle. The system also undergoes a transcritical bifurcation at α=α^. If r−d/a>Nmax, a saddle-node bifurcation occurs at α=α˜. Therefore, a bistability and hysteresis phenomenon appear in between α^ and α˜ for any dome-shaped catch rate.


## 3.1. Strictly Increasing Catch Rate

AssuminggN=N for which the resulting functional response, fN, becomes a saturating functional response and taking ϕk,η,P as (26), we particularise model (7) as(27)dNdt=rη+1−η1+kPN−dN−aN2−c1−αNP1+c1−αhN,dPdt=θc1−αNP1+c1−αhN−mP,with initial population size N0=N0 and P0=P0. It has been already proved that the model (27) is positively invariant and uniformly bounded for any time t≥0. Two boundary equilibrium points that always exist are given by (i) the trivial equilibrium E0=0,0 and (ii) the predator-free equilibrium E1=r−d/a,0. Since gN=N is strictly increasing and unbounded, exactly one coexistence equilibrium, E2N∗,P∗, exists provided θ>mh, and N∗<r−d/a. From (15) and (16), we obtain N∗=m/c1−αθ−mh and P∗ is the positive root of the following quadratic equation:(28)mkθP2+kdN∗+aN∗2+mθ−rkηN∗P+dN∗+aN∗2−rN∗=0.The conditionN∗<r−d/a assures that the above equation has exactly one positive root. On the other hand, if N∗>r−d/a, kdN∗+aN∗2+m/θ−rkηN∗=kdN∗+aN∗2−rN∗+m/θ+rk1−ηN∗>0; hence, no positive P∗ exists. From the existence conditions of E2, we can define an existence region A1 in α−θ parametric plane:(29)A1=0<α<1,0<θ<1:α<α˜,θ>θ˜,where α˜=1−am/cr−dθ−mh and θ˜=mh+am/cr−d.The trivial equilibrium pointE0 is always an unstable saddle. Predator-free equilibrium point E1 will be locally asymptotically stable if r−d/a<m/c1−αθ−mh, which can be deduced directly from (22). Interestingly, stability of E1 ensures the nonexistence of E2. Thus, the stability region of E1 in α−θ parametric plane is given by(30)A2=0<α<1,0<θ<1:α>α˜,θ>θ˜.In the previous section, it has been established that stability of coexistence steady state for strictly increasing catch rate exclusively depends on the sign ofTrJE2, where JE2 is the Jacobian matrix for E2. Substituting gN∗=N∗ in (23) and using the equation of nontrivial predator nullcline, we obtain(31)TrJE2=−aN∗+c21−α2hP∗N∗1+ch1−αN∗2.Hence,E2 is locally asymptotically stable if TrJE2<0. When k=0, one can easily derive that system (27) will be locally asymptotically stable around E2 in the parametric region:(32)A3=0<α<1,0<θ<1:α′<α<α˜,θ>θ′,c>ahr−d,where α′=1−aθ+mh/chr−dθ−mh and θ′=mh1+chr−d/a/chr−d/a−1. In presence of fear, that is, when k>0, we can find α∗ such that α∗<α′ and E2 is locally asymptotically stable for α∗<α′<α<α˜. Hence, fear in prey promotes the possibility of coexistence in environment.


### 3.1.1. Local Bifurcations: System (27)

Here, we closely observe the change of system behaviour through local bifurcations with respect to the parametersα and k in which various dynamical consequences, such as stability exchange, occurrence of new steady state, and appearance of limit cycle, may happen. Such local bifurcations can also be analyzed from the positions of nontrivial nullclines.(a) Transcritical Bifurcation. In system (27), coexistence steady state E2 exchanges its stability with axial equilibrium point E1 with respect to the bifurcating parameter α; that is, a transcritical bifurcation between E1 and E2 occurs. Earlier, it has been already proved that when α∗<α<α˜, the coexistence state E2 is locally asymptotically stable, but the axial equilibrium E1 is unstable. Whenever α=α˜, two equilibria coincide and exchange their stability. For α>α˜, E2 becomes unstable (in this case, E2 is biologically infeasible) and E1 becomes stable. Here, we prove this bifurcation scenario using Sotomayor’s theorem [43]. For α=α˜, the system exhibits only one axial equilibrium point E1. At E1, the Jacobian matrix is evaluated as(33)JE1;α=α˜=−r−d−rk1−ηr−da−mθ00.One of the eigenvalues ofJE1;α=α˜ is zero. The eigenvectors of JE1;α=α˜ and JE1;α=α˜T corresponding to the zero eigenvalue are, respectively, given by V=v11 and W=01, where v1=−rk1−η/a−m/θr−d. Using the expression for V and W and adopting the similar notations from Perko (2013) [43], one can easily calculate the transversality conditions for transcritical bifurcation as follows:(34)WTFαE1,α=α˜=0,WTDFαE1,α=α˜V=−cr−dθ−mh2aθ≠0,WTD2FE1,α=α˜V,V=−2amθ−mhθr−drk1−ηa+mθr−d≠0.Thus, from Sotomayor’s theorem [43], it can be concluded that the model system undergoes a transcritical bifurcation as parameter α passes through the threshold value α˜.(b) Hopf Bifurcation. We have already observed that coexistence steady state E2 (if exists) may lose their stability when the trace of the corresponding Jacobian matrix JE2 becomes positive. We show that this stability loss occurs via Hopf bifurcation. Taking α as a variable parameter, we assume the characteristic equation of JE2 is λ2−Tαλ+Dα=0, where Tα=TrJE2 and Dα=detJE2. At α=α∗, Tα∗=0 and Dα∗>0. Then, the characteristic equation has purely imaginary roots, which are λ1=iDα∗ and λ2=−iDα∗. Also, it can be shown that the transversality condition for Hopf bifurcation is fulfilled; that is, d/dαReλ|α∗≠0. Therefore, system (27) undergoes a Hopf bifurcation when α=α∗.For quantitative study of the system and to explore the effect of fear and habitat complexity on population dynamics, we execute some numerical study. In this regard, we therefore considerk and α as variable parameters. Other parametric values are hypothetically fixed at r1=2.8,η=0.4,d=0.9,a=0.008,c=0.15,h=0.08,θ=0.3,m=0.6. Figures 2(a)–2(f) are the bifurcation diagrams of system (27) for the bifurcation parameter α. In Figures 2(a) and 2(b), we observe that in absence of fear, both prey and predator oscillate for lower degree of habitat complexity. But as α crosses the lower threshold value αh=α∗=0.515, a stable behaviour is observed around the coexistence steady state and a Hopf bifurcation occurs at α∗=0.515. Numerically, it can be shown that the first Lyapunov coefficient is negative, and hence the Hopf bifurcation is supercritical and the bifurcating limit cycle is stable. For the intermediate values of α, precisely, if 0.515<α<0.9332, the coexistence steady state remains stable; in fact, it is globally asymptotically stable. A transcritical bifurcation occurs at the higher threshold value αtc=α˜=0.9332, and beyond this value, the predator population goes extinct, and the prey population settles to its carrying capacity r−d/a. Now, if we slightly increase the value of fear level to k=0.1, the Hopf bifurcation point drops to αh=α∗=0.165 (Figures 2(c) and 2(d)). Moreover, for higher level of fear, Hopf bifurcation disappears; that is, no oscillatory behaviour occurs (Figures 2(e) and 2(f)). Therefore, we can conclude that the higher level of fear promotes the coexistence in environment.Figure 2
Bifurcation diagram of system (27) for the bifurcation parameter α. (a, b) The change of population biomass for k=0, and in that case, the system undergoes a supercritical Hopf bifurcation at αh=0.515 and a transcritical bifurcation at αtc=0.9332. (c, d) The change of population biomass for k=0.1, in which case the values of αh drops significantly. Moreover, for higher level of fear, that Hopf bifurcation disappears; that is, no oscillatory behaviour is observed in (e, f). In other words, higher level of fear promotes the coexistence. Blue-coloured curve and red-coloured curve indicate the prey/predator biomass of E1=N1,P1=r−d/a,0 and E2=N∗,P∗, respectively. Solid line refers to the stable behaviour, and dashed line refers to the unstable behaviour of the system around any equilibrium point. The dashed-dot curve indicates the max and min of prey/predator population. The set of parameter values taken here is r1=2.8,η=0.4,d=0.9,a=0.008,c=0.15,h=0.08,θ=0.3,m=0.6. (a) Change of prey biomass for the bifurcation parameter α and fixed k = 0. (b) Change of predator biomass for the bifurcation parameter α and fixed k = 0. (c) Change of prey biomass for the bifurcation parameter α and fixed k = 0.1. (d) Change of predator biomass for the bifurcation parameter α and fixed k = 0.1. (e) Change of prey biomass for the bifurcation parameter α and fixed k = 1. (f) Change of predator biomass for the bifurcation parameter α and fixed k = 1.
(a)(b)(c)(d)(e)(f)In Figure3, we have plotted the stability regions of different steady states emerging from system (27) in α−k parametric plane. Here, E2 is unstable and exhibits a stable limit cycle below the Hopf bifurcation curve (red-coloured curve), and it is asymptotically stable in between Hopf bifurcation curve and a transcritical bifurcation curve (blue-coloured curve). Both equilibrium points, E0 and E1, are saddle points below the transcritical bifurcation curve, α=α˜=1−am/cr−dθ−mh, which is independent of fear parameter, k. Hence, if the degree of habitat complexity increases to α˜, whatever the value of k is, predator species always goes extinct, and the prey species survives to its environmental carrying capacity, r−d/a. Also, note that for α>α˜, E2 is biologically infeasible and E0 is a unstable saddle, and that suggests that E1 is globally asymptotically stable. From each region, we take a typical values of α,k and plot the phase trajectories to show the representative behaviour of the solutions of system (27) in that region (Figure 4). Phase trajectories are drawn for different initial biomass 100,20, 130,10, 80,60, and 150,45 to validate that the observed results are global.Figure 3
Stability regions of different equilibrium points of the model system (27) in α−k parametric plane. For lower values of α and k, both prey and predator populations oscillate and settle to a stable limit cycle. As the degree of habitat complexity crosses the lower threshold curve “Hopf bifurcation curve,” the system experiences a supercritical Hopf bifurcation. The coexistence steady state, E2, remains stable; in fact, globally asymptotically stability for the values of α,k belongs to the light blue shaded region. A transcritical bifurcation occurs as the degree of habitat complexity crosses the critical threshold α˜=αtc=0.9332. Other set of parameter values is taken as r1=2.8,η=0.4,d=0.9,a=0.008,c=0.15,h=0.08,θ=0.3,m=0.6..Figure 4
Phase portrait for different representative values ofα and k taken from each region of Figure 3. In each figure, all the steady states, specified by small red-coloured filled circles, are the intersection of the prey nullclines (red-coloured curves) and the predator nullclines (green-coloured curves). Trajectories (blue-coloured curves) are drawn for different initial conditions and all the initial conditions are marked by small blue-coloured circles. The black-coloured orbit represents a stable limit cycle that occurs due to supercritical Hopf bifurcation. Other set of parameter values is r1=2.8,η=0.4,d=0.9,a=0.008,c=0.15,h=0.08,θ=0.3,m=0.6.. (a) E0 (0,0) and E1 (237.5,0) are unstable; limit cycle occurs around E2 (17.7,8.6). (b) E0 (0,0) and E1 (237.5,0) are unstable; E2(39.7,6.2) is globally asymptotically stable. (c) E0 (0,0) is unstable; E1 (237.5,0) is globally asymptotically stable.
(a)(b)(c)


### 3.1.2. Impact of Fear and Habitat Complexity in Population Biomass

The prey biomass of the coexistence steady stateE2N∗,P∗ is independent of fear parameter k, and so change of fear level does not affect the prey biomass. On the other hand, predator biomass in coexistence steady state that satisfies (28) explicitly depends on k. Also,(35)dP∗dk=−rθ1−ηN∗P∗m1+kP∗2+rkθ1−ηN∗<0,which indicates that predator biomass decreases due to increased level of fear. Biologically, this means that due to the scare of being caught by predator, prey species forage less and that reduces the predation process and consequently predator biomass decreases because of starvation.Now, we differentiateN∗ and P∗ with respect to α to obtain the proper behavioural fluctuations of species biomass in respect of habitat complexity. Derivative of N∗ with respect to α is(36)dN∗dα=mc1−α2θ−mh>0.Hence, the prey population increases as the degree of habitat complexity increases. Now, the differentiation ofP∗ with respect to α gives(37)dP∗dα=mP∗/rθN∗2−a/rm/rθN∗+1−ηk/1+kP∗2dN∗dα.Therefore, variation of predator biomass depends on the sign ofmP∗/rθN∗2−a/r. In Figures 2(a)–2(f), we have plotted the change of species biomass with respect to α and observed that prey population increases as α>αh increases and attains to its carrying capacity at α=αtc. On the other hand, predator biomass initially increases for α>αh, but after a certain value of α, predator biomass decreases steeply and disappears.We can now identify whether the dynamical characteristics emerging from this particular type of system are applicable for model (7) with a general strictly increasing catch rate, gN, that shapes the functional response saturating form. Suppose gN is unbounded increasing function, then the coexistence steady state occurs under the parametric conditions θ>mh and N∗<r−d/a, but for bounded and increasing gN with supgN=M, we must have another additional requirement m/c1−αθ−mh<M for existence of coexistence steady state. Therefore, one can observe different dynamical behaviour, in particular, existence versus nonexistence of interior equilibrium point (coexistence state) for two similar shaped functional responses. However, if the coexistence steady state occurs, all other qualitative behaviour, such as stability, bifurcations, occurrence of limit cycles, impact of fear, and habitat complexity, will be similar to this particular type of model.


## 3.1.1. Local Bifurcations: System (27)

Here, we closely observe the change of system behaviour through local bifurcations with respect to the parametersα and k in which various dynamical consequences, such as stability exchange, occurrence of new steady state, and appearance of limit cycle, may happen. Such local bifurcations can also be analyzed from the positions of nontrivial nullclines.(a) Transcritical Bifurcation. In system (27), coexistence steady state E2 exchanges its stability with axial equilibrium point E1 with respect to the bifurcating parameter α; that is, a transcritical bifurcation between E1 and E2 occurs. Earlier, it has been already proved that when α∗<α<α˜, the coexistence state E2 is locally asymptotically stable, but the axial equilibrium E1 is unstable. Whenever α=α˜, two equilibria coincide and exchange their stability. For α>α˜, E2 becomes unstable (in this case, E2 is biologically infeasible) and E1 becomes stable. Here, we prove this bifurcation scenario using Sotomayor’s theorem [43]. For α=α˜, the system exhibits only one axial equilibrium point E1. At E1, the Jacobian matrix is evaluated as(33)JE1;α=α˜=−r−d−rk1−ηr−da−mθ00.One of the eigenvalues ofJE1;α=α˜ is zero. The eigenvectors of JE1;α=α˜ and JE1;α=α˜T corresponding to the zero eigenvalue are, respectively, given by V=v11 and W=01, where v1=−rk1−η/a−m/θr−d. Using the expression for V and W and adopting the similar notations from Perko (2013) [43], one can easily calculate the transversality conditions for transcritical bifurcation as follows:(34)WTFαE1,α=α˜=0,WTDFαE1,α=α˜V=−cr−dθ−mh2aθ≠0,WTD2FE1,α=α˜V,V=−2amθ−mhθr−drk1−ηa+mθr−d≠0.Thus, from Sotomayor’s theorem [43], it can be concluded that the model system undergoes a transcritical bifurcation as parameter α passes through the threshold value α˜.(b) Hopf Bifurcation. We have already observed that coexistence steady state E2 (if exists) may lose their stability when the trace of the corresponding Jacobian matrix JE2 becomes positive. We show that this stability loss occurs via Hopf bifurcation. Taking α as a variable parameter, we assume the characteristic equation of JE2 is λ2−Tαλ+Dα=0, where Tα=TrJE2 and Dα=detJE2. At α=α∗, Tα∗=0 and Dα∗>0. Then, the characteristic equation has purely imaginary roots, which are λ1=iDα∗ and λ2=−iDα∗. Also, it can be shown that the transversality condition for Hopf bifurcation is fulfilled; that is, d/dαReλ|α∗≠0. Therefore, system (27) undergoes a Hopf bifurcation when α=α∗.For quantitative study of the system and to explore the effect of fear and habitat complexity on population dynamics, we execute some numerical study. In this regard, we therefore considerk and α as variable parameters. Other parametric values are hypothetically fixed at r1=2.8,η=0.4,d=0.9,a=0.008,c=0.15,h=0.08,θ=0.3,m=0.6. Figures 2(a)–2(f) are the bifurcation diagrams of system (27) for the bifurcation parameter α. In Figures 2(a) and 2(b), we observe that in absence of fear, both prey and predator oscillate for lower degree of habitat complexity. But as α crosses the lower threshold value αh=α∗=0.515, a stable behaviour is observed around the coexistence steady state and a Hopf bifurcation occurs at α∗=0.515. Numerically, it can be shown that the first Lyapunov coefficient is negative, and hence the Hopf bifurcation is supercritical and the bifurcating limit cycle is stable. For the intermediate values of α, precisely, if 0.515<α<0.9332, the coexistence steady state remains stable; in fact, it is globally asymptotically stable. A transcritical bifurcation occurs at the higher threshold value αtc=α˜=0.9332, and beyond this value, the predator population goes extinct, and the prey population settles to its carrying capacity r−d/a. Now, if we slightly increase the value of fear level to k=0.1, the Hopf bifurcation point drops to αh=α∗=0.165 (Figures 2(c) and 2(d)). Moreover, for higher level of fear, Hopf bifurcation disappears; that is, no oscillatory behaviour occurs (Figures 2(e) and 2(f)). Therefore, we can conclude that the higher level of fear promotes the coexistence in environment.Figure 2
Bifurcation diagram of system (27) for the bifurcation parameter α. (a, b) The change of population biomass for k=0, and in that case, the system undergoes a supercritical Hopf bifurcation at αh=0.515 and a transcritical bifurcation at αtc=0.9332. (c, d) The change of population biomass for k=0.1, in which case the values of αh drops significantly. Moreover, for higher level of fear, that Hopf bifurcation disappears; that is, no oscillatory behaviour is observed in (e, f). In other words, higher level of fear promotes the coexistence. Blue-coloured curve and red-coloured curve indicate the prey/predator biomass of E1=N1,P1=r−d/a,0 and E2=N∗,P∗, respectively. Solid line refers to the stable behaviour, and dashed line refers to the unstable behaviour of the system around any equilibrium point. The dashed-dot curve indicates the max and min of prey/predator population. The set of parameter values taken here is r1=2.8,η=0.4,d=0.9,a=0.008,c=0.15,h=0.08,θ=0.3,m=0.6. (a) Change of prey biomass for the bifurcation parameter α and fixed k = 0. (b) Change of predator biomass for the bifurcation parameter α and fixed k = 0. (c) Change of prey biomass for the bifurcation parameter α and fixed k = 0.1. (d) Change of predator biomass for the bifurcation parameter α and fixed k = 0.1. (e) Change of prey biomass for the bifurcation parameter α and fixed k = 1. (f) Change of predator biomass for the bifurcation parameter α and fixed k = 1.
(a)(b)(c)(d)(e)(f)In Figure3, we have plotted the stability regions of different steady states emerging from system (27) in α−k parametric plane. Here, E2 is unstable and exhibits a stable limit cycle below the Hopf bifurcation curve (red-coloured curve), and it is asymptotically stable in between Hopf bifurcation curve and a transcritical bifurcation curve (blue-coloured curve). Both equilibrium points, E0 and E1, are saddle points below the transcritical bifurcation curve, α=α˜=1−am/cr−dθ−mh, which is independent of fear parameter, k. Hence, if the degree of habitat complexity increases to α˜, whatever the value of k is, predator species always goes extinct, and the prey species survives to its environmental carrying capacity, r−d/a. Also, note that for α>α˜, E2 is biologically infeasible and E0 is a unstable saddle, and that suggests that E1 is globally asymptotically stable. From each region, we take a typical values of α,k and plot the phase trajectories to show the representative behaviour of the solutions of system (27) in that region (Figure 4). Phase trajectories are drawn for different initial biomass 100,20, 130,10, 80,60, and 150,45 to validate that the observed results are global.Figure 3
Stability regions of different equilibrium points of the model system (27) in α−k parametric plane. For lower values of α and k, both prey and predator populations oscillate and settle to a stable limit cycle. As the degree of habitat complexity crosses the lower threshold curve “Hopf bifurcation curve,” the system experiences a supercritical Hopf bifurcation. The coexistence steady state, E2, remains stable; in fact, globally asymptotically stability for the values of α,k belongs to the light blue shaded region. A transcritical bifurcation occurs as the degree of habitat complexity crosses the critical threshold α˜=αtc=0.9332. Other set of parameter values is taken as r1=2.8,η=0.4,d=0.9,a=0.008,c=0.15,h=0.08,θ=0.3,m=0.6..Figure 4
Phase portrait for different representative values ofα and k taken from each region of Figure 3. In each figure, all the steady states, specified by small red-coloured filled circles, are the intersection of the prey nullclines (red-coloured curves) and the predator nullclines (green-coloured curves). Trajectories (blue-coloured curves) are drawn for different initial conditions and all the initial conditions are marked by small blue-coloured circles. The black-coloured orbit represents a stable limit cycle that occurs due to supercritical Hopf bifurcation. Other set of parameter values is r1=2.8,η=0.4,d=0.9,a=0.008,c=0.15,h=0.08,θ=0.3,m=0.6.. (a) E0 (0,0) and E1 (237.5,0) are unstable; limit cycle occurs around E2 (17.7,8.6). (b) E0 (0,0) and E1 (237.5,0) are unstable; E2(39.7,6.2) is globally asymptotically stable. (c) E0 (0,0) is unstable; E1 (237.5,0) is globally asymptotically stable.
(a)(b)(c)


## 3.1.2. Impact of Fear and Habitat Complexity in Population Biomass

The prey biomass of the coexistence steady stateE2N∗,P∗ is independent of fear parameter k, and so change of fear level does not affect the prey biomass. On the other hand, predator biomass in coexistence steady state that satisfies (28) explicitly depends on k. Also,(35)dP∗dk=−rθ1−ηN∗P∗m1+kP∗2+rkθ1−ηN∗<0,which indicates that predator biomass decreases due to increased level of fear. Biologically, this means that due to the scare of being caught by predator, prey species forage less and that reduces the predation process and consequently predator biomass decreases because of starvation.Now, we differentiateN∗ and P∗ with respect to α to obtain the proper behavioural fluctuations of species biomass in respect of habitat complexity. Derivative of N∗ with respect to α is(36)dN∗dα=mc1−α2θ−mh>0.Hence, the prey population increases as the degree of habitat complexity increases. Now, the differentiation ofP∗ with respect to α gives(37)dP∗dα=mP∗/rθN∗2−a/rm/rθN∗+1−ηk/1+kP∗2dN∗dα.Therefore, variation of predator biomass depends on the sign ofmP∗/rθN∗2−a/r. In Figures 2(a)–2(f), we have plotted the change of species biomass with respect to α and observed that prey population increases as α>αh increases and attains to its carrying capacity at α=αtc. On the other hand, predator biomass initially increases for α>αh, but after a certain value of α, predator biomass decreases steeply and disappears.We can now identify whether the dynamical characteristics emerging from this particular type of system are applicable for model (7) with a general strictly increasing catch rate, gN, that shapes the functional response saturating form. Suppose gN is unbounded increasing function, then the coexistence steady state occurs under the parametric conditions θ>mh and N∗<r−d/a, but for bounded and increasing gN with supgN=M, we must have another additional requirement m/c1−αθ−mh<M for existence of coexistence steady state. Therefore, one can observe different dynamical behaviour, in particular, existence versus nonexistence of interior equilibrium point (coexistence state) for two similar shaped functional responses. However, if the coexistence steady state occurs, all other qualitative behaviour, such as stability, bifurcations, occurrence of limit cycles, impact of fear, and habitat complexity, will be similar to this particular type of model.


## 3.2. Dome-Shaped Catch Rate

Here, we investigate the stability analysis, and the impact of fear and habitat complexity of model (7) in the presence of a dome-shaped catch rate. Particularly, gN is taken as gN=N/1+N2. Then,(38)limN⟶0g′N=1,limN⟶∞g′N=0,which signify that catch rate at lower prey biomass is proportional to the available prey biomass but at higher prey biomass, group defense leads to lower predation success. The maximum catch rate occurs at Nmax=1, then gNmax=1/2; we call it M. Incorporating the above catch rate and the particular fear function, ϕk,η,P, as (26), in the general model (7), we get(39)dNdt=rη+1−η1+kPN−dN−aN2−c1−αNP1+c1−αhN+N2,dPdt=θc1−αNP1+c1−αhN+N2−mP,with initial population size N0=N0 and P0=P0. By Theorem 1, we can say that model (39) is positively invariant and uniformly bounded for any time t≥0. Two boundary steady states (i) E0=0,0 and (ii) E1r−d/a,0 always exist. Since gN=N/1+N2 is dome-shaped, multiple coexistence steady states (in particular, at most two) may occur. The coexistence steady state E2N∗,P∗ satisfies the following two equations:(40)gN∗=N∗1+N∗2=mc1−αθ−mh,(41)1+kηP∗1+kP∗=dr+aN∗r+mP∗rθN∗.Two positiveN∗, say N1∗ and N2∗, exist if m/c1−αθ−mh<gNmax=1/2; that is, α<1−2m/cθ−mh=α˜ (say) and θ>mh+2/c. For each positive N∗, positive P∗ exists if N∗<r−d/a. Therefore, exactly one coexistence steady state E2N∗,P∗ exists if α<1−m/cθ−mhr−d/a+a/r−d=α^ (say) and θ>mh+m/cr−d/a+a/r−d. For α>α^, two possibilities occur: either no E2 exists or exactly two E2E21N1∗,P1∗,E22N2∗,P2∗ exist. Whenever r−d/a<Nmax, no E2 exists and for r−d/a>Nmax, precisely two coexistence state occur if α^<α<α˜. We summarize all of these existence conditions and the number of coexistence states in Table 2.Table 2
Different conditions for existence of interior steady states. Here,α˜=1−2m/cθ−mh and α^=1−m/cθ−mhr−d/a+a/r−d.
Cond. 1Cond. 2Cond. 3No. of coexistence statesr−d/a>Nmaxα˜<α<1—00<α<α˜0<α<α^1α^<α<α˜2r−d/a<Nmaxα˜<α<1—00<α<α˜0<α<α^1α^<α<α˜0The trivial equilibriumE0 is always unstable saddle. The axial equilibrium point E1 is locally asymptotically stable if condition (22) holds, which implies α>α^. Note that whenever exactly one coexistence state occurs, E1 becomes unstable. Now, we investigate the local stability of the coexistence steady state. We can easily get that(42)g′Ni∗=1−Ni∗21+Ni∗22,which is a crucial expression for stability of the nontrivial equilibrium. In particular, a necessary condition for stability is that g′Ni∗>0, which is always true if exactly one coexistence steady state exists. If two coexistence steady states, say E21∗N1∗,P1∗ and E22∗N2∗,P2∗ with N1∗<1<N2∗, exist, then g′N2∗<0; that is, E22∗N2∗,P2∗ must be saddle point. Other E21∗N1∗,P1∗ will be asymptotically stable if TrJE21=−aN1∗+mθ−c1−αθ−mh2g′N1∗P1∗/θ2N1∗<0; that is, α<1−m−aθN1∗θ/cθ−mh2g′N1∗P1∗. Hence, bistability occurs between E21 and E1. On the other hand, if only one coexistence steady state, namely, E2N∗,P∗, exists, then it will be asymptotically stable if TrJE2<0. In Figure 5(a), precisely one coexistence steady state E20.13,1.11 exists and as TrJE2<0, that E2 is locally asymptotically stable. For the initial size 0.1,1.0 or 0.3,1.3, population biomass approaches the coexistence state E20.13,1.11 for long time interval. If we slightly change the initial size to 0.5,1.5 or 2,2, population biomass exhibits some limit cycle behaviour. This happens due to the occurrence of subcritical Hopf bifurcation, which will be discussed in the next section. In Figure 5(c), multiple coexistence states E210.38,1.66 and E222.62,0.53 occur, where E21 is locally asymptotically stable and E22 is unstable saddle. Also, the predator-free steady state E14.44,0 is locally asymptotically stable. For initial size 1,3 or 2,4, population biomass approaches the coexistence steady state E210.38,1.66, and for 1,0.2 or 5,0.3, predator population goes extinct; that is, Nt,Pt tends to E14.44,0 for long time interval. Hence, in this case, bistability between E1 and E21 occurs. In Figure 5(e), no interior steady state exists; hence, the predator-free steady state E14.44,0 is globally asymptotically stable. Different trajectories for initial population biomass 0.1,3.5, 0.5,1.5, 3,4, and 4,2 are drawn in Figure 5. Parameter values are described in Figure 5.Figure 5
Phase portrait for different representative values ofα and k taken from each region of Figure 8. In each figure, all the steady states, specified by small red-coloured filled circles, are the intersection of the prey nullclines (red-coloured curves) and the predator nullclines (green-coloured curves). Trajectories (blue-coloured curves) are drawn for different initial conditions, and all the initial conditions are marked by small blue-coloured circles. The black dotted orbit represents an unstable limit cycle that occurs due to subcritical Hopf bifurcation. Other parameter values are the same as in Figure 6. (a) E0(0,0) and E1(4.44,0) are unstable; E2(0.13,1.11) is locally asymptotically stable. (b) E0(0,0) and E1(4.44,0) are unstable; limit cycle occurs around E2(0.17,2.19). (c) E0(0,0) is unstable; E1(4.44,0) and E210.38,1.66 are locally asymptotically stable; E222.62,0.53 is unstable. (d) E0(0,0) is unstable; E1(4.44,0) is locally asymptotically stable; limit cycle occurs around E210.27,2.41; E223.73,0.66 is unstable. (e) E0(0,0) is unstable; E1(4.44,0) is globally asymptotically stable. No coexistence is possible.
(a)(b)(c)(d)(e)


### 3.2.1. Local Bifurcations: System (39)

Now, we observe that system (39) undergoes various local bifurcations as the parameters k and α cross through critical thresholds.(a) Transcritical Bifurcation. We prove that the axial equilibrium point E1 exchanges its stability with coexistence steady state E2 as the bifurcating parameter α crosses the threshold value α=α^; that is, a transcritical bifurcation occurs between E1 and E2. For α=α^, the Jacobian matrix of system (39) at E1 is evaluated as follows:(43)JE1;α=α^=−r−d−rk1−ηr−da−mθ00.Clearly one of the eigenvalues is zero. The eigenvectors ofJE1;α=α^ and JE1;α=α^T corresponding to the zero eigenvalue are given by V=v11 and W=01 respectively, where v1=−rk1−η/a−m/θr−d. Using the expressions for V and W, we compute the transversality conditions as follows:(44)WTFαE1,α=α^=0,WTDFαE1,α=α^V=−acθ−mh2r−dθa2+r−d2≠0,WTD2FE1,α=α^V,V=−2mθ−mha2−r−d2θa2+r−d2rk1−ηa+mθr−d≠0,provided r−d/a≠1. All the transversality conditions for transcritical bifurcation are satisfied, and hence by Sotomayor’s theorem, we can conclude that system (39) exhibits transcritical bifurcation for α=α^. In numerical demonstration, we call this α^ as αtc. We conclude this result in the following theorem.Theorem 6.
The predator-free equilibrium pointE1 exhibits a transcritical bifurcation when the bifurcation parameter α crosses the threshold value α^=1−m/cθ−mhr−d/a+a/r−d.(b) Saddle-Node Bifurcation. Let r−d/a>1 and choose α as a bifurcation parameter. Then, the equation obtained from (40), that is, N2−c1−αθ−mh/mN+1=0, has a double root N∗ if α=1−2m/cθ−mh=α˜. Geometrically, this means that for the parametric value α=α˜, the curve y=gN=N/1+N2 touches the straight line y=m/c1−αθ−mh at N=1=N∗, and essentially gN will be the maximum at this point. Hence, g′N∗=0. Now, as r−d/a>1, then corresponding this N∗=1, exactly one P∗ exists. Therefore, for α=α˜, exactly one coexistence steady state 1,P∗ exists with multiplicity 2. The Jacobian matrix of system (39) at E21,P∗ is given by(45)JE2;α=α˜=rη+1−η1+kP∗−d−2a−rk1−η1+kP∗2−mθ00.Hence, one eigenvalue ofJE2;α=α˜ is zero. The eigenvectors of JE2;α=α˜ and JE2;α=α˜T corresponding to the zero eigenvalue are computed as V=v11 and W=01, respectively, where v1=rk1−η/1+kP∗+m/θ/rη+1−η/1+kP∗−d−2a, provided the denominator is not zero. For zero denominator, Bogdanov–Takens bifurcation occurs. Two transversality conditions for saddle-node bifurcation are as follows:(46)WTFαE2,α=α˜=−cθ−mh2P∗2θ≠0,WTD2FE2,α=α˜V,V=−mθ−mhθv12≠0.Therefore, by Sotomayor’s theorem, we can conclude that system (39) undergoes a saddle-node bifurcation for α=α˜. In numerical investigation, we call this α˜ as αsn.Theorem 7.
System (39) undergoes a saddle-node bifurcation as the bifurcation parameter α crosses threshold value α˜=1−2m/cθ−mh.(c) Hopf Bifurcation. We know that the coexistence steady state E2 (if exists) lost its stability when the trace of the corresponding Jacobian matrix JE2 becomes positive. Also, note that for multiple coexistence states, only E21 may lose its stability as E22 is always unstable saddle. We show that this instability occurs via Hopf bifurcation. Let the characteristic equation of JE2 for a variable parameter α is given by λ2−Tαλ+Dα=0, where Tα=TrJE2 and Dα=detJE2. For the suitable choice of other parameter values, we may find a critical value of α, say α∗ from the expression of Tα, such that Tα∗=0 and Dα∗>0. Then, the characteristic equation has purely imaginary eigenvalues, which are λ1=iDα∗ and λ2=−iDα∗. The transversality condition for Hopf bifurcation d/dαReλ|α∗≠0 is also satisfied. Therefore, system (39) undergoes a Hopf bifurcation, when α=α∗.For better visualization of the bifurcation analysis, we execute some numerical study and demonstrate how the fear levelk and the strength of habitat complexity α control the system dynamics. In this regard, we consider the following set of parametric values hypothetically: r=0.6,η=0.4,d=0.2,a=0.09,c=0.5,h=0.007,θ=0.4,m=0.01. Figures 6(a)–6(f) illustrate the qualitative change in system dynamics for the variable parameter α in presence of different levels of fear. In Figures 6(a) and 6(b), we have plotted the change of population biomass and the corresponding bifurcation point for the bifurcation parameter α and the lower fear level, fixed at k=0.1. We can observe that for sufficiently lower strength of habitat complexity, coexistence steady state is locally asymptotically stable. Further, if we increase the value of α to α∗=αh=0.6293, system (39) exhibits limit cycle behaviour and undergoes a Hopf bifurcation around the coexistence state. Also, we calculate the first Lyapunov coefficient numerically which is positive in sign. Hence, the Hopf bifurcation is subcritical, and the bifurcating limit cycle is unstable. Potentially this subcritical case is much more dramatic and dangerous in system biology. Once the strength of habitat complexity crosses the threshold value αh, the trajectories of this system must bounce to the nearest attractor (in this case, it must be either a limit cycle with high amplitude or predator-free steady state) which is far away from the coexistence steady state. Whenever this attractor is a limit cycle, it is impossible to get back the stable coexistence steady state by reversing the value of α (see Figure 7(a)). For either case, that is, when the attractor is the predator-free steady state, we have to reverse the α value adequately to get back the coexistence state (see Figure 7(b)). In other words, population biomass is either trapped in a high amplitude limit cycle through the subcritical Hopf bifurcation or exhibits a looping behaviour. This phenomenon is known as hysteresis. Now at α^=αtc=0.7665, the system undergoes a transcritical bifurcation and the predator-free equilibrium becomes locally asymptotically stable. A saddle-node bifurcation also occurs at a higher degree of habitat complexity, α˜=αsn=0.9, through which the coexistence steady states disappear. Multiple coexistences (E21 and E22) and bistability phenomenon occur in between α^ and α˜, which is specified by the grey-shaded region. Figure 6(a) demonstrates that prey equilibrium biomass increases continuously due to increasing value of habitat complexity, whereas Figure 6(b) shows that predator equilibrium biomass initially increases, but for higher value of habitat complexity, predator species goes extinct.Figure 6
Bifurcation diagram of codimension-1 for the bifurcation parameterα. (a, b) The change of population biomass for the lower fear level, k=0.1, in which case a subcritical Hopf bifurcation at αh=0.6293, a saddle-node bifurcation at αsn=0.9, and a transcritical bifurcation at αtc=0.7665 occur. (c, d) The change of population biomass for the fear level k=0.8, in which case another Hopf bifurcation at αh2=0.89 occurs and that Hopf bifurcation is supercritical. But for larger fear level k=1.5, both Hopf bifurcations disappear in (e, f). Multiple coexistence occurs in the grey-shaded region. Blue-coloured curve and red-coloured curve indicate the prey/predator biomass of E1=N1,P1=r−d/a,0 and E2=N∗,P∗ (or E21N1∗,P1∗ and E22N2∗,P2∗), respectively. Solid lines and dashed lines refer to the stable and unstable behaviour of the system respectively. Dashed-dot curves indicate the max and min of prey/predator population. Other parameter values are fixed at r=0.6,η=0.4,d=0.2,a=0.09,c=0.5,h=0.007,θ=0.4,m=0.01. (a) Change of prey biomass for the bifurcation parameter α and fixed k = 0.1. (b) Change of predator biomass for the bifurcation parameter α and fixed k = 0.1. (c) Change of prey biomass for the bifurcation parameter α and fixed k = 0.8. (d) Change of predator biomass for the bifurcation parameter α and fixed k = 0.8. (e) Change of prey biomass for the bifurcation parameter α and fixed k = 1.5. (f) Change of predator biomass for the parameter α and fixed k = 1.5.
(a)(b)(c)(d)(e)(f)Figure 7
Evidence of hysteresis phenomenon. Using the black arrows, we show how the predator biomass is either trapped in a high amplitude limit cycle through the subcritical Hopf bifurcation (a) or exhibits a hysteresis loop (b). In (a), we takek=0.8, and in (b), k=1.5. Other parameter values are taken from Figure 6. Grey-shaded regions refer to the existence of multiple coexistence equilibrium points. Diagrams for prey populations are qualitatively similar.
(a)(b)To observe the impact of fear in population dynamics, we analyze the model numerically for slightly increasing value ofk. In Figures 6(c) and 6(d), we have plotted the bifurcation diagram and change of population biomass for the bifurcation parameter α and a fixed fear level at k=0.8. In this case, the subcritical Hopf bifurcation occurs at αh1=0.84, which is larger than the previous one αh=0.6293. Also, another Hopf bifurcation occurs at αh2=0.89 through which the coexistence steady state E21 becomes locally asymptotically stable. Numerically, it can be shown that the first Lyapunov coefficient is negative, and hence this Hopf bifurcation is supercritical. For higher fear level, k=1.5, both Hopf bifurcations disappear (see Figures 6(e) and 6(f)) and only transcritical and saddle-node bifurcations remained unchanged. Hence, a higher level of fear promotes the coexistence in environment provided the habitat complexity is sufficiently low. Note that α^=αtc=1−m/cθ−mhr−d/a+a/r−d and α˜=αsn=1−2m/cθ−mh are free from fear parameter, and, therefore, bistability phenomenon must occur between one coexistence steady state and the predator-free steady state. Also, for sufficiently higher degree of habitat complexity, coexistence in environment fails to exist and the system settles into the predator-free steady state.(d) Bogdanov–Takens Bifurcation. For a suitable set of parameter values, system (39) exhibits a saddle-node bifurcation, as well as Hopf bifurcation. For example, we have already proved that for r−d/a>1, two coexistence steady states E21 and E22 appear through saddle-node bifurcation below the threshold value of α˜=1−2m/cθ−mh. Here, E22 is always unstable saddle, and the stability of E21 depends on the sign of TrJE21. Therefore, we can find another critical value of α=α∗ at which TrJE21=0 and E21 becomes unstable through Hopf bifurcation. Geometrically, Bogdanov–Takens (BT) bifurcation is a point in two-parametric bifurcation plane, where the saddle-node bifurcation curve meets the Hopf bifurcation curve, and at that steady state, the corresponding Jacobian matrix has a zero eigenvalue with algebraic multiplicity two.(e) Bautin Bifurcation. The coexistence equilibrium point E21, if it exists, may change its stability through a supercritical or a subcritical Hopf bifurcation. For a supercritical case, a stable limit cycle occurs through Hopf bifurcation and the first Lyapunov coefficient l1>0, whereas for subcritical case, an unstable limit cycle appears and the first Lyapunov coefficient, l1<0. Now, if l1=0, system (39) exhibits a Bautin bifurcation or generalized Hopf (GH) bifurcation in two-parametric bifurcation plane. This bifurcation point separates the subcritical and supercritical Hopf bifurcations, and for neighbouring parametric values, the system exhibits one stable and one unstable limit cycles, which collide and vanish through a saddle-node bifurcation of periodic orbits.In Figure8, we have plotted the stability regions of various steady states emerging from model (39) in α−k parametric plane fixing all other parameter values at r=0.6,η=0.4,d=0.2,a=0.09,c=0.5,h=0.007,θ=0.4,m=0.01 and inspect how the system dynamics depends on the fear level, as well as the habitat complexity. The figure includes a transcritical bifurcation curve (magenta-coloured), a saddle-node bifurcation curve (blue-coloured), and a Hopf bifurcation curve (black-coloured). Both transcritical curve, α=α^=1−m/cθ−mhr−d/a+a/r−d, and saddle-node curve, α=α˜=1−2m/cθ−mh, are free from fear parameter k. Hence, the system always experiences these two local bifurcations for bifurcating parameter α without any restriction on k. One can observe that these bifurcation curves divide the whole parametric region into five small subregions R1−R5. Particularly in R1, unique coexistence steady state occurs which is locally asymptotically stable. That coexistence state eventually turns into an unstable state via Hopf bifurcation, if the parametric values α,k enter R2 region. Numerically, one can prove that this Hopf bifurcation is subcritical. Now, if the degree of habitat complexity crosses the transcritical bifurcation curve and lies between α^ and α˜, multiple coexistence steady state occurs, of which one is unstable saddle and the other one either is locally asymptotically stable (in R3) or exhibits stable/unstable limit cycle (in R4). Further, in these parametric regions, the predator-free steady state is also locally asymptotically stable. Hence, bistability phenomenon occurs in both R3 and R4. Further, for higher value of α, coexistence states disappear through a saddle-node bifurcation and all the trajectories settle into the predator-free steady state (in R5). At αGH,kGH=0.87,0.9, two Hopf bifurcations—one is subcritical and the other one is supercritical—collide and disappear through saddle-node bifurcation of periodic orbits, for which a generalized Hopf bifurcation (GH) or Bautin bifurcation occurs. Also, a Bogdanov–Takens (BT) bifurcation occurs at αBT,kBT=0.9,0.44, where the Hopf bifurcation curve touches the saddle-node bifurcation curve. Five phase portraits for five representative values of α and k, taken from each region R1−R5, are plotted in Figures 5(a)–5(e) and display the typical behaviour of the solution trajectories of system (39) in that region.Figure 8
Stability region of different equilibrium points and two-parametric bifurcation diagram of the model system (39) in α−k parametric plane. One transcritical bifurcation curve (magenta-coloured line), one saddle-node bifurcation curve (blue-coloured line), and one Hopf bifurcation curve (black-coloured curve) split up the whole region into five small subregions R1−R5. In R1, exactly one coexistence steady state E2 occurs which is locally asymptotically stable. In R2, that E2 becomes unstable through subcritical Hopf bifurcation. Multiple coexistence (E21 and E22) occurs whenever α and k together belong to R3 and R4. In R3, E21 is locally asymptotically stable and E22 is unstable saddle. Further, E21 becomes unstable in R4 through Hopf bifurcation. Also, in these two regions, predator-free equilibrium point E1 is locally asymptotically stable. Hence, bistablity phenomenon occurs in R3 and R4. In R5, no coexistence occurs and E1 is globally asymptotically stable. A two-parametric bifurcation Bogdanov–Takens (BT) and generalized Hopf (GH) appear in this diagram. Other parameter values are the same as in Figure 6.


### 3.2.2. Impact of Fear and Habitat Complexity in Population Biomass

The prey biomass of the coexistence steady stateE2N∗,P∗ is independent of fear parameter k, and so we only discuss here the effects of fear upon the predator biomass. Differentiating (41) with respect to k, we obtain(47)dP∗dk=−rθ1−ηN∗P∗m1+kP∗2+rkθ1−ηN∗<0,which indicates that predator biomass decreases due to increased level of fear; that means due to the scare of being caught by predator, the interaction among predator and prey reduces significantly and consequently predator biomass decreases because of starvation.Now, to obtain the variation of species biomass, we differentiateN∗ and P∗ with respect to α. Here,(48)dN∗dα=mc1−α2θ−mh1−N∗21+N∗22.Hence, ifN∗<1, the prey population increases as the degree of habitat complexity increases, but for N∗>1, prey population decreases. In this particular model, if N∗>1, gN∗ lies in the decreasing branch of gN and hence the corresponding coexistence steady state is unstable saddle, therefore, from stability view point, prey population increases for growing habitat complexity. Now,(49)dP∗dα=mP∗/rθN∗2−a/rm/rθN∗+1−ηk/1+kP∗2dN∗dα.Therefore, variation of predator biomass depends on the sign ofmP∗/rθN∗2−a/r. From Figures 6(a)–6(f), one can observe the variation of species biomass with respect to α. Prey population increases as α increases and attains to its carrying capacity at α=αtc. On the other hand, predator biomass initially increases, but after a certain value of α, predator biomass goes extinct.After some detailed analysis of this particular type of model, we can now recognize the principal characteristics which can be generalized in model (7), where gN is a general dome-shaped catch rate. Let, gN attain its maximum at Nmax and gNmax=M. The number of coexistence steady states depends on the different restrictions of α, which are described in Table 2, with α˜=1−m/cMθ−mh and α^=1−m/cθ−mhgr−d/a. Whenever exactly one coexistence steady state E2 occurs, stability of that E2 depends on the sign of TrJE2. The system experiences a Hopf bifurcation around E2, if TrJE2 becomes zero for some threshold value of α=α∗ or k=k∗. For multiple cases, one coexistence state is either stable or becomes unstable via Hopf bifurcation and the other one is unstable saddle. The system also undergoes a transcritical bifurcation at α=α^. If r−d/a>Nmax, a saddle-node bifurcation occurs at α=α˜. Therefore, a bistability and hysteresis phenomenon appear in between α^ and α˜ for any dome-shaped catch rate.


## 3.2.1. Local Bifurcations: System (39)

Now, we observe that system (39) undergoes various local bifurcations as the parameters k and α cross through critical thresholds.(a) Transcritical Bifurcation. We prove that the axial equilibrium point E1 exchanges its stability with coexistence steady state E2 as the bifurcating parameter α crosses the threshold value α=α^; that is, a transcritical bifurcation occurs between E1 and E2. For α=α^, the Jacobian matrix of system (39) at E1 is evaluated as follows:(43)JE1;α=α^=−r−d−rk1−ηr−da−mθ00.Clearly one of the eigenvalues is zero. The eigenvectors ofJE1;α=α^ and JE1;α=α^T corresponding to the zero eigenvalue are given by V=v11 and W=01 respectively, where v1=−rk1−η/a−m/θr−d. Using the expressions for V and W, we compute the transversality conditions as follows:(44)WTFαE1,α=α^=0,WTDFαE1,α=α^V=−acθ−mh2r−dθa2+r−d2≠0,WTD2FE1,α=α^V,V=−2mθ−mha2−r−d2θa2+r−d2rk1−ηa+mθr−d≠0,provided r−d/a≠1. All the transversality conditions for transcritical bifurcation are satisfied, and hence by Sotomayor’s theorem, we can conclude that system (39) exhibits transcritical bifurcation for α=α^. In numerical demonstration, we call this α^ as αtc. We conclude this result in the following theorem.Theorem 6.
The predator-free equilibrium pointE1 exhibits a transcritical bifurcation when the bifurcation parameter α crosses the threshold value α^=1−m/cθ−mhr−d/a+a/r−d.(b) Saddle-Node Bifurcation. Let r−d/a>1 and choose α as a bifurcation parameter. Then, the equation obtained from (40), that is, N2−c1−αθ−mh/mN+1=0, has a double root N∗ if α=1−2m/cθ−mh=α˜. Geometrically, this means that for the parametric value α=α˜, the curve y=gN=N/1+N2 touches the straight line y=m/c1−αθ−mh at N=1=N∗, and essentially gN will be the maximum at this point. Hence, g′N∗=0. Now, as r−d/a>1, then corresponding this N∗=1, exactly one P∗ exists. Therefore, for α=α˜, exactly one coexistence steady state 1,P∗ exists with multiplicity 2. The Jacobian matrix of system (39) at E21,P∗ is given by(45)JE2;α=α˜=rη+1−η1+kP∗−d−2a−rk1−η1+kP∗2−mθ00.Hence, one eigenvalue ofJE2;α=α˜ is zero. The eigenvectors of JE2;α=α˜ and JE2;α=α˜T corresponding to the zero eigenvalue are computed as V=v11 and W=01, respectively, where v1=rk1−η/1+kP∗+m/θ/rη+1−η/1+kP∗−d−2a, provided the denominator is not zero. For zero denominator, Bogdanov–Takens bifurcation occurs. Two transversality conditions for saddle-node bifurcation are as follows:(46)WTFαE2,α=α˜=−cθ−mh2P∗2θ≠0,WTD2FE2,α=α˜V,V=−mθ−mhθv12≠0.Therefore, by Sotomayor’s theorem, we can conclude that system (39) undergoes a saddle-node bifurcation for α=α˜. In numerical investigation, we call this α˜ as αsn.Theorem 7.
System (39) undergoes a saddle-node bifurcation as the bifurcation parameter α crosses threshold value α˜=1−2m/cθ−mh.(c) Hopf Bifurcation. We know that the coexistence steady state E2 (if exists) lost its stability when the trace of the corresponding Jacobian matrix JE2 becomes positive. Also, note that for multiple coexistence states, only E21 may lose its stability as E22 is always unstable saddle. We show that this instability occurs via Hopf bifurcation. Let the characteristic equation of JE2 for a variable parameter α is given by λ2−Tαλ+Dα=0, where Tα=TrJE2 and Dα=detJE2. For the suitable choice of other parameter values, we may find a critical value of α, say α∗ from the expression of Tα, such that Tα∗=0 and Dα∗>0. Then, the characteristic equation has purely imaginary eigenvalues, which are λ1=iDα∗ and λ2=−iDα∗. The transversality condition for Hopf bifurcation d/dαReλ|α∗≠0 is also satisfied. Therefore, system (39) undergoes a Hopf bifurcation, when α=α∗.For better visualization of the bifurcation analysis, we execute some numerical study and demonstrate how the fear levelk and the strength of habitat complexity α control the system dynamics. In this regard, we consider the following set of parametric values hypothetically: r=0.6,η=0.4,d=0.2,a=0.09,c=0.5,h=0.007,θ=0.4,m=0.01. Figures 6(a)–6(f) illustrate the qualitative change in system dynamics for the variable parameter α in presence of different levels of fear. In Figures 6(a) and 6(b), we have plotted the change of population biomass and the corresponding bifurcation point for the bifurcation parameter α and the lower fear level, fixed at k=0.1. We can observe that for sufficiently lower strength of habitat complexity, coexistence steady state is locally asymptotically stable. Further, if we increase the value of α to α∗=αh=0.6293, system (39) exhibits limit cycle behaviour and undergoes a Hopf bifurcation around the coexistence state. Also, we calculate the first Lyapunov coefficient numerically which is positive in sign. Hence, the Hopf bifurcation is subcritical, and the bifurcating limit cycle is unstable. Potentially this subcritical case is much more dramatic and dangerous in system biology. Once the strength of habitat complexity crosses the threshold value αh, the trajectories of this system must bounce to the nearest attractor (in this case, it must be either a limit cycle with high amplitude or predator-free steady state) which is far away from the coexistence steady state. Whenever this attractor is a limit cycle, it is impossible to get back the stable coexistence steady state by reversing the value of α (see Figure 7(a)). For either case, that is, when the attractor is the predator-free steady state, we have to reverse the α value adequately to get back the coexistence state (see Figure 7(b)). In other words, population biomass is either trapped in a high amplitude limit cycle through the subcritical Hopf bifurcation or exhibits a looping behaviour. This phenomenon is known as hysteresis. Now at α^=αtc=0.7665, the system undergoes a transcritical bifurcation and the predator-free equilibrium becomes locally asymptotically stable. A saddle-node bifurcation also occurs at a higher degree of habitat complexity, α˜=αsn=0.9, through which the coexistence steady states disappear. Multiple coexistences (E21 and E22) and bistability phenomenon occur in between α^ and α˜, which is specified by the grey-shaded region. Figure 6(a) demonstrates that prey equilibrium biomass increases continuously due to increasing value of habitat complexity, whereas Figure 6(b) shows that predator equilibrium biomass initially increases, but for higher value of habitat complexity, predator species goes extinct.Figure 6
Bifurcation diagram of codimension-1 for the bifurcation parameterα. (a, b) The change of population biomass for the lower fear level, k=0.1, in which case a subcritical Hopf bifurcation at αh=0.6293, a saddle-node bifurcation at αsn=0.9, and a transcritical bifurcation at αtc=0.7665 occur. (c, d) The change of population biomass for the fear level k=0.8, in which case another Hopf bifurcation at αh2=0.89 occurs and that Hopf bifurcation is supercritical. But for larger fear level k=1.5, both Hopf bifurcations disappear in (e, f). Multiple coexistence occurs in the grey-shaded region. Blue-coloured curve and red-coloured curve indicate the prey/predator biomass of E1=N1,P1=r−d/a,0 and E2=N∗,P∗ (or E21N1∗,P1∗ and E22N2∗,P2∗), respectively. Solid lines and dashed lines refer to the stable and unstable behaviour of the system respectively. Dashed-dot curves indicate the max and min of prey/predator population. Other parameter values are fixed at r=0.6,η=0.4,d=0.2,a=0.09,c=0.5,h=0.007,θ=0.4,m=0.01. (a) Change of prey biomass for the bifurcation parameter α and fixed k = 0.1. (b) Change of predator biomass for the bifurcation parameter α and fixed k = 0.1. (c) Change of prey biomass for the bifurcation parameter α and fixed k = 0.8. (d) Change of predator biomass for the bifurcation parameter α and fixed k = 0.8. (e) Change of prey biomass for the bifurcation parameter α and fixed k = 1.5. (f) Change of predator biomass for the parameter α and fixed k = 1.5.
(a)(b)(c)(d)(e)(f)Figure 7
Evidence of hysteresis phenomenon. Using the black arrows, we show how the predator biomass is either trapped in a high amplitude limit cycle through the subcritical Hopf bifurcation (a) or exhibits a hysteresis loop (b). In (a), we takek=0.8, and in (b), k=1.5. Other parameter values are taken from Figure 6. Grey-shaded regions refer to the existence of multiple coexistence equilibrium points. Diagrams for prey populations are qualitatively similar.
(a)(b)To observe the impact of fear in population dynamics, we analyze the model numerically for slightly increasing value ofk. In Figures 6(c) and 6(d), we have plotted the bifurcation diagram and change of population biomass for the bifurcation parameter α and a fixed fear level at k=0.8. In this case, the subcritical Hopf bifurcation occurs at αh1=0.84, which is larger than the previous one αh=0.6293. Also, another Hopf bifurcation occurs at αh2=0.89 through which the coexistence steady state E21 becomes locally asymptotically stable. Numerically, it can be shown that the first Lyapunov coefficient is negative, and hence this Hopf bifurcation is supercritical. For higher fear level, k=1.5, both Hopf bifurcations disappear (see Figures 6(e) and 6(f)) and only transcritical and saddle-node bifurcations remained unchanged. Hence, a higher level of fear promotes the coexistence in environment provided the habitat complexity is sufficiently low. Note that α^=αtc=1−m/cθ−mhr−d/a+a/r−d and α˜=αsn=1−2m/cθ−mh are free from fear parameter, and, therefore, bistability phenomenon must occur between one coexistence steady state and the predator-free steady state. Also, for sufficiently higher degree of habitat complexity, coexistence in environment fails to exist and the system settles into the predator-free steady state.(d) Bogdanov–Takens Bifurcation. For a suitable set of parameter values, system (39) exhibits a saddle-node bifurcation, as well as Hopf bifurcation. For example, we have already proved that for r−d/a>1, two coexistence steady states E21 and E22 appear through saddle-node bifurcation below the threshold value of α˜=1−2m/cθ−mh. Here, E22 is always unstable saddle, and the stability of E21 depends on the sign of TrJE21. Therefore, we can find another critical value of α=α∗ at which TrJE21=0 and E21 becomes unstable through Hopf bifurcation. Geometrically, Bogdanov–Takens (BT) bifurcation is a point in two-parametric bifurcation plane, where the saddle-node bifurcation curve meets the Hopf bifurcation curve, and at that steady state, the corresponding Jacobian matrix has a zero eigenvalue with algebraic multiplicity two.(e) Bautin Bifurcation. The coexistence equilibrium point E21, if it exists, may change its stability through a supercritical or a subcritical Hopf bifurcation. For a supercritical case, a stable limit cycle occurs through Hopf bifurcation and the first Lyapunov coefficient l1>0, whereas for subcritical case, an unstable limit cycle appears and the first Lyapunov coefficient, l1<0. Now, if l1=0, system (39) exhibits a Bautin bifurcation or generalized Hopf (GH) bifurcation in two-parametric bifurcation plane. This bifurcation point separates the subcritical and supercritical Hopf bifurcations, and for neighbouring parametric values, the system exhibits one stable and one unstable limit cycles, which collide and vanish through a saddle-node bifurcation of periodic orbits.In Figure8, we have plotted the stability regions of various steady states emerging from model (39) in α−k parametric plane fixing all other parameter values at r=0.6,η=0.4,d=0.2,a=0.09,c=0.5,h=0.007,θ=0.4,m=0.01 and inspect how the system dynamics depends on the fear level, as well as the habitat complexity. The figure includes a transcritical bifurcation curve (magenta-coloured), a saddle-node bifurcation curve (blue-coloured), and a Hopf bifurcation curve (black-coloured). Both transcritical curve, α=α^=1−m/cθ−mhr−d/a+a/r−d, and saddle-node curve, α=α˜=1−2m/cθ−mh, are free from fear parameter k. Hence, the system always experiences these two local bifurcations for bifurcating parameter α without any restriction on k. One can observe that these bifurcation curves divide the whole parametric region into five small subregions R1−R5. Particularly in R1, unique coexistence steady state occurs which is locally asymptotically stable. That coexistence state eventually turns into an unstable state via Hopf bifurcation, if the parametric values α,k enter R2 region. Numerically, one can prove that this Hopf bifurcation is subcritical. Now, if the degree of habitat complexity crosses the transcritical bifurcation curve and lies between α^ and α˜, multiple coexistence steady state occurs, of which one is unstable saddle and the other one either is locally asymptotically stable (in R3) or exhibits stable/unstable limit cycle (in R4). Further, in these parametric regions, the predator-free steady state is also locally asymptotically stable. Hence, bistability phenomenon occurs in both R3 and R4. Further, for higher value of α, coexistence states disappear through a saddle-node bifurcation and all the trajectories settle into the predator-free steady state (in R5). At αGH,kGH=0.87,0.9, two Hopf bifurcations—one is subcritical and the other one is supercritical—collide and disappear through saddle-node bifurcation of periodic orbits, for which a generalized Hopf bifurcation (GH) or Bautin bifurcation occurs. Also, a Bogdanov–Takens (BT) bifurcation occurs at αBT,kBT=0.9,0.44, where the Hopf bifurcation curve touches the saddle-node bifurcation curve. Five phase portraits for five representative values of α and k, taken from each region R1−R5, are plotted in Figures 5(a)–5(e) and display the typical behaviour of the solution trajectories of system (39) in that region.Figure 8
Stability region of different equilibrium points and two-parametric bifurcation diagram of the model system (39) in α−k parametric plane. One transcritical bifurcation curve (magenta-coloured line), one saddle-node bifurcation curve (blue-coloured line), and one Hopf bifurcation curve (black-coloured curve) split up the whole region into five small subregions R1−R5. In R1, exactly one coexistence steady state E2 occurs which is locally asymptotically stable. In R2, that E2 becomes unstable through subcritical Hopf bifurcation. Multiple coexistence (E21 and E22) occurs whenever α and k together belong to R3 and R4. In R3, E21 is locally asymptotically stable and E22 is unstable saddle. Further, E21 becomes unstable in R4 through Hopf bifurcation. Also, in these two regions, predator-free equilibrium point E1 is locally asymptotically stable. Hence, bistablity phenomenon occurs in R3 and R4. In R5, no coexistence occurs and E1 is globally asymptotically stable. A two-parametric bifurcation Bogdanov–Takens (BT) and generalized Hopf (GH) appear in this diagram. Other parameter values are the same as in Figure 6.


## 3.2.2. Impact of Fear and Habitat Complexity in Population Biomass

The prey biomass of the coexistence steady stateE2N∗,P∗ is independent of fear parameter k, and so we only discuss here the effects of fear upon the predator biomass. Differentiating (41) with respect to k, we obtain(47)dP∗dk=−rθ1−ηN∗P∗m1+kP∗2+rkθ1−ηN∗<0,which indicates that predator biomass decreases due to increased level of fear; that means due to the scare of being caught by predator, the interaction among predator and prey reduces significantly and consequently predator biomass decreases because of starvation.Now, to obtain the variation of species biomass, we differentiateN∗ and P∗ with respect to α. Here,(48)dN∗dα=mc1−α2θ−mh1−N∗21+N∗22.Hence, ifN∗<1, the prey population increases as the degree of habitat complexity increases, but for N∗>1, prey population decreases. In this particular model, if N∗>1, gN∗ lies in the decreasing branch of gN and hence the corresponding coexistence steady state is unstable saddle, therefore, from stability view point, prey population increases for growing habitat complexity. Now,(49)dP∗dα=mP∗/rθN∗2−a/rm/rθN∗+1−ηk/1+kP∗2dN∗dα.Therefore, variation of predator biomass depends on the sign ofmP∗/rθN∗2−a/r. From Figures 6(a)–6(f), one can observe the variation of species biomass with respect to α. Prey population increases as α increases and attains to its carrying capacity at α=αtc. On the other hand, predator biomass initially increases, but after a certain value of α, predator biomass goes extinct.After some detailed analysis of this particular type of model, we can now recognize the principal characteristics which can be generalized in model (7), where gN is a general dome-shaped catch rate. Let, gN attain its maximum at Nmax and gNmax=M. The number of coexistence steady states depends on the different restrictions of α, which are described in Table 2, with α˜=1−m/cMθ−mh and α^=1−m/cθ−mhgr−d/a. Whenever exactly one coexistence steady state E2 occurs, stability of that E2 depends on the sign of TrJE2. The system experiences a Hopf bifurcation around E2, if TrJE2 becomes zero for some threshold value of α=α∗ or k=k∗. For multiple cases, one coexistence state is either stable or becomes unstable via Hopf bifurcation and the other one is unstable saddle. The system also undergoes a transcritical bifurcation at α=α^. If r−d/a>Nmax, a saddle-node bifurcation occurs at α=α˜. Therefore, a bistability and hysteresis phenomenon appear in between α^ and α˜ for any dome-shaped catch rate.


## 4. Discussion

Constructing a predator-prey interaction in modelling framework and identifying the crucial components that determine the ecological complexity are the main ecological research objectives. In particular, functional response is one of the fundamental components that determines the dynamics of a predator-prey interaction. In this study, we have derived a functional response based on timescale separation arguments, which includes the catch rate,gN, and depending upon the monotonic property of gN, various functional responses emerge. In particular, if gN increases strictly, a saturating functional response occurs; however, if gN increases initially and after some certain prey biomass due to group defense gN decreases, the resulting functional response will be in a dome-shaped form. Various experimental evidences suggest that increasing structural complexity of habitat reduces the foraging efficiency of predators, and consequently, the attack rate or encounter rate decreases significantly. It is also reasonably true that the presence of predator may affects the physiological behaviour of prey species to such an extent that it can be more effective in reducing prey biomass than direct predation. Incorporating these biological aspects, we have investigated the predator-prey model in terms of two different shaped catch rates and recognized how the choice of two different shaped functional responses affect the qualitative behaviour of a system.Our proposed model is biologically feasible in the sense that any positive solution initiating in the first quadrant remains positively invariant and uniformly bounded for any timet≥0. We have established the uniform persistent condition for the general model that ensures the existence of the positive equilibrium point of the system. The local asymptotic stability of different equilibrium points has been investigated. By taking a particular form of fear function, we have executed some extensive and comparative analyses of the model dynamics with two different shaped functional responses.For a model with saturating functional response, only one interior equilibrium point can exist, but for dome-shaped functional response, up to two coexistence equilibria can occur and, in general, the complexity of the dynamical behaviour is increased. The geometrical explanations for existence of steady state(s) for both saturating and dome-shaped functional response are discussed in Figure1. We have implemented bifurcation analysis for the degree of habitat complexity α revealing that higher degree of habitat complexity increases the extinction possibility of the predator and the carrying capacity of the prey will be the only stable attractor, no matter how sensitive the prey population is to potential dangers the predator population. Also, interestingly, we cannot obtain a stable predator-free steady state by purely increasing the level of fear.For the saturating case, lowering the degree of habitat complexity gives rise to stable coexistence and further deceasing may produce oscillatory behaviour (see Figures2(a)–2(d)). However, this oscillations can be abolished if the level of fear increases sufficiently (see Figures 2(e)–2(f)). But the whole dynamics will collapse, if the coexistence steady state fails to exit, and this dynamical discrepancy may appear for two different but similar saturating-shaped functional responses. Therefore, one can observe the structural sensitivity in this type of model systems. Although the structural sensitivity fails to predict the accurate dynamical behaviour, from our investigations, we can anticipate the overall possibilities.On the other hand, for dome-shaped functional response, coexistence steady state will be locally asymptotically stable, even if the degree of habitat complexity decreases to zero. In that case, the coexistence state is surrounded by an unstable limit cycle which is emerged due to subcritical Hopf bifurcation (see Figures6(a)–6(d)). Once the system undergoes a subcritical Hopf bifurcation for the bifurcation parameter α, the trajectories of that system must bounce to the nearest attractor which is far away from the coexistence steady state, and this dramatic behaviour can be overcome through a larger alteration of α, or sometimes it is impossible to overcome (hysteresis phenomena). For the intermediate values of α, the system exhibits multiple coexistence steady states and a bistability occurs between predator-free steady state and a coexistence steady state (see Figures 6(a)–6(d)). Now, if we adequately increase the level of fear, oscillatory behaviour will be disappeared which agree with dynamics for saturating case except the bistable phenomenon (see Figures 6(e) and 6(f)).Thus the present study describes different possible dynamical behaviour of the predator-prey system with a comparison between different shaped function responses. Hence, our study will be helpful to understand the interplay between predator and prey with more ecological impacts.



---

*Source: 6427864-2021-09-11.xml*