# Modeling and Analysis of Predator-Prey Model with Fear Effect in Prey and Hunting Cooperation among Predators and Harvesting

Authors: Basaznew Belew; Dawit Melese
Journal: Journal of Applied Mathematics
 (2022)

---

## Abstract

In this paper, we present and analyze predator-prey system where prey population is linearly harvested and affected by fear and the prey population has grown logistically in the absence of predators. The predator population follows hunting cooperation, and it predates the prey population in the Holling type II functional responses. Based on those assumptions, a two-dimensional mathematical model is derived. The positivity, boundedness, and extinction of both prey and predator populations of the solution of the system are discussed. The existence, stability (local and global), and the Hopf bifurcation analysis of the biologically feasible equilibrium points are investigated. The aim of this research is to explore the effect of fear on the prey population and hunting cooperation on the predator population, and both prey and predator populations are harvested linearly and taken as control parameters of the model. If the values ofc1>1, then both prey and predator populations are extinct and also fear parameter has a stabilizing effect on system 4. From the numerical simulation, it was found that the fear effect, hunting cooperation, prey harvesting, and predator harvesting change the dynamics of system 4.


---

## Body

## 1. Introduction

One of the most challenging tasks in ecological model is understanding ecosystem dynamics. Population ecology is concerned with the developing theory and insight about the persistence, structure, and dynamics of biological communities. The predator-prey dynamics are one of the most important topics in theoretical ecology due to their universal existence and importance [1]. One must carefully construct a biologically meaningful and mathematically tractable population model [2]. The most important element in a population model is a predator-prey model which describes the number (density) of prey consumed per predator per unit time for given quantities (density) of prey and predator communities [3].The harvesting of the prey species or the predator species or both is another important issue from an ecological and economic point of view. It generally has a strong impact on the dynamics of biological resources. Sometimes, harvesting of predators is used to control the predator population so as to prevent the extinction of the prey species. Concerning the conservation for the long-term benefits of humanity, there is a wide range of interest in the use of bioeconomic modeling to gain insight in the scientific management of renewable resources like fisheries and forestries [4].The emerging view is that indirect effect (fear effect) on prey population is even more powerful than direct effect [5–7]. Due to fear of predation, prey animals can change its grazing zone to a safer place and sacrifice their highest intake rate areas [8–11]. At the last time, many researchers have proposed and studied fear effect of predator-prey model. Here, due to fear of predation, growth rate of prey reduces, i.e., increase fear effect implies decrease growth rate of prey [7, 11]. Recently, however, [10] conducted a manipulation on song sparrows during an entire breeding season to determine whether perceived predation risk could affect reproduction even in the absence of direct killing. For example, [7] have investigated the modeling of the fear effect in predator-prey model; it is considered as mortality rate of predator and the Holling type I functional responses and showed the local and global stability of the feasible equilibria and the existence of the Hopf bifurcation.The role of mutual beneficial interactions (e.g., cooperation) among social animals is rapidly growing in the research field in population dynamics and conversion biology cooperative hunting among vertebrates (lions [12], African wild dogs [13], chimpanzees [14], wolves [15], birds [16], primates [17], etc.). Different animals adopt different strategies for making a group; for example, wolves have the strategies to get close to the prey and catch it for the acquisition of large and faster prey [18]. In general, animals of the same species compute with each other for resources, but because some specific mechanism is at work, they help each other for mutual benefit and show cooperative behavior. The hunting cooperation in predator-prey model has been investigated by many researchers [10, 11, 19–24].In recent years, more attention has been paid to the study of the dynamics of hunting cooperation and fear in a predator-prey interaction [24–26]. For example, [26] design the mathematical model of predator-prey model with fear and hunting cooperation with the Lessie-Gower functional responses, and they investigated the local stability and bifurcation analysis. However, most of the works focus on the nonharvesting of hunting cooperation and fear effect of a predator-prey system. Thus, the main aim of this paper is to study the impact of fear effect in prey, hunting cooperation among predator, and prey and predator harvesting of a predator-prey system.The organization of this paper is as follows: Section2 is where the mathematical model formulation is discussed; Section 3 deals with the analysis, existence, and boundedness of the model, extinction criteria, the local and global stability analysis of the biologically feasible equilibrium points, and bifurcation analysis of the system; Section 4 deals with the numerical simulation results that are discussed; Section 5, lastly, deals with the conclusion part.The novelty of this paper is the consideration of fear effect in prey, hunting cooperation among predator, and both prey and predator harvesting with nonlinear rate with the Holling type II functional responses.


## 2. Mathematical Model Formulation

In this paper, the following assumptions are taken and a prey population is denoted byNT and predator population is denoted by PT.
(1)
The prey species grows logistically in the absence of predator with carrying capacityK and internsic growth rate r(2)
The predator species feed their prey in the Holling type II functional response(3)
The birth rate of prey population is reduced due to fear induced by predators(4)
The predator population uses hunting cooperation among themselves for capturing the prey population(5)
The predator natural death rate is assumed to be constantm1(6)
The preyNT and predator PT are assumed to be harvested by an external force according to the linear type of functional harvesting (proportionate harvesting)Based on the above assumptions and parameter definitions in Table1, the model equation is given by the following system of nonlinear differential equations.
(1)dNdT=rN1+τP1−NK−λ+aPNP1+htλ+aPN−H1Q1N,(2)dPdT=eλ+aPNP1+htλ+aPN−m1P−H2Q2P.N0≥0 and P0≥0. All these parameters in equation (1) are assumed to be positive. To reduce the number of parameters, we nondimensionalized the model system using the following scaling:
(3)n=NK,p=λPr,t=Tr.Table 1
Parameter definition.
ParametersBiological meaningrIntrinsic growth rate of preykEnvironmental carrying capacitym1Natural death rate of predatoreConversion factorτFear parameteraPredator cooperation in huntingλAttack rate per predator and preyhtHandling timeQ1The combined external effort devoted to nonselective harvesting of the prey by the external harvester (not by predator)Q2The combined external effort devoted to nonselective harvesting of the predator by the external harvesterH1The prey catchability coefficientH2Predator catchability coefficientThe nondimensionalized form of equation (1) can be written as follows:
(4)dndt=n1−n1+βp−1+αppn1+H1+αpn−c1n=N11,dpdt=δ1+αpnp1+H1+αpn−mp−c2p=P11,where
(5)β=τrλ,α=arλ2,H=htλk,c1=H1Q1r,δ=eλkr,m=m1r,c2=H2Q2r,with n0≥0,p0≥0.


## 3. The Model Analysis

### 3.1. Existence and Positive in Variance

Fort>0, let Y=nt,ptT and FY=F1Y,F2YT, where F1Y=n1−n/1+βp−1+αppn/1+H1+αpn1+αppn/1+H1+αpn−c1n,F2Y=δ1+αpnp/1+H1+αpn−mp−c2p. Then, system (4) can be written as dY/dt=FY, where F:C+⟶R+2 with Y0=Y0∈R+2. Here, Fi∈C∞R for i=1,2. Thus, the function F is locally Lipschitzian and completely continuous on R+2. Therefore, the solution of system (4) with nonnegative initial condition exists and is unique. Moreover, it can be shown that these solutions exist for all t>0 and stay nonnegative. Hence, the region R+2 is an invariant domain of system (4).


### 3.2. Boundedness

In theoretical ecology model, boundedness of a system implies that the system is biologically well-behaved. The following theorem ensures the boundedness of system (4).Theorem 1.
The solutionst=nt,pt of system (4) which initiate in interior of R2 is uniformly bounded.
Proof. The first equation of system (4) that represents the prey equation is bounded through, dn/dt≤n1−n. So by comparison theorem, we obtain nt≤1+Ce−t, for all t≤0, where C=1−n0/n0. Thus, supt⟶∞nt=1, which implies that nt≤1 for sufficiently large t. Since x is bounded, without any loss of generality it follows that the solutions of y is bounded. To show this reality, we assure this mathematically as follows: we define a positive definite function as follows:
(6)Wt=nt+1δpt.
From definition ofWt, Wt is differential in some maximal interval 0,tb. Thus, differentiating Wt with respect to time and substituting the reaction terms of system (4) yields
(7)dWtdt=n1−n1+βp−1+αpnp1+H1+αpn−c1n+1+αpnp1+H1+αpn−1δmp+c2p.=n1−n1+βp−c1n−1δmp+c2p.≤n1−n−c1n−pδm+c2≤14−c1−pδm+c2,sinceMaxn1−n=14.≤14−γwtwhereγ=minc1,m+c2.
Here, we have(8)dwtdt+γwt≤14.
Applying the theory of differential inequality, we obtain(9)0≤wt≤14γ1−e−γt+w0e−γt.
Fort⟶∞, we have 0<wt<1/4γ. Hence, all the solution of system (4) starting in R+2 for any ε>0 are contained in the region R=n,p∈R2,n+p/δ≤1/4γ+ε.A system is said to be dissipative if all population initiating inR2+=fn,p:n>0;p>0 are uniformly limited by their environment. Thus, the following theorem establishes the dissipativeness of system (4).Theorem 2.
System (4) is dissipative.
Proof. System (4) is uniformly bounded. Therefore, system (4) is dissipative.


### 3.3. Extinction Criteria

Theorem 3.
Ifc1>1, then limt⟶∞nt=0.
Proof. We have(10)dndt=−c1n−np1+pα1+Hn1+pα+1−nn1+pβ.≤n1−c1.
Integrating both sides and it becomes(11)nt≤n0e∫0t1−c1dt.
Therefore,limt⟶∞nt=0, if c1>1.Theorem 4.
Ifδ<Hm+c2, then limt⟶∞pt=0.
Proof. We have(12)dpdt=−c2p−mp+np1+pαδ1+Hn1+pα≤δpH−mp−c2p.
Since1+H1+αpn>H1+αpn, we have
(13)np1+pαδ1+Hn1+pα<np1+pαδHn1+pα.
Integrating both sides and it becomes(14)pt≤p0e∫0tδ−Hm+c2/Hdt.
Therefore, we havelimt⟶∞pt=0, if δ<Hm+c2.


### 3.4. Model Analysis

#### 3.4.1. Equilibrium Points

System (4) has the following biologically feasible equilibrium points.
(1)
Trivial equilibrium pointE00,0(2)
Axial equilibrium pointE11−c1,0 which exists for c1<1(3)
The positive interior equilibrium pointE2n2,p2 exists for
(15)δ>Hc2+m,1−c1Hc2+m−δ>1,α>21−c1,c1Hmβδ3<H2c2+m22α+β−Hc2+m2β+α4+c1βδ+β+α2+c1βδ2,where
(16)n2=c2+m1+pαδ−c2+mH.Ifn2 exists, δ>c2+mH, and p2 is the unique positive root of the polynomial equation
(17)a4p4+a3p3+a2p2+a1p+a0=0,with
(18)a0=δ1+H1−c1Hc2+m−δ,a1=H2c2+m2+c1δ2α+β+c1c2Hβ+αδ+Hδc2+mα−2+c1α,a2=H2c2+m22α+β−Hc2+m2β+α4+c1βδ+β+α2+c1βδ2−c1Hmβδ3,a3=αα+2β−Hc2+m+δ2,a4=βHc2+mα−αδ2.One can see thatp2 exists for using the Discartes rule of sign; change equation (17) that has a unique positive solution; the following conditions are satisfied.
(19)a4>0,a3>0,a0>0for1−c1Hc2+m−δ>1,a1>0forα>21−c1,a2<0forc1Hmβδ3<H2c2+m22α+β−Hc2+m2β+α4+c1βδ+β+α2+c1βδ2.


### 3.5. Stability Analysis

In order to investigate the local stability property of system (4), we shall calculate the Jacobian matrix at each equilibrium. The Jacobian matrix of system (4) at any arbitrary equilibrium En,p is given as
(20)JE=a11a12a21a22,where
(21)a11=−c1−p1+pα1+Hn+npα2+1−2n1+pβ,a12=−1+2pα+Hn1+pα21+Hn+npα2+−1+nβ1+pβ2,a21=p1+pαδ1+Hn+npα2,a22=−c2−m+n1+2pα+Hn1+pα2δ1+Hn+Hnpα2.Theorem 5.
The trivial equilibrium pointE00,0 is stable if c1>1.
Proof. The eigenvalues of the Jacobian matrixJ of system (4) at E00,0 are 1−c1,−c2−m. All the eigenvalues are negative if and only if c1>1. Therefore, system (4) at the equilibrium point E00,0 is locally asymptotically stable.Theorem 6.
The equilibrium pointE11−c1,0 of system (4) is locally asymptotically stable if c1<1 and c2+m>−1+c1δ/−1+−1+c1H.
Proof. The eigenvalues of the Jacobian matrix of system (4) at the axial equilibrium point E11−c1,0 are −1+c1,−c2−m+−1+c1δ/−1+−1+c1H. Therefore, the equilibrium point E11−c1,0 is locally asymptotically stable if c1<1 and c2+m>−1+c1δ/−1+−1+c1H.Theorem 7.
The coexistence equilibrium pointE2 is locally asymptotically stable if
(22)c1+c2+m>1−2n1+pβ+−p1+pα+n1+2pα+Hn1+pα2δ1+Hn+Hnpα2,a11a22−a12a21>0,where
(23)a11=−c1−p21+p2α1+Hn2+n2p2α2+1−2n21+p2β,(24)a12=n−1+2pα+Hn1+pα21+Hn+npα2+−1+nβ1+pβ2,(25)a21=p21+p2αδ1+Hn2+n2p2α2,a22=−c2−m+n21+2p2α+Hn21+p2α2δ1+Hn2+Hn2p2α2.
Proof. The Jacobian matrixJ2 of system (4) around the interior equilibrium point E2 is given by Jij,i,j=1,2 where aij are given by (23). Let λi,i=1,2 be the roots of the characteristics equation of J2 which is given by
(26)λ2+A1λ+A2=0,where
(27)A1=−a11−a22,A2=a11a22−a12a21.
By the Routh-Hurwitz criterion,E2 is locally asymptotically stable if and only if A1>0 and A2>0. Then, we have A1>0, for c1+c2+m>1−2n/1+pβ+−p1+pα+n1+2pα+Hn1+pα2δ/1+Hn+Hnpα2. Therefore, the coexistence equilibrium point E2 is locally asymptotically stable if the conditions in (22) hold.


### 3.6. Global Stability Analysis of the Coexistence Equilibrium Point

Theorem 8.
System (4) does not admit any periodic solution around the positive equilibrium point E2 provided
(28)−1+H1−c1<0,H1−c2+m+δ<0.
Proof. Letnt,pt be solution of system (4). Define the Dulac function
(29)Q=1+H1+αpnnp.
Then, we have(30)∂N11Q∂n+∂P11Q∂p=−αpc1+c2+m+2n−1H−βpαpc1+c2+m+c1Hpβp+1.+H−c1−2n+1+αδpβp+1−1pβp+1.=−1+H1−c1−H2n+c1αp+2nαp+c1βpp1+pβ.+αpδ+H1−c2+m+αβp2δ−c1+c2+mHp1+pβ.<0, provided (28) holds.
Thus, by the Dulac criterion, system (4) does not have a nontrivial periodic solution. Therefore, the local asymptotically stability of system (4) around E2 and the boundedness of the solutions of system (4) ensure its globally asymptotically stability around the coexistence equilibrium point.


### 3.7. Hopf Bifurcation

Characteristic equation of system (4) at E2n2,p2 is given by
(31)λ2+A1αλ+A2α=0,where
(32)A1=−a11−a22,A2=a11a22−a12a21,with aij,i,j=1,2 defined as in (23).Theorem 9.
System (31) possesses a Hopf bifurcation around E2 when α passes through α∗ provided that A1α=0.
Proof. Atα=α∗, the characteristic of equation (4) can be written as λ2+A2α=0, which has roots λ1=iA2 and λ2=−iA2. Thus, there exist a pair of purely complex eigenvalues. Also A1 and A2 are the smooth functions of α. So for α neighborhood of α∗, the roots have the form λ1=q1α+iq2α and λ2=q2α−iq2α, where qiα are real functions for i=1,2. Next, we verify the transversality condition.
(33)ddαReλiα=α∗≠0,i=1,2.
Now, substitutingq1α+iq2α=λK into (31), derivating it along α, and comparing real and imaginary parts, respectively, we get
(34)2q1q1′−2q2q2′+A1′q1+A1q1′+A2′=0,2q1q2′+2q2q1′+A1′q2+A1q2′i=0.
Simplifying it and it becomes(35)q1′X1−q2′X2+X3=0,q2′X1+q1′X2+X4=0.where
(36)X1=2q1+A1≠2q1,X2=2q2,X3=A1′q1+A2,X4=A1′q2.
Thus, we have(37)q1′=−X3X1+X4X2X12+X22.
Atα=α∗, we get
(38)ddαReλiα=α∗=−X3X1+X4X2X12+X22≠0.
Hence, the transversality condition holds. Therefore, the Hopf Bifurcation occurs atα=α∗.


## 3.1. Existence and Positive in Variance

Fort>0, let Y=nt,ptT and FY=F1Y,F2YT, where F1Y=n1−n/1+βp−1+αppn/1+H1+αpn1+αppn/1+H1+αpn−c1n,F2Y=δ1+αpnp/1+H1+αpn−mp−c2p. Then, system (4) can be written as dY/dt=FY, where F:C+⟶R+2 with Y0=Y0∈R+2. Here, Fi∈C∞R for i=1,2. Thus, the function F is locally Lipschitzian and completely continuous on R+2. Therefore, the solution of system (4) with nonnegative initial condition exists and is unique. Moreover, it can be shown that these solutions exist for all t>0 and stay nonnegative. Hence, the region R+2 is an invariant domain of system (4).


## 3.2. Boundedness

In theoretical ecology model, boundedness of a system implies that the system is biologically well-behaved. The following theorem ensures the boundedness of system (4).Theorem 1.
The solutionst=nt,pt of system (4) which initiate in interior of R2 is uniformly bounded.
Proof. The first equation of system (4) that represents the prey equation is bounded through, dn/dt≤n1−n. So by comparison theorem, we obtain nt≤1+Ce−t, for all t≤0, where C=1−n0/n0. Thus, supt⟶∞nt=1, which implies that nt≤1 for sufficiently large t. Since x is bounded, without any loss of generality it follows that the solutions of y is bounded. To show this reality, we assure this mathematically as follows: we define a positive definite function as follows:
(6)Wt=nt+1δpt.
From definition ofWt, Wt is differential in some maximal interval 0,tb. Thus, differentiating Wt with respect to time and substituting the reaction terms of system (4) yields
(7)dWtdt=n1−n1+βp−1+αpnp1+H1+αpn−c1n+1+αpnp1+H1+αpn−1δmp+c2p.=n1−n1+βp−c1n−1δmp+c2p.≤n1−n−c1n−pδm+c2≤14−c1−pδm+c2,sinceMaxn1−n=14.≤14−γwtwhereγ=minc1,m+c2.
Here, we have(8)dwtdt+γwt≤14.
Applying the theory of differential inequality, we obtain(9)0≤wt≤14γ1−e−γt+w0e−γt.
Fort⟶∞, we have 0<wt<1/4γ. Hence, all the solution of system (4) starting in R+2 for any ε>0 are contained in the region R=n,p∈R2,n+p/δ≤1/4γ+ε.A system is said to be dissipative if all population initiating inR2+=fn,p:n>0;p>0 are uniformly limited by their environment. Thus, the following theorem establishes the dissipativeness of system (4).Theorem 2.
System (4) is dissipative.
Proof. System (4) is uniformly bounded. Therefore, system (4) is dissipative.


## 3.3. Extinction Criteria

Theorem 3.
Ifc1>1, then limt⟶∞nt=0.
Proof. We have(10)dndt=−c1n−np1+pα1+Hn1+pα+1−nn1+pβ.≤n1−c1.
Integrating both sides and it becomes(11)nt≤n0e∫0t1−c1dt.
Therefore,limt⟶∞nt=0, if c1>1.Theorem 4.
Ifδ<Hm+c2, then limt⟶∞pt=0.
Proof. We have(12)dpdt=−c2p−mp+np1+pαδ1+Hn1+pα≤δpH−mp−c2p.
Since1+H1+αpn>H1+αpn, we have
(13)np1+pαδ1+Hn1+pα<np1+pαδHn1+pα.
Integrating both sides and it becomes(14)pt≤p0e∫0tδ−Hm+c2/Hdt.
Therefore, we havelimt⟶∞pt=0, if δ<Hm+c2.


## 3.4. Model Analysis

### 3.4.1. Equilibrium Points

System (4) has the following biologically feasible equilibrium points.
(1)
Trivial equilibrium pointE00,0(2)
Axial equilibrium pointE11−c1,0 which exists for c1<1(3)
The positive interior equilibrium pointE2n2,p2 exists for
(15)δ>Hc2+m,1−c1Hc2+m−δ>1,α>21−c1,c1Hmβδ3<H2c2+m22α+β−Hc2+m2β+α4+c1βδ+β+α2+c1βδ2,where
(16)n2=c2+m1+pαδ−c2+mH.Ifn2 exists, δ>c2+mH, and p2 is the unique positive root of the polynomial equation
(17)a4p4+a3p3+a2p2+a1p+a0=0,with
(18)a0=δ1+H1−c1Hc2+m−δ,a1=H2c2+m2+c1δ2α+β+c1c2Hβ+αδ+Hδc2+mα−2+c1α,a2=H2c2+m22α+β−Hc2+m2β+α4+c1βδ+β+α2+c1βδ2−c1Hmβδ3,a3=αα+2β−Hc2+m+δ2,a4=βHc2+mα−αδ2.One can see thatp2 exists for using the Discartes rule of sign; change equation (17) that has a unique positive solution; the following conditions are satisfied.
(19)a4>0,a3>0,a0>0for1−c1Hc2+m−δ>1,a1>0forα>21−c1,a2<0forc1Hmβδ3<H2c2+m22α+β−Hc2+m2β+α4+c1βδ+β+α2+c1βδ2.


## 3.4.1. Equilibrium Points

System (4) has the following biologically feasible equilibrium points.
(1)
Trivial equilibrium pointE00,0(2)
Axial equilibrium pointE11−c1,0 which exists for c1<1(3)
The positive interior equilibrium pointE2n2,p2 exists for
(15)δ>Hc2+m,1−c1Hc2+m−δ>1,α>21−c1,c1Hmβδ3<H2c2+m22α+β−Hc2+m2β+α4+c1βδ+β+α2+c1βδ2,where
(16)n2=c2+m1+pαδ−c2+mH.Ifn2 exists, δ>c2+mH, and p2 is the unique positive root of the polynomial equation
(17)a4p4+a3p3+a2p2+a1p+a0=0,with
(18)a0=δ1+H1−c1Hc2+m−δ,a1=H2c2+m2+c1δ2α+β+c1c2Hβ+αδ+Hδc2+mα−2+c1α,a2=H2c2+m22α+β−Hc2+m2β+α4+c1βδ+β+α2+c1βδ2−c1Hmβδ3,a3=αα+2β−Hc2+m+δ2,a4=βHc2+mα−αδ2.One can see thatp2 exists for using the Discartes rule of sign; change equation (17) that has a unique positive solution; the following conditions are satisfied.
(19)a4>0,a3>0,a0>0for1−c1Hc2+m−δ>1,a1>0forα>21−c1,a2<0forc1Hmβδ3<H2c2+m22α+β−Hc2+m2β+α4+c1βδ+β+α2+c1βδ2.


## 3.5. Stability Analysis

In order to investigate the local stability property of system (4), we shall calculate the Jacobian matrix at each equilibrium. The Jacobian matrix of system (4) at any arbitrary equilibrium En,p is given as
(20)JE=a11a12a21a22,where
(21)a11=−c1−p1+pα1+Hn+npα2+1−2n1+pβ,a12=−1+2pα+Hn1+pα21+Hn+npα2+−1+nβ1+pβ2,a21=p1+pαδ1+Hn+npα2,a22=−c2−m+n1+2pα+Hn1+pα2δ1+Hn+Hnpα2.Theorem 5.
The trivial equilibrium pointE00,0 is stable if c1>1.
Proof. The eigenvalues of the Jacobian matrixJ of system (4) at E00,0 are 1−c1,−c2−m. All the eigenvalues are negative if and only if c1>1. Therefore, system (4) at the equilibrium point E00,0 is locally asymptotically stable.Theorem 6.
The equilibrium pointE11−c1,0 of system (4) is locally asymptotically stable if c1<1 and c2+m>−1+c1δ/−1+−1+c1H.
Proof. The eigenvalues of the Jacobian matrix of system (4) at the axial equilibrium point E11−c1,0 are −1+c1,−c2−m+−1+c1δ/−1+−1+c1H. Therefore, the equilibrium point E11−c1,0 is locally asymptotically stable if c1<1 and c2+m>−1+c1δ/−1+−1+c1H.Theorem 7.
The coexistence equilibrium pointE2 is locally asymptotically stable if
(22)c1+c2+m>1−2n1+pβ+−p1+pα+n1+2pα+Hn1+pα2δ1+Hn+Hnpα2,a11a22−a12a21>0,where
(23)a11=−c1−p21+p2α1+Hn2+n2p2α2+1−2n21+p2β,(24)a12=n−1+2pα+Hn1+pα21+Hn+npα2+−1+nβ1+pβ2,(25)a21=p21+p2αδ1+Hn2+n2p2α2,a22=−c2−m+n21+2p2α+Hn21+p2α2δ1+Hn2+Hn2p2α2.
Proof. The Jacobian matrixJ2 of system (4) around the interior equilibrium point E2 is given by Jij,i,j=1,2 where aij are given by (23). Let λi,i=1,2 be the roots of the characteristics equation of J2 which is given by
(26)λ2+A1λ+A2=0,where
(27)A1=−a11−a22,A2=a11a22−a12a21.
By the Routh-Hurwitz criterion,E2 is locally asymptotically stable if and only if A1>0 and A2>0. Then, we have A1>0, for c1+c2+m>1−2n/1+pβ+−p1+pα+n1+2pα+Hn1+pα2δ/1+Hn+Hnpα2. Therefore, the coexistence equilibrium point E2 is locally asymptotically stable if the conditions in (22) hold.


## 3.6. Global Stability Analysis of the Coexistence Equilibrium Point

Theorem 8.
System (4) does not admit any periodic solution around the positive equilibrium point E2 provided
(28)−1+H1−c1<0,H1−c2+m+δ<0.
Proof. Letnt,pt be solution of system (4). Define the Dulac function
(29)Q=1+H1+αpnnp.
Then, we have(30)∂N11Q∂n+∂P11Q∂p=−αpc1+c2+m+2n−1H−βpαpc1+c2+m+c1Hpβp+1.+H−c1−2n+1+αδpβp+1−1pβp+1.=−1+H1−c1−H2n+c1αp+2nαp+c1βpp1+pβ.+αpδ+H1−c2+m+αβp2δ−c1+c2+mHp1+pβ.<0, provided (28) holds.
Thus, by the Dulac criterion, system (4) does not have a nontrivial periodic solution. Therefore, the local asymptotically stability of system (4) around E2 and the boundedness of the solutions of system (4) ensure its globally asymptotically stability around the coexistence equilibrium point.


## 3.7. Hopf Bifurcation

Characteristic equation of system (4) at E2n2,p2 is given by
(31)λ2+A1αλ+A2α=0,where
(32)A1=−a11−a22,A2=a11a22−a12a21,with aij,i,j=1,2 defined as in (23).Theorem 9.
System (31) possesses a Hopf bifurcation around E2 when α passes through α∗ provided that A1α=0.
Proof. Atα=α∗, the characteristic of equation (4) can be written as λ2+A2α=0, which has roots λ1=iA2 and λ2=−iA2. Thus, there exist a pair of purely complex eigenvalues. Also A1 and A2 are the smooth functions of α. So for α neighborhood of α∗, the roots have the form λ1=q1α+iq2α and λ2=q2α−iq2α, where qiα are real functions for i=1,2. Next, we verify the transversality condition.
(33)ddαReλiα=α∗≠0,i=1,2.
Now, substitutingq1α+iq2α=λK into (31), derivating it along α, and comparing real and imaginary parts, respectively, we get
(34)2q1q1′−2q2q2′+A1′q1+A1q1′+A2′=0,2q1q2′+2q2q1′+A1′q2+A1q2′i=0.
Simplifying it and it becomes(35)q1′X1−q2′X2+X3=0,q2′X1+q1′X2+X4=0.where
(36)X1=2q1+A1≠2q1,X2=2q2,X3=A1′q1+A2,X4=A1′q2.
Thus, we have(37)q1′=−X3X1+X4X2X12+X22.
Atα=α∗, we get
(38)ddαReλiα=α∗=−X3X1+X4X2X12+X22≠0.
Hence, the transversality condition holds. Therefore, the Hopf Bifurcation occurs atα=α∗.


## 4. Numerical Simulation

In the present study, the rate of prey harvestingc1, the rate of predator harvesting c2, the fear factor coefficient α, and hunting cooperation coefficient β are the key parameters which will be taken as control parameters. The numerical simulation is carried out using the MATLAB software package for the set of parameter values given in Table 2.Table 2
Parameter value of system (4).
ParameterNumeric valueSourceH0.9Assumedδ0.5Assumedm0.1Assumedc1Variable—c2Variables—αVariable—βVariable—


### 4.1. Effect of Varying the Prey Harvesting Ratec1

Let us fix the variables in Table2 as c2=0.05,α=0.25,β=0.2.It is observed that system (4) approaches to the trivial equilibrium point for c1≥1, to the predator-free equilibrium point (axial equilibrium point) for 0.589<c1<1, and to the coexistence equilibrium point for c1<0.5891.The illustration of the coexistence equilibrium point E2 exists and is locally asymptotically stable if c1<0.5891, when c1 lies in the range of 0.589<c1<1; then, the coexistence equilibrium point annihilates with the predator-free equilibrium point (axial equilibrium point) E1, and E1 changes its stability and becomes locally asymptotically stable. Further increasing c1 lies in the range of c1≥1 which leads to the predator-free equilibrium point that annihilates with the trivial equilibrium point E0, and E0 changes its stability and becomes locally asymptotically stable.Figures1(a) and 1(b), respectively, show that the time series and parametric plot of the trajectory of system (4), i.e., it approaches locally asymptotically to the coexistence equilibrium point E2 for c1=0.3. Figure 2(a) shows that the trajectories of system (4) approach locally asymptotically to the predator-free equilibrium point (axial equilibrium point) E1 for c1=0.7, i.e., one can see that the prey species can coexist whereas the predator species goes to extinction as the value of 0.589<c1<1. Figure 2(b) shows that the trajectories of system (4) approach locally asymptotically to the trivial equilibrium point E0 stability of E0 for c1=1, i.e., one can see that both species are extinct as the value of c1≥1. Figure 3(a) shows that an increase in the prey harvesting rate c1 leads to an increase in the prey population, whereas Figure 3(b) illustrates that as the rate of harvesting c1 increases, the density of the predator population decreases.Figure 1
(a) Time series solution of system (4) (b) and parametric plot of the equilibrium point of system (4) approach asymptotically to the coexistence equilibrium point E2 with parametric values as in Table 2 except c2=0.05,c1=0.3,andα=0.25,β=0.2.
(a)(b)Figure 2
(a) The time series solution of system (4) approaches locally asymptotically to the equilibrium point E1 with parametric value as in Table 2 except c2=0.05,c1=0.7,α=0.25,andβ=0.2, (b) The time series solution of system (4) approaches to the equilibrium point E0 except c2=0.05,c1=1,α=0.25,andβ=0.2.
(a)(b)Figure 3
(a) The density of prey population for the parameter values as in Table2 except c2=0.05,α=0.25,andβ=0.2 and c1=0.02,0.1,0.3,0.5, (b) The density of predator population for the parameter values as in Table 2 except c2=0.05,α=0.25,andβ=0.2 and c1=0.02,0.1,0.3,0.5.
(a)(b)


### 4.2. Effect of Predator Harvesting Ratec2

For the parametric values as in Table2 with c1=0.12,α=0.25,andβ=0.2, the predator-free equilibrium point (axial equilibrium point) E1 and the coexistence equilibrium point E2 exists for 0.145536<c2<1 and c2<0.145536, respectively. The illustration of the coexistence equilibrium point E2 exists and is locally asymptotically stable in c2<0.145536. And c2 lies between in the range 0.145536<c2<1; the coexistence equilibrium point annihilates with the predator-free equilibrium point (axial equilibrium point) E1, and E1 changes its stability and becomes locally asymptotically stable.Figures4(a) and 4(b) show the stability of the time series and parametric plot of system (4), respectively, i.e., it approaches to the coexistence equilibrium point E2 for c2=0.1. Figure 5(a) shows the stability of system (4) approaches to locally asymptotically to the predator-free equilibrium point (axial equilibrium point) E1 for c2=0.2. Furthermore, one can see that the prey species coexist whereas the predator species goes to extinction whenever 0.145536<c2<1. From Figures 6(a) and 6(b), it can be observed that an increase in the harvesting rate of predator leads to increase in the prey population density and a decrease predator population density, respectively.Figure 4
(a) Time series solution and (b) parametric plot of system (4) about the equilibrium point E2 with parametric value as in Table 2 except c2=0.1,c1=0.12,α=0.25,andβ=0.2.
(a)(b)Figure 5
(a) Time series solution of system (4) approaches asymptotically to the equilibrium point E1 with parametric value as in Table 2 except c2=0.2,,α=0.25,andβ=0.2.Figure 6
(a) The density of prey population for the parameter values as in Table2 except c1=0.12,α=0.25,andβ=0.2 and c2=0.01,0.07,0.1,0.14 and (b) the density of predator population for the parameter values as in Table 2 except c1=0.12,α=0.25,andβ=0.2 and c2=0.01,0.07,0.1,0.14.
(a)(b)


### 4.3. Effect of Hunting Cooperationα

The figures show that the population size of both prey and predator oscillate, but they are not periodic. The amplitudes of the oscillation of both populations approach the equilibrium point of the system.Figures7(a) and 7(b) show the stability of the time series and parametric plot of system (4), respectively, i.e., it approaches asymptotically to the coexistence equilibrium point E2 for α=1. Figures 8(a) and 8(b) show the instability of time series and parametric plot of the system, respectively, i.e., it does not approach asymptotically to the equilibrium point E2 for α=2. From Figures 9(a) and 9(b), it can be observed that an increase in the hunting cooperation of predator leads to a decrease in the prey population density and predator population density. It is observed that for the parameter values given in Table 2 with c1=0.12,c2=0.05,andα=0.2, the trajectory of system (4) approaches the coexistence equilibrium point for α<1.81905 and the Hopf bifurcation occurs at α=1.81905. Otherwise, it is unstable.Figure 7
(a) Time series solution and (b) parametric plot of system (4) about the equilibrium point E2 with parametric value as in Table 2 except c2=0.05,c1=0.12,α=1,andβ=0.2.
(a)(b)Figure 8
(a) Time series solution and (b) parametric plot of system (4) about the equilibrium point E2 with parametric value as in Table 2, respectively, except c2=0.05,c1=0.12,α=2,andβ=0.2.
(a)(b)Figure 9
(a) The density of prey population for the parameter values as in Table2 except c1=0.12,α=0.06,0.4,1.4,1.7,β=0.2 and c2=0.05. (b) The density of predator population for the parameter values as in Table 2 except c1=0.12,α=0.06,0.4,1.4,1.7,β=0.2 and c2=.05.
(a)(b)


### 4.4. Effect of Fearβ

Now, to investigate the impact of fear effect on the dynamics on system (4), for the parametric values with c1=0.35,α=0.5,H=0.5,δ=0.5,m=0.1,andc2=0.1, the coexistence equilibrium point E2 exists for 0.159042<β<9.966312; otherwise, the predator-free equilibrium point exists.The coexistence equilibrium pointE2 exists and is locally asymptotically stable in 0.159042<β<9.966312. And β lies within the range 0.159042<β<9.966312; the coexistence equilibrium point annihilates with the predator-free equilibrium point (axial equilibrium point) E1, and E1 changes its stability and becomes locally asymptotically stable.Figures10(a) and 10(b) show the stability of the time series and parametric plot of system (4), i.e., it approaches asymptotically to the coexistence equilibrium point E2 for β=0.95. Figure 11(a) shows that the stability of system (4) approaches asymptotically to the predator-free equilibrium point (axial equilibrium point) E1 for β=10. Furthermore, one can see that the prey species coexist whereas the predator species goes to extinction as the value of β is out of the range 0.159042<β<9.966312. From Figures 12(a) and 12(b), it can be observed that an increase in the fear strength leads to an increase in the prey population density and a decrease in the predator population density.Figure 10
(a) Time series solution and (b) parametric plot of system (4) about the equilibrium point E2 with parametric value as in Table 2 except H=0.5,c2=0.1,c1=0.35,α=0.5,andβ=0.95.
(a)(b)Figure 11
Time series solution of system (4) about the equilibrium point E2 with parametric value as in Table 2 except H=0.5,c2=0.1,c1=0.35,α=0.5,andβ=10.Figure 12
(a) The density of prey population for the parameter values as in Table2 except c1=0.35,α=0.5,andβ=0.7,2,5,7 and c2=0.1, and (b) the density of predator population for the parameter values as in Table 2 except c1=0.35,α=0.5,andβ=0.7,2,5,7 and c2=0.1.
(a)(b)Here, we choose different values ofα; system (2) shows stable focuses and periodic oscillation. If we incorporate with fear strength, a fear phenomenon into the system for a fixed hunting cooperation α>1.81904, the unstable system becomes stable for a fixed values of α=1.9 other than α=0.25. Figures 13(a) and 13(b) show the instability of system (4) about the coexistence equilibrium point E2 whenever β=0 for a fixed value of α=1.9. Figure 14 shows the stability of system (4) about the coexistence equilibrium point E2 whenever β=0.75 for a fixed value of α=1.9. We can observe that fear factor (β) has a stabilizing effect on the dynamics of the system (c.f. Figures 13(a) and 13(b) and Figures 14(a) and 14(b)).Figure 13
(a) Time series solution and (b) parametric plot of system (4) about the equilibrium point E2 with parametric value as in Table 2 except c2=0.05,c1=0.12,α=1.9,andβ=0.
(a)(b)Figure 14
(a) Time series solution and (b) parametric plot of system (4) about the equilibrium point E2 with parametric value as in Table 2 except c2=0.05,c1=0.12,α=1.9,andβ=0.75.
(a)(b)Further, we investigate the global stability of system (4) for c1=0.12,c2=0.05,α=0.5,β=0.2,H=0.9,m=0.1,andδ=0.5 with different initial conditions I1,I2,I3,andI4, which shows the global converges of the trajectories of system (2) to the coexistence equilibrium point E2=0.361,0.549 (c.f. Figure 15). This shows the globally asymptotically stability of the interior equilibrium point E2.Figure 15
Phase plan trajectory of system (4) corresponding to different initial conditions I1=0.569,0.398,I2=0.478,0.367,I3=0.695,0.266,I4=0.713,0.158,andI5=0.297,0.512.On the analysis of system (4), one can see that, for the parametric values as in Table 2 and c1=0.12c2=0.05,andβ=0.2, the endemic equilibrium point E2 is locally asymptotically stable for α<1.81905 and unstable for α>1.81905. The existence of the Hopf bifurcation occurs at α=1.81905 (see Figure 16).Figure 16
Hopf bifurcation diagram of system (4) around the coexistence equilibrium point E2 for the data set in Table 2 except c1=0.12c2=0.05,andβ=0.2.


## 4.1. Effect of Varying the Prey Harvesting Ratec1

Let us fix the variables in Table2 as c2=0.05,α=0.25,β=0.2.It is observed that system (4) approaches to the trivial equilibrium point for c1≥1, to the predator-free equilibrium point (axial equilibrium point) for 0.589<c1<1, and to the coexistence equilibrium point for c1<0.5891.The illustration of the coexistence equilibrium point E2 exists and is locally asymptotically stable if c1<0.5891, when c1 lies in the range of 0.589<c1<1; then, the coexistence equilibrium point annihilates with the predator-free equilibrium point (axial equilibrium point) E1, and E1 changes its stability and becomes locally asymptotically stable. Further increasing c1 lies in the range of c1≥1 which leads to the predator-free equilibrium point that annihilates with the trivial equilibrium point E0, and E0 changes its stability and becomes locally asymptotically stable.Figures1(a) and 1(b), respectively, show that the time series and parametric plot of the trajectory of system (4), i.e., it approaches locally asymptotically to the coexistence equilibrium point E2 for c1=0.3. Figure 2(a) shows that the trajectories of system (4) approach locally asymptotically to the predator-free equilibrium point (axial equilibrium point) E1 for c1=0.7, i.e., one can see that the prey species can coexist whereas the predator species goes to extinction as the value of 0.589<c1<1. Figure 2(b) shows that the trajectories of system (4) approach locally asymptotically to the trivial equilibrium point E0 stability of E0 for c1=1, i.e., one can see that both species are extinct as the value of c1≥1. Figure 3(a) shows that an increase in the prey harvesting rate c1 leads to an increase in the prey population, whereas Figure 3(b) illustrates that as the rate of harvesting c1 increases, the density of the predator population decreases.Figure 1
(a) Time series solution of system (4) (b) and parametric plot of the equilibrium point of system (4) approach asymptotically to the coexistence equilibrium point E2 with parametric values as in Table 2 except c2=0.05,c1=0.3,andα=0.25,β=0.2.
(a)(b)Figure 2
(a) The time series solution of system (4) approaches locally asymptotically to the equilibrium point E1 with parametric value as in Table 2 except c2=0.05,c1=0.7,α=0.25,andβ=0.2, (b) The time series solution of system (4) approaches to the equilibrium point E0 except c2=0.05,c1=1,α=0.25,andβ=0.2.
(a)(b)Figure 3
(a) The density of prey population for the parameter values as in Table2 except c2=0.05,α=0.25,andβ=0.2 and c1=0.02,0.1,0.3,0.5, (b) The density of predator population for the parameter values as in Table 2 except c2=0.05,α=0.25,andβ=0.2 and c1=0.02,0.1,0.3,0.5.
(a)(b)


## 4.2. Effect of Predator Harvesting Ratec2

For the parametric values as in Table2 with c1=0.12,α=0.25,andβ=0.2, the predator-free equilibrium point (axial equilibrium point) E1 and the coexistence equilibrium point E2 exists for 0.145536<c2<1 and c2<0.145536, respectively. The illustration of the coexistence equilibrium point E2 exists and is locally asymptotically stable in c2<0.145536. And c2 lies between in the range 0.145536<c2<1; the coexistence equilibrium point annihilates with the predator-free equilibrium point (axial equilibrium point) E1, and E1 changes its stability and becomes locally asymptotically stable.Figures4(a) and 4(b) show the stability of the time series and parametric plot of system (4), respectively, i.e., it approaches to the coexistence equilibrium point E2 for c2=0.1. Figure 5(a) shows the stability of system (4) approaches to locally asymptotically to the predator-free equilibrium point (axial equilibrium point) E1 for c2=0.2. Furthermore, one can see that the prey species coexist whereas the predator species goes to extinction whenever 0.145536<c2<1. From Figures 6(a) and 6(b), it can be observed that an increase in the harvesting rate of predator leads to increase in the prey population density and a decrease predator population density, respectively.Figure 4
(a) Time series solution and (b) parametric plot of system (4) about the equilibrium point E2 with parametric value as in Table 2 except c2=0.1,c1=0.12,α=0.25,andβ=0.2.
(a)(b)Figure 5
(a) Time series solution of system (4) approaches asymptotically to the equilibrium point E1 with parametric value as in Table 2 except c2=0.2,,α=0.25,andβ=0.2.Figure 6
(a) The density of prey population for the parameter values as in Table2 except c1=0.12,α=0.25,andβ=0.2 and c2=0.01,0.07,0.1,0.14 and (b) the density of predator population for the parameter values as in Table 2 except c1=0.12,α=0.25,andβ=0.2 and c2=0.01,0.07,0.1,0.14.
(a)(b)


## 4.3. Effect of Hunting Cooperationα

The figures show that the population size of both prey and predator oscillate, but they are not periodic. The amplitudes of the oscillation of both populations approach the equilibrium point of the system.Figures7(a) and 7(b) show the stability of the time series and parametric plot of system (4), respectively, i.e., it approaches asymptotically to the coexistence equilibrium point E2 for α=1. Figures 8(a) and 8(b) show the instability of time series and parametric plot of the system, respectively, i.e., it does not approach asymptotically to the equilibrium point E2 for α=2. From Figures 9(a) and 9(b), it can be observed that an increase in the hunting cooperation of predator leads to a decrease in the prey population density and predator population density. It is observed that for the parameter values given in Table 2 with c1=0.12,c2=0.05,andα=0.2, the trajectory of system (4) approaches the coexistence equilibrium point for α<1.81905 and the Hopf bifurcation occurs at α=1.81905. Otherwise, it is unstable.Figure 7
(a) Time series solution and (b) parametric plot of system (4) about the equilibrium point E2 with parametric value as in Table 2 except c2=0.05,c1=0.12,α=1,andβ=0.2.
(a)(b)Figure 8
(a) Time series solution and (b) parametric plot of system (4) about the equilibrium point E2 with parametric value as in Table 2, respectively, except c2=0.05,c1=0.12,α=2,andβ=0.2.
(a)(b)Figure 9
(a) The density of prey population for the parameter values as in Table2 except c1=0.12,α=0.06,0.4,1.4,1.7,β=0.2 and c2=0.05. (b) The density of predator population for the parameter values as in Table 2 except c1=0.12,α=0.06,0.4,1.4,1.7,β=0.2 and c2=.05.
(a)(b)


## 4.4. Effect of Fearβ

Now, to investigate the impact of fear effect on the dynamics on system (4), for the parametric values with c1=0.35,α=0.5,H=0.5,δ=0.5,m=0.1,andc2=0.1, the coexistence equilibrium point E2 exists for 0.159042<β<9.966312; otherwise, the predator-free equilibrium point exists.The coexistence equilibrium pointE2 exists and is locally asymptotically stable in 0.159042<β<9.966312. And β lies within the range 0.159042<β<9.966312; the coexistence equilibrium point annihilates with the predator-free equilibrium point (axial equilibrium point) E1, and E1 changes its stability and becomes locally asymptotically stable.Figures10(a) and 10(b) show the stability of the time series and parametric plot of system (4), i.e., it approaches asymptotically to the coexistence equilibrium point E2 for β=0.95. Figure 11(a) shows that the stability of system (4) approaches asymptotically to the predator-free equilibrium point (axial equilibrium point) E1 for β=10. Furthermore, one can see that the prey species coexist whereas the predator species goes to extinction as the value of β is out of the range 0.159042<β<9.966312. From Figures 12(a) and 12(b), it can be observed that an increase in the fear strength leads to an increase in the prey population density and a decrease in the predator population density.Figure 10
(a) Time series solution and (b) parametric plot of system (4) about the equilibrium point E2 with parametric value as in Table 2 except H=0.5,c2=0.1,c1=0.35,α=0.5,andβ=0.95.
(a)(b)Figure 11
Time series solution of system (4) about the equilibrium point E2 with parametric value as in Table 2 except H=0.5,c2=0.1,c1=0.35,α=0.5,andβ=10.Figure 12
(a) The density of prey population for the parameter values as in Table2 except c1=0.35,α=0.5,andβ=0.7,2,5,7 and c2=0.1, and (b) the density of predator population for the parameter values as in Table 2 except c1=0.35,α=0.5,andβ=0.7,2,5,7 and c2=0.1.
(a)(b)Here, we choose different values ofα; system (2) shows stable focuses and periodic oscillation. If we incorporate with fear strength, a fear phenomenon into the system for a fixed hunting cooperation α>1.81904, the unstable system becomes stable for a fixed values of α=1.9 other than α=0.25. Figures 13(a) and 13(b) show the instability of system (4) about the coexistence equilibrium point E2 whenever β=0 for a fixed value of α=1.9. Figure 14 shows the stability of system (4) about the coexistence equilibrium point E2 whenever β=0.75 for a fixed value of α=1.9. We can observe that fear factor (β) has a stabilizing effect on the dynamics of the system (c.f. Figures 13(a) and 13(b) and Figures 14(a) and 14(b)).Figure 13
(a) Time series solution and (b) parametric plot of system (4) about the equilibrium point E2 with parametric value as in Table 2 except c2=0.05,c1=0.12,α=1.9,andβ=0.
(a)(b)Figure 14
(a) Time series solution and (b) parametric plot of system (4) about the equilibrium point E2 with parametric value as in Table 2 except c2=0.05,c1=0.12,α=1.9,andβ=0.75.
(a)(b)Further, we investigate the global stability of system (4) for c1=0.12,c2=0.05,α=0.5,β=0.2,H=0.9,m=0.1,andδ=0.5 with different initial conditions I1,I2,I3,andI4, which shows the global converges of the trajectories of system (2) to the coexistence equilibrium point E2=0.361,0.549 (c.f. Figure 15). This shows the globally asymptotically stability of the interior equilibrium point E2.Figure 15
Phase plan trajectory of system (4) corresponding to different initial conditions I1=0.569,0.398,I2=0.478,0.367,I3=0.695,0.266,I4=0.713,0.158,andI5=0.297,0.512.On the analysis of system (4), one can see that, for the parametric values as in Table 2 and c1=0.12c2=0.05,andβ=0.2, the endemic equilibrium point E2 is locally asymptotically stable for α<1.81905 and unstable for α>1.81905. The existence of the Hopf bifurcation occurs at α=1.81905 (see Figure 16).Figure 16
Hopf bifurcation diagram of system (4) around the coexistence equilibrium point E2 for the data set in Table 2 except c1=0.12c2=0.05,andβ=0.2.


## 5. Conclusion

In this paper, it can be conclude that the formulated model is mathematically meaningfully valid and biologically well posed by providing boundedness positivity and existence of the solution of the model. Trivial, axial, and coexistence equilibrium points are investigated. Moreover, it is observed that in our model, trivial equilibrium point is always locally asymptotically stable forc1>1. Axial equilibrium point is locally asymptotically stable if and only if the variable satisfies the condition c1<1 and c2+m>−1+c1δ/−1+−1+c1H. The local and global stability condition for the coexistence equilibrium point of system (4) are obtained. The prey population will be extinct if the prey growth rate is less than the prey harvesting rate (see Theorem 3) and the predator population extincts (Theorem 4 holds). A bifurcation analysis of system 4 undergoes a Hopf bifurcation under a certain condition.Numerical simulations verify our analytical findings. The following main results are obtained for the parametric values given in Table2.Increasing prey harvestingc1 lies in the range 0.589<c1<1 which leads to the extinction of predator and it is locally asymptotically stable to the axial equilibrium point (predator-free equilibrium point) (see Figure 2(a)). Again, further increase c1 to the range c1>1results in the extinction of both the prey and predator populations and it is locally asymptotically stable to the trivial equilibrium point (see Figure 2(b)). Finally increase c1 lies in the range c1<0.589 and leads to increase in prey population and decreased in predator population; then, system 4 is locally asymptotically stable to the coexistence equilibrium point (see Figures 3(a) and 3(b)), respectively.Increase predator harvestingc2 lies in the range between 0.145536<c2<1 which leads to the extinction of predator population and it is locally asymptotically stable to the axial equilibrium point (predator-free equilibrium point) (see Figure 5(a)). Again, increase c2 lies in the range c2<0.145536 which increases the prey population and decreases the predator population and it is locally asymptotically stable to the coexistence equilibrium point (see Figures 6(a) and 6(b)).Increase hunting cooperation of predator populationα lies in the range α<1.81905 which leads to decrease both the prey and predator populations and it is locally asymptotically stable to the coexistence equilibrium point (see Figures 9(a) and 9(b)). The coexistence equilibrium point loses its stability, and the Hopf bifurcation occurs at α=1.81905; otherwise, it is unstable.Increase fear effect of prey populationβ lies in the range out of 0.159042<β<9.966312 which leads to the extinction of predator population and it is locally asymptotically stable to the axial equilibrium point (predator-free equilibrium point) (see Figure 11(a)). Again, further β lies in the range 0.159042<β<9.966312 which increases the prey population and decreases the predator population and it is locally asymptotically stable to the coexistence equilibrium point (see Figures 12(a) and 12(b)). Finally, we can observe that fear factor β has a stabilizer effect on the dynamics of system 4 (see Figures 13(a) and 13(b) and Figures 14(a) and 14(b)).



---

*Source: 2776698-2022-12-03.xml*