# Dynamics of a Predator-Prey Population in the Presence of Resource Subsidy under the Influence of Nonlinear Prey Refuge and Fear Effect

Authors: Sudeshna Mondal; G. P. Samanta; Juan J. Nieto
Journal: Complexity
 (2021)

---

## Abstract

In this work, our aim is to investigate the impact of a non-Kolmogorov predator-prey-subsidy model incorporating nonlinear prey refuge and the effect of fear with Holling type II functional response. The model arises from the study of a biological system involving arctic foxes (predator), lemmings (prey), and seal carcasses (subsidy). The positivity and asymptotically uniform boundedness of the solutions of the system have been derived. Analytically, we have studied the criteria for the feasibility and stability of different equilibrium points. In addition, we have derived sufficient conditions for the existence of local bifurcations of codimension 1 (transcritical and Hopf bifurcation). It is also observed that there is some time lag between the time of perceiving predator signals through vocal cues and the reduction of prey’s birth rate. So, we have analyzed the dynamical behaviour of the delayed predator-prey-subsidy model. Numerical computations have been performed using MATLAB to validate all the analytical findings. Numerically, it has been observed that the predator, prey, and subsidy can always exist at a nonzero subsidy input rate. But, at a high subsidy input rate, the prey population cannot persist and the predator population has a huge growth due to the availability of food sources.


---

## Body

## 1. Introduction

In the ecological system, the predator-prey interaction is one of the most significant tools which is comparatively easy to observe in the field. But fear of the predator felt by the prey (indirect effect) also plays a vital role since its effect is stronger than direct predation [1, 2]. The cost of fear can reduce the reproduction rate of prey because it affects the physiological condition of prey population. As a result, the prey species may get a long-term loss. In support of this, it is mentioned that, in the Greater Yellowstone Ecosystem, wolves (Canis lupus) affect the reproductive physiology of elk (Cervus elaphus) [3]. When the prey species recognize the predator signal (chemical/vocal), they spend more time to become keenly watchful to detect danger rather than in foraging. So, the birth rate of the scared prey reduces and adopts some survival mechanisms like starvation [1, 2]. For examples, some birds react to the sound of predator with antipredator defenses [1, 2] and they flee from their nests at the first sign of danger [2]. This antipredator behaviour may affect survival and reproduction of the birds [2]. It has been experimentally investigated that, in the absence of direct killing, the reproduction of the offspring of song sparrows (Melospiza melodia) could be reduced by 40% as a result of impact of feeling fear created by the predator [4]. So, this reduction caused by the antipredator behaviour affects the birth rate and survival of offspring. Thus, the cost of fear (apart from direct predation) should be introduced in a predator-prey interaction. Mathematical formulation of the impact of fear on the two species prey-predator system has been initiated by Wang et al. [5] in 2016 introducing fear factor: fk,y=1/1+ky. It involves a parameter k denoting the level of fear to represent the antipredator behaviour of the prey. Some research works have already been done on the ecological system under the influence of predation fear felt by prey species [6–17]. Moreover, the impact of fear in a two-species predator-prey model with prey refuge was analyzed by many researchers [11, 18, 19].In evolutionary biology, prey refuge is a concept which helps an organism to protect themselves from predation by hiding in an area inaccessible to the predator, for example, in a wolf-ungulate system, ungulates may seek refuge by migrating to areas outside the core territories of wolves. Also, it has many significant roles on the dynamics of predator-prey interactions: prey refuge may decrease the chance of extinction of prey. Researchers have mainly used the dynamic nature of predator-prey model with linear prey refuge (that is,mx amounts of prey are unavailable to the predator, where m∈0,1 is the coefficient of refuge and x is the biomass of prey species) with Volterra response [20–22]. Recently, Mondal and Samanta have studied the dynamics of the predator-prey system with prey refuge dependent on both species (that is, mxy amounts of prey are free from the predator risk, where 0<1−my<1, m is prey refuge coefficient, x is the biomass of prey population, and y is the biomass of predator population) in the presence of additional food (for details, see [23]). In 2020, Mondal and Samanta [11] have also analyzed the dynamics of predator-prey interaction having nonlinear prey refuge function Φx,y=mxy/a+y which is the amount of prey that are free from predation, where a is half saturation constant and y is the biomass of predator population.Many experimental studies suspect that the introduction of resource subsidies may disrupt otherwise stable food web linkages [24–26]. Such concept is significant for resource management purposes. It is learnt that reintroduced wolves in Yellowstone Park switch to bison when their preferred ungulate prey, namely, elk, are rare in the concerned ecosystem [27]. Mathematically, the influence of resource subsidy on the predator-prey model has been initiated in the work of Nevai and Van Gorder [28]. They have discussed how different subsidy input rate may affect the prey and predator population to persist in the ecosystem.Generalist predator can consume more than one food source: either multiple prey population or a combination of prey population and resource subsidy. There are many rich theoretical research on ecological systems involving generalist predator [29–31]. Also, there are a variety of real-life applications for such systems [32–34]. From literature surveys, it has been shown that generalist predator can persist in an ecosystem even if one particular prey species is going towards extinction [29–31].In 2012, Nevai and Van Gorder [28] first extended the Kolmogorov model to a non-Kolmogorov predator-prey-subsidy model. It has been observed that the predator-prey-subsidy model occurs in the arctic foxes (predator), lemmings (prey), and seal carcasses (subsidy). Motivated by the works of Das and Samanta [35], Nevai and Van Gorder [28], and Xu et al. [36], we have analyzed the dynamical behaviour of a mathematical model of non-Kolmogorov form that includes the three components (predator, prey, and subsidy) with the impacts of nonlinear prey refuge function and the fear effect felt by the prey in the presence of the predator. To the best of our knowledge, there does not exist any mathematical model to explore the impact of fear effect incorporating nonlinear prey refuge function in predator-prey-subsidy interaction.The organization of this work is structured as follows: in Section2, a mathematical model has been formulated with the influences of nonlinear prey refuge and fear effect. Section 3 shows that the proposed model is well-behaved. In Section 4, feasibility criteria and stability of all the equilibria of the proposed system (in absence of delay) have been studied. The equilibria can change their stability nature through transcritical and Hopf bifurcation which are also analyzed in this section. Generally, the reduction of prey’s birth rate due to the effect of fear will not be an instantaneous biological process but deviated through some time lag, so the study of time-delay τ is very meaningful to obtain the more realistic dynamics. So, Section 5 deals with the dynamic behaviour of the delayed system for two equilibrium points E3 (subsidy free) and E∗ (interior), respectively. Section 6 provides the numerical computations which support the analytical calculations. Section 7 provides a brief conclusion about the system dynamics.


## 2. Model Formulation

In 2020, Mondal and Samanta [11] analyzed the dynamics of a delayed predator-prey interaction incorporating nonlinear prey refuge function under the influence of fear effect and additional food. Motivated by the work of Mondal and Samanta [11], we have first considered a predator-prey-subsidy model with nonlinear prey refuge function where the prey and subsidy occur in the same habitat and they are both consumed by a single generalist predator according to the following differential equations:(1)dxdt=rx−d1x−a1x2−a21−my/a+yxyb1+w+1−my/a+yx,x0>0,dwdt=A−d2w−a3wyb1+w+1−my/a+yx,w0>0,dydt=c1a21−my/a+yxyb1+w+1−my/a+yx+c2a3wyb1+w+1−my/a+yx−d3y,y0>0,where x is prey population, w denotes the population of the subsidy, and y is the generalist predator which exploits both the prey and subsidy. For example, wolves (predator) consume both deer (prey) and salmon carcasses (subsidy) [37].The term1−my/a+yx represents the quantity of prey available to the predator, i.e., mxy/a+y amounts of prey are free from predation risk where mxy/a+y is designated as nonlinear prey refuge function. Also, we have modeled the dynamics of a generalist predator with Holling type II [38–41] response function in the presence of nonlinear prey refuge function.All parameters are positive (exceptA≥0) and biologically meaningful. Parameters are described in Table 1.Table 1
Description of biologically meaningful parameters.
ParametersDescriptionsrBirth rate of the preyd1Natural death rate of the preyd2The subsidy decay ratea1Mortality rate due to intraspecific competition among the individuals of the prey populationa2Consumption rate of the predatora3Maximum rate at which the predator consumes the subsidyb1Handling time assumed to be uniform over all food sourcesc1Conversion rate of the energy that the predator obtains from the target preyc2Conversion rate of the energy that the predator obtains from the subsidyASubsidy input ratemCoefficient of prey refuge (m∈0,1)d3Mortality rate of the predatoraHalf saturation constant for refuge functionApart from direct consumption, feeling of fear among the individuals of the prey species in presence of predator is very common in predator-prey interaction which changes life-history, behavioural responses, and reproduction capability of prey species. In ecology, effect of fear is a common factor, but there does not exist any considerable attention to introduce the impact of fear in the mathematical modeling. Experimental studies indicate that the feeling of fear among the individuals of the prey species in presence of predator reduces the prey’s birth rate. So, birth rate of prey speciesr is multiplied by a monotone decreasing function fk,y=1/1+ky, where k≥0 is a level of fear [5]. The fear function fk,y satisfies the following conditions:(1)
f0,y=1: when there is no fear effect on the prey species, the birth rate of the prey is not reduced(2)
fk,0=1: when there is no predator, the birth rate of the prey species is not reduced in the presence of fear effect(3)
∂fk,y/∂k<0: when fear effect increases, the birth rate of the prey reduces(4)
∂fk,y/∂y<0: when predator species increases, prey population reducesOur main focus is to analyze the dynamic nature of the predator-prey-subsidy model with the influence of nonlinear prey refuge and fear effect. So, system (1) can be modified in the following aspects:(2)dxdt=rx1+ky−d1x−a1x2−a21−my/a+yxyb1+w+1−my/a+yx,dwdt=A−d2w−a3wyb1+w+1−my/a+yx,dydt=c1a21−my/a+yxyb1+w+1−my/a+yx+c2a3wyb1+w+1−my/a+yx−d3y,with initial conditions:(3)x0>0,w0>0,y0>0.Throughout the analysis of this work, we have takenc1>c2 which is biologically meaningful.


## 3. Positivity and Uniform Boundedness

Theorem 1.
Every solution of system (2) with (3) uniquely exists and is positive for all t≥0.Proof.
Solution (xt,wt,yt) of (2) with (3) exists and is unique on 0,ξ, where 0<ξ≤+∞ [42].
From (2) with (3),(4)xt=x0exp∫0tr1+kyθ−d1−a1xθ−a21−myθ/a+yθyθb1+wθ+1−myθ/a+yθxθdθ>0,yt=y0exp∫0tc1a21−myθ/a+yθxθb1+wθ+1−myθ/a+yθxθ+c2a3wθb1+wθ+1−myθ/a+yθxθ−d3dθ>0.
Now, we claim thatwt>0 for all t∈0,ξ. If it does not hold then there exists t1∈0,ξ such that wt1=0, w˙t≤0, and wt>0 on 0,t1. From the second equation of (2),(5)wt1=w0exp−∫0t1d2+a3yθb1+wθ+1−myθ/a+yθxθdθ+∫0t1Aexp∫t1ud2+a3yθb1+wθ+1−myθ/a+yθxθdθdu>0,a contradiction with wt1=0. So, wt>0,∀t∈0,ξ.
Hence, solutions of (2) stay positive for all t≥0.Theorem 2.
All solutions of system (2) which start in ℝ+3 are asymptotically uniformly bounded.Proof.
Case 1: ifr>d1, from the first equation of (2),(6)dxdt≤rx1+ky−d1x−a1x2≤rx−d1x−a1x2sincerx1+ky≤rx=r−d1x1−xr−d1/a1,∴limsupt⟶+∞xt≤r−d1a1,sincer>d1.Let us takeP=x+w+y/c1.Differentiating both sides with respect tot, we obtain(7)dPdt=dxdt+dwdt+1c1dydt=rx1+ky−d1x−a1x2+A−d2w−a3wya+yb1+wa+y+a+1−myx+c2/c1a3wya+yb1+wa+y+a+1−myx−d3yc1≤rx−d1x−a1x2+A−d2w−d3yc1−a3wya+yb1+wa+y+a+1−myx1−c2c1≤2r−d1x−a1x2+A−r−d1x−d2w−d3yc1sincec1>c2=−a1x−r−d1a12+r−d12a1+A−r−d1x+d2w+d3yc1≤r−d12a1+A−r−d1x+d2w+d3yc1,∴dPdt≤r−d12a1+A−r−d1x+d2w+d3yc1.Let(8)η=minr−d1,d2,d3,whenr>d1.Then,(9)dPdt+ηP≤r−d12a1+A.Using the Gronwall inequality, we obtain(10)0<Pxt,wt,yt≤r−d12/a1+Aη1−e−ηt+e−ηtPx0,w0,y0,∴0<Pxt,wt,yt≤r−d12/a1+Aη,ast⟶∞.Thus, all solutions of system (2) enter into the region:(11)Ω=x,w,y∈ℝ+3:0<xt≤r−d1a1;0<Pxt,wt,yt≤r−d12/a1+Aη.Case 2: ifr<d1, from the first equation of (2) we obtain limt⟶∞xt=0.
Now, from the second and third equations of (2), we have(12)dwdt+1c2dydt=A−d2w+c1/c2a21−my/a+yxyb1+w+1−my/a+yx−d3yc2.
For larget,(13)dwdt+1c2dydt=A−d2w−d3yc2∵limt⟶∞xt=0.
Let(14)η′=mind2,d3.
Then,(15)ddtw+yc2+η′w+yc2≤A.
Using Gronwall inequality, we obtain(16)0<wt+ytc2≤Aη′1−e−η′t+e−η′tw0+y0c2,∴0<wt+ytc2≤Aη′,ast⟶∞.
Hence, the theorem.


## 4. Equilibrium Points and Stability Analysis

### 4.1. Equilibria

#### 4.1.1. Trivial Equilibrium Point

Extinction:E00,0,0.


#### 4.1.2. Axial Equilibrium Points

(i)
Subsidy only:E10,A/d2,0(ii)
Prey only:E2r−d1/a1,0,0 exists if r>d1 and A=0


#### 4.1.3. Planer Equilibrium Points

(i)
Subsidy free:E3x¯,0,y¯ exists if A=0, x¯>0, and y¯>0 where x¯ and y¯ can be obtained by solving the equations:(17)r1+ky−d1−a1x−a21−my/a+yyb1+1−my/a+yx=0,c1a21−my/a+yxb1+1−my/a+yx−d3=0,and we get(18)y¯=ad3b1−x¯c1a2−d31−mc1a2−d3x¯−d3b1,m≠1,wherex¯ is a positive root of the equation:(19)B0x4+B1x3+B2x2+B3x+B4=0.Here,(20)B0=a11−m1−m−akc1a2−d32,B1=a1b1d31−mak−1c1a2−d3−a1b1d31−m−akc1a2−d3+d11−m1−m−akc1a2−d32−r1−m2c1a2−d32,B2=b1d1d31−mak−1c1a2−d3−a1b12d32ak−1−b1d1d31−m−akc1a2−d3−d3c11−m−akc1a2−d32+2b1d3r1−mc1a2−d3,B3=−b12d32r+b1d32c11−m−akc1a2−d3−b1d32c1ak−1c1a2−d3−b12d1d32ak−1,B4=b12d33c1ak−1,andy¯=ad3b1−x¯c1a2−d3/1−mc1a2−d3x¯−d3d1 exists if 1−mc1a2−d3x¯<d3b1<c1a2−d3x¯ and c1a2>d3.(ii)
Prey free:E40,w^=d3b1/c2a3−d3,y^=A−d2w^b1+w^/a3 exists if c2a3>d3 and y^>0.(iii)
Predator free:E5r−d1/a1,A/d2,0 exists if r>d1.


#### 4.1.4. Interior (Coexistence) Equilibrium Point

Solving the following system of equations,(21)r1+ky−d1−a1x−a21−my/a+yyb1+w+1−my/a+yx=0,A−d2w−a3wyb1+w+1−my/a+yx=0,c1a21−my/a+yxb1+w+1−my/a+yx+c2a3wb1+w+1−my/a+yx−d3=0,we can obtain E∗x∗,w∗,y∗ using the software MATHEMATICA with the following existence conditions:(1)
d3<minc1a2,c2a3(2)
0<w∗<minA/d2,d3b1/c2a3−d3(3)
r>d1+a1x∗+b1d3−c2a3−d3w∗A−d2w∗/a3x∗w∗c1a2−d3 (otherwise, predator population goes into extinction)


### 4.2. Stability Analysis

Now, we will study the stability conditions of all equilibria for the proposed system (2).The Jacobian matrixJ00,0,0 at E00,0,0 is given by(22)J00,0,0=r−d1000−d2000−d3.The eigenvalues ofJ00,0,0 are r−d1, −d2<0, and −d3<0. Then. we have stated the following theorem.Theorem 3.
Trivial equilibrium pointE00,0,0 is locally asymptotically stable (LAS) if r<d1 and unstable if r>d1.The Jacobian matrixJ10,A/d2,0 at E10,A/d2,0 is as follows:(23)J10,Ad2,0=r−d1000−d2−a3A/d2b1+A/d200−d3+c2a3A/d2b1+A/d2.We see that(24)r−d1,−d2<0,−d3+c2a3A/d2b1+A/d2,are the eigenvalues of the matrix J10,A/d2,0. Thus, we have the following theorem.Theorem 4.
Axial equilibrium point (subsidy only)E10,A/d2,0 is locally asymptotically stable if r<d1 and c2a3A/d2/b1+A/d2<d3 and unstable if either r>d1,c2a3A/d2/b1+A/d2<d3 or r<d1,c2a3A/d2/b1+A/d2>d3 or r>d1,c2a3A/d2/b1+A/d2>d3.The Jacobian matrixJ2r−d1/a1,0,0 at E2r−d1/a1,0,0 is given by(25)J2r−d1a1,0,0=−r+d10r−d1a1−rk−a2b1+r−d1/a10−d2000−d3+c1a2r−d1/a1b1+r−d1/a1.The eigenvalues ofJ2r−d1/a1,0,0 are −r+d1<0, −d2<0 or −d3+c1a2r−d1/a1/b1+r−d1/a1. Thus. we have stated the following theorem.Theorem 5.
Axial equilibrium point (prey only)E2r−d1/a1,0,0 is locally asymptotically stable if(26)c1a2r−d1/a1b1+r−d1/a1<d3⟹0<r<d1+a1b1d3c1a2−d3withc1a2>d3.
Ifr>d1+a1b1d3/c1a2−d3 with c1a2>d3, then E2r−d1/a1,0,0 is unstable.The Jacobian matrixJ3x¯,0,y¯ at E3x¯,0,y¯ is as follows:(27)J3x¯,0,y¯=b11b12b13b21b22b23b31b32b33,where(28)b11=x¯−a1+a2y¯a+1−my¯2b1a+y¯+a+1−my¯x¯2,b12=a2x¯y¯a+1−my¯a+y¯b1a+y¯+a+1−my¯x¯2,b13=x¯−rk1+ky¯2−a2a+21−my¯b1a+y¯+a+1−my¯x¯+a2y¯a+1−my¯b1+1−mx¯b1a+y¯+a+1−my¯x¯2,b21=0,b22=−d2−a3y¯b1+1−my¯/a+y¯x¯,b23=0,b31=y¯c1a2a+1−my¯b1a+y¯+a+1−my¯x¯−c1a2a+1−my¯2x¯b1a+y¯+a+1−my¯x¯2,b32=c2a3y¯a+y¯b1a+y∗+a+1−my¯x¯,b33=y¯c1a21−mx¯b1a+y¯+a+1−my¯x¯−c1a2a+1−my¯x¯b1+1−mx¯b1a+y¯+a+1−my¯x¯2.The characteristic equation corresponding toJ3x¯,0,y¯ is expressed as(29)λ3+B1λ2+B2λ+B3=0,where B1=−b11+b22+b33, B2=b22b33+b11b33−b13b31+b11b22, and B3=−b22b11b33−b13b31.Theorem 6.
Subsidy-free equilibrium pointE3x¯,0,y¯ is locally asymptotically stable if b11<0, b33<0, and b13<0.The Jacobian matrixJ40,w^,y^ corresponding to E40,w^,y^ is given by(30)J40,w^,y^=c11c12c13c21c22c23c31c32c33,where(31)c11=r1+ky^−d1−a2a+1−my^y^b1+w^a+y^,c12=0,c13=0,c21=a31−ma+y^w^y^2b1+w^a+y^2,c22=−d2−a3y^a+y^b1+w^a+y^,c23=−a3w^a+2y^b1+w^a+y^+a3y^w^a+y^b1+w^b1+w^a+y^2,c31=c1a2y^a+1−my^b1+w^a+y^,c32=c2a3y^a+y^b1+w^a+y^−c2a3w^y^a+y^2b1+w^a+y^2=c2a3y^b1b1+w^a+y^2,c33=c2a3w^y^b1+w^a+y^−c2a3w^y^a+y^b1+w^b1+w^a+y^2=0.The characteristic equation corresponding toJ40,w^,y^ is expressed as(32)λ3+C1λ2+C2λ+C3=0,where(33)C1=−c11+c22,C2=−c23c32+c11c22,C3=+c11c23c32.Theorem 7.
Prey-free equilibrium pointE40,w^,y^ is locally asymptotically stable if c11<0 and c23<0.The Jacobian matrixJ5r−d1/a1,A/d2,0 corresponding to E5r−d1/a1,A/d2,0 is given by(34)J5r−d1a1,Ad2,0=d11d12d13d21d22d23d31d32d33,where(35)d11=−r+d1,d12=0,d13=r−d1a1−rk−a2b1+A/d2+r−d1/a1,d21=0,d22=−d2,d23=−a3A/d2b1+A/d2+r−d1/a1,d31=d32=0,d33=c1a2r−d1/a1b1+A/d2+r−d1/a1+c2a3A/d2b1+A/d2+r−d1/a1−d3.The characteristic equation corresponding toJ5r−d1/a1,A/d2,0 is expressed as(36)λ3+D1λ2+D2λ+D3=0,where(37)D1=−d11+d22+d33,D2=d22d33+d11d33+d11d22,D3=−d11d22d33.Theorem 8.
Predator free equilibrium pointE5r−d1/a1,A/d2,0 is locally asymptotically stable if(38)c1a2r−d1/a1b1+A/d2+r−d1/a1+c2a3A/d2b1+A/d2+r−d1/a1<d3⟹0<r<d1+a1b1d3−c2a3−d3A/d2c1a2−d3,provided d3<minc1a2,c2a3 and 0<A/d2<b1d3/c2a3−d3.The Jacobian matrixJ∗x∗,w∗,y∗ corresponding to E∗x∗,w∗,y∗ is as follows:(39)J∗x∗,w∗,y∗=a11a12a13a21a22a23a31a32a33,where(40)a11=x∗−a1+a2a+1−my∗2y∗b1+w∗a+y∗+a+1−my∗x∗2,a12=a2a+1−my∗a+y∗x∗y∗b1+w∗a+y∗+a+1−my∗x∗2,a13=−rkx∗1+ky∗2−a2a+21−my∗x∗b1+w∗a+y∗+a+1−my∗x∗+a2a+1−my∗x∗y∗b1+w∗+1−mx∗b1+w∗a+y∗+a+1−my∗x∗2,a21=a3a+y∗a+1−my∗w∗y∗b1+w∗a+y∗+a+1−my∗x∗2,a22=−d2−a3a+y∗y∗b1+w∗a+y∗+a+1−my∗x∗+a3a+y∗2w∗y∗b1+w∗a+y∗+a+1−my∗x∗2,a23=−a3a+2y∗w∗b1+w∗a+y∗+a+1−my∗x∗+a3a+y∗b1+w∗+1−mx∗b1+w∗a+y∗+a+1−my∗x∗2,a31=c1a2a+1−my∗y∗b1+w∗a+y∗+a+1−my∗x∗−c1a2a+1−my∗2x∗y∗b1+w∗a+y∗+a+1−my∗x∗2−c2a3w∗y∗a+y∗a+1−my∗b1+w∗a+y∗+a+1−my∗x∗2,a32=−c1a2a+1−my∗x∗y∗a+y∗b1+w∗a+y∗+a+1−my∗x∗2+c2a3a+y∗y∗b1+w∗a+y∗+a+1−my∗x∗−c2a3w∗y∗a+y∗2b1+w∗a+y∗+a+1−my∗x∗2,a33=c1a21−mx∗y∗+c2a3w∗y∗b1+w∗a+y∗+a+1−my∗x∗−c1a2a+1−my∗x∗y∗+c2a3w∗y∗a+y∗b1+w∗+1−mx∗b1+w∗a+y∗+a+1−my∗x∗2.The characteristic equation corresponding toJ∗x∗,w∗,y∗ is expressed as(41)λ3+A1λ2+A2λ+A3=0,where(42)A1=−a11+a22+a33,A2=a22a33−a23a32+a11a33−a13a31+a11a22−a12a21,A3=−a11a22a33−a23a32+a12a23a31−a21a33+a13a21a32−a22a31.Theorem 9.
The coexistence equilibriumE∗x∗,w∗,y∗ is locally asymptotically stable if A1>0, A3>0, and A1A2>A3, where A1,A2,andA3 are stated in (41).


### 4.3. Local Bifurcations of Codimension 1

#### 4.3.1. Transcritical Bifurcation

Theorem 10.
System (2) undergoes a transcritical bifurcation around E10,A/d2,0 if d1TC=r and c2a3A/d2/b1+A/d2<d3 (TC stands for transcritical bifurcation).Proof.
We apply Sotomayor’s theorem [43] to prove the occurrence of a transcritical bifurcation around E1 with d1 as bifurcation parameter. For applicability of Sotomayor’s theorem, exactly one of the eigenvalues of the Jacobian matrix at E1 must be zero and other eigenvalues must have negative real parts. So, we need to fulfill the condition c2a3A/d2/b1+A/d2<d3.
The eigenvectors ofJE1=vij and JE1T corresponding to the zero eigenvalue of E10,A/d2,0 are obtained as V=v1,v2,v3T and W=1,0,0T, respectively, where v1=1, v2=v3=0, v11=0, v12=0, v13=0, v21=0, v22=−d2, v23=0, v31=0, v32=0, v33=−d3+c2a3A/d2/b1+A/d2.
ComputeΔ1, Δ2, and Δ3 as follows:(43)Δ1=WT⋅Fd10,Ad2,0;d1TC=1,0,0⋅∂F1∂d1∂F2∂d1∂F3∂d1E1=1,0,0⋅−x00E1=0,where F=F1,F2,F3T and F1, F2, and F3 are given by(44)F1=rx1+ky−d1x−a1x2−a21−my/a+yxyb1+w+1−my/a+yx,F2=A−d2w−a3wyb1+w+1−my/a+yx,F3=c1a21−my/a+yxyb1+w+1−my/a+yx+c2a3wyb1+w+1−my/a+yx−d3y,Δ2=WT⋅DFd10,Ad2,0;d1TCV=1,0,0⋅∂2F1∂x∂d1∂2F1∂w∂d1∂2F1∂y∂d1∂2F2∂x∂d1∂2F2∂w∂d1∂2F2∂y∂d1∂2F3∂x∂d1∂2F3∂w∂d1∂2F3∂y∂d1E1⋅100=1,0,0⋅−100000000E1⋅100=−1≠0,Δ3=WT⋅D2F0,Ad2,0;d1TCV,V=1,0,0⋅D∂F1∂xv1+∂F1∂wv2+∂F1∂yv3∂F2∂xv1+∂F2∂wv2+∂F2∂yv3∂F3∂xv1+∂F3∂wv2+∂F3∂yv3E1⋅v1v2v3=1,0,0⋅∂2F1∂x2v12+∂2F1∂w2v22+∂2F1∂y2v32+2∂2F1∂x∂wv1v2+2∂2F1∂x∂yv1v3+2∂2F1∂w∂yv2v3∂2F2∂x2v12+∂2F2∂w2v22+∂2F2∂y2v32+2∂2F2∂x∂wv1v2+2∂2F2∂x∂yv1v3+2∂2F2∂w∂yv2v3∂2F3∂x2v12+∂2F3∂w2v22+∂2F3∂y2v32+2∂2F3∂x∂wv1v2+2∂2F3∂x∂yv1v3+2∂2F3∂w∂yv2v3E1=−2a1≠0.
Therefore, by Sotomayor’s theorem [43], system (2) undergoes a transcritical bifurcation at d1=d1TC around the axial equilibrium point E1.Theorem 11.
System (2) exhibits a transcritical bifurcation around E2r−d1/a1,0,0 if(45)d3TC=c1a2r−d1/a1b1+r−d1/a1.Proof.
Let us apply Sotomayor’s theorem [43] to prove the occurrence of a transcritical bifurcation around E2 with d3 as bifurcation parameter. For applicability of Sotomayor’s theorem, exactly one of the eigenvalues of the Jacobian matrix at E2 must be zero and other eigenvalues must have negative real parts.
The eigenvectors ofJE2=tij and JE2T corresponding to the zero eigenvalue of E2r−d1/a1,0,0 are obtained as V=v1,v2,v3T and W=0,0,1T, respectively, where v1=1/a1−rk−a2/b1+r−d1/a1<0, v2=0, v3=1, t11=−r+d1, t12=0, t13=r−d1/a1−rk−a2/b1+r−d1/a1, t21=0, t22=−d2, t23=0, t31=0, t32=0, t33=0.
ComputeΔ1, Δ2, and Δ3 as follows:(46)Δ1=WT⋅Fd1r−d1a1,0,0;d3TC=0,0,1⋅∂F1∂d3∂F2∂d3∂F3∂d3E2=0,0,1⋅00−yE2=0,where F=F1,F2,F3T and F1, F2, and F3 are given by(47)F1=rx1+ky−d1x−a1x2−a21−my/a+yxyb1+w+1−my/a+yx,F2=A−d2w−a3wyb1+w+1−my/a+yx,F3=c1a21−my/a+yxyb1+w+1−my/a+yx+c2a3wyb1+w+1−my/a+yx−d3y,Δ2=WT⋅DFd3r−d1a1,0,0;d3TCV=0,0,1⋅∂2F1∂x∂d3∂2F1∂w∂d3∂2F1∂y∂d3∂2F2∂x∂d3∂2F2∂w∂d3∂2F2∂y∂d3∂2F3∂x∂d3∂2F3∂w∂d3∂2F3∂w∂d3E2⋅v1v2v3=0,0,1⋅00000000−1E2⋅v1v2v3=−v3=−1≠0,Δ3=WT⋅D2Fr−d1a1,0,0;d3TCV,V=0,0,1⋅D∂F1∂xv1+∂F1∂wv2+∂F1∂yv3∂F2∂xv1+∂F2∂wv2+∂F2∂yv3∂F3∂xv1+∂F3∂wv2+∂F3∂yv3E2⋅v1v2v3=0,0,1⋅∂2F1∂x2v12+∂2F1∂w2v22+∂2F1∂y2v32+2∂2F1∂x∂wv1v2+2∂2F1∂x∂yv1v3+2∂2F1∂w∂yv2v3∂2F2∂x2v12+∂2F2∂w2v22+∂2F2∂y2v32+2∂2F2∂x∂wv1v2+2∂2F2∂x∂yv1v3+2∂2F2∂w∂yv2v3∂2F3∂x2v12+∂2F3∂w2v22+∂2F3∂y2v32+2∂2F3∂x∂wv1v2+2∂2F3∂x∂yv1v3+2∂2F3∂w∂yv2v3E2=2a2b1c1v1−m/ab1+r−d1/a12≠0.
Therefore, by Sotomayor’s theorem [43], system (2) exhibits a transcritical bifurcation at d3=d3TC around the axial equilibrium point E2.Theorem 12.
System (2) undergoes a transcritical bifurcation around E40,w^,y^ if(48)d1TC=r1+ky^−a2a+1−my^y^b1+y^a+y^,c2a3w^y^b1+w^a+y^<c2a3w^y^a+y^b1+w^b1+w^a+y^2.Proof.
Proof is the same as in Theorem10.Theorem 13.
System (2) undergoes a transcritical bifurcation around E5r−d1/a1,A/d2,0 if(49)d3TC=c1a2r−d1/a1b1+A/d2+r−d1/a1+c2a3A/d2b1+A/d2+r−d1/a1.Proof.
Proof is the same as in Theorem11.


#### 4.3.2. Hopf Bifurcation aroundE∗x∗,w∗,y∗

Let us considerk as a bifurcation parameter of system (2) where the characteristic equation at E∗ is(50)λ3+A1kλ2+A2kλ+A3k=0.Then, Hopf bifurcation theorem is stated as follows.Theorem 14.
(Hopf bifurcation theorem [44]). If A1k, A2k, and A3k are the smooth functions of k in NεkH, ε>0, kH∈ℝ for which the characteristic equation (50) has the following:(i)
A pair of imaginary eigenvaluesλ=q1k±iq2k with q1k and q2k∈ℝ so that they become purely complex at k=kH and dq1/dk|k=kH≠0(ii)
The other eigenvalue is negative atk=kH; then, a Hopf bifurcation appears around E∗ at k=kHTheorem 15.
IfE∗x∗,w∗,y∗ is locally asymptotically stable, then a Hopf bifurcation is exhibited around E∗x∗,w∗,y∗ when k passes through its critical value kH provided A1kH>0, A3kH>0, and A1kHA2kH=A3kH (kH is a positive root of equation A1kA2k−A3k=0).Proof.
Atk=kH, we can write equation (50) as(51)λ2+A2λ+A1=0.
The roots of equation (51) are λ1=iA2, λ2=−iA2, and λ3=−A1. Also A1, A2, and A3 are the smooth functions of k. So, the roots of equation (59) have the form λ1=p1k+ip2k, λ2=p1k−ip2k, and λ3=p3k where pik are real functions of k in an open neighborhood of kH for i=1,2,3. Next, we verify the transversality condition:(52)ddkReλik|k=kH≠0,i=1,2.
Puttingλk=p1k+ip2k in (59), we get(53)p1k+ip2k3+A1kp1k+ip2k2+A2kp1k+ip2k+A3k=0.
Differentiating both sides with respect tok, we have(54)3p1k+ip2k2p˙1k+ip˙2k+2A1kp1k+ip2kp˙1k+ip˙2k+A˙1kp1k+ip2k2+A2p˙1k+ip˙2k+A˙2kp1k+ip2k+A˙3k=0.
Comparing real and imaginary parts from both sides, we obtain(55)X1p˙1−X2p˙2+X3=0,(56)X2p˙1+X1p˙2+X4=0,where(57)X1=3p12−p22+2A1p1+A2,X2=6p1p2+2A1p2,X3=A˙1p12−p22+A˙2p1+A˙3≠A˙1p12−p22+A˙2p1+A1A˙2+A2A˙1,sinceA3≠A1A2in a deleted neighbour hood ofkH,and X4=2A˙1p1p2+A˙2p2.
Multiplying (55) by X1 and (56) by X2 and then adding, we get(58)X12+X22p˙1+X1X3+X2X4=0⟹p˙1=−X1X3+X2X4X12+X22.
Atk=kH,Case 1:p1=0, p2=A2. Then, X1=−2A2, X2=2A1A2, X3≠A1A˙2, and X4=A˙2A2. ∴X1X3+X2X4≠0.Case 2:p1=0, p2=−A2. Then, X1=−2A2, X2=−2A1A2, X3≠A1A˙2, and X4=−A˙2A2. ∴X1X3+X2X4≠0.
Also,λ3=−A1kH<0.
Hence, this theorem is proved by virtue of Theorem14.


#### 4.3.3. Hopf Bifurcation aroundE3x¯,0,y¯

Let us considerk as a bifurcation parameter of system (2) where the characteristic equation of E3 is(59)λ3+B1kλ2+B2kλ+B3k=0,and then Hopf bifurcation theorem is stated as follows.Theorem 16.
(Hopf bifurcation theorem [44]). If B1k, B2k, and B3k are the smooth functions of k in Nεk∗, ε>0, k∗∈ℝ for which the characteristic equation (59) has the following:(i)
A pair of imaginary eigenvaluesλ=p1′k±ip2′k with p1′k and p2′k∈ℝ so that they become purely imaginary at k=k∗ and dp1′/dk|k=k∗≠0(ii)
The other eigenvalue is negative atk=k∗; then a Hopf bifurcation occurs around E3x¯,0,y¯ at k=k∗Theorem 17.
IfE3x¯,0,y¯ is locally asymptotically stable, then a Hopf bifurcation appears around subsidy-free equilibrium E3x¯,0,y¯ when k passes through its critical value k∗ provided B1k∗>0, B3k∗>0, and B1k∗B2k∗=B3k∗ (k∗ is a positive root of equation B1kB2k−B3k=0).Proof.
Proof is the same as in Theorem15.


## 4.1. Equilibria

### 4.1.1. Trivial Equilibrium Point

Extinction:E00,0,0.


### 4.1.2. Axial Equilibrium Points

(i)
Subsidy only:E10,A/d2,0(ii)
Prey only:E2r−d1/a1,0,0 exists if r>d1 and A=0


### 4.1.3. Planer Equilibrium Points

(i)
Subsidy free:E3x¯,0,y¯ exists if A=0, x¯>0, and y¯>0 where x¯ and y¯ can be obtained by solving the equations:(17)r1+ky−d1−a1x−a21−my/a+yyb1+1−my/a+yx=0,c1a21−my/a+yxb1+1−my/a+yx−d3=0,and we get(18)y¯=ad3b1−x¯c1a2−d31−mc1a2−d3x¯−d3b1,m≠1,wherex¯ is a positive root of the equation:(19)B0x4+B1x3+B2x2+B3x+B4=0.Here,(20)B0=a11−m1−m−akc1a2−d32,B1=a1b1d31−mak−1c1a2−d3−a1b1d31−m−akc1a2−d3+d11−m1−m−akc1a2−d32−r1−m2c1a2−d32,B2=b1d1d31−mak−1c1a2−d3−a1b12d32ak−1−b1d1d31−m−akc1a2−d3−d3c11−m−akc1a2−d32+2b1d3r1−mc1a2−d3,B3=−b12d32r+b1d32c11−m−akc1a2−d3−b1d32c1ak−1c1a2−d3−b12d1d32ak−1,B4=b12d33c1ak−1,andy¯=ad3b1−x¯c1a2−d3/1−mc1a2−d3x¯−d3d1 exists if 1−mc1a2−d3x¯<d3b1<c1a2−d3x¯ and c1a2>d3.(ii)
Prey free:E40,w^=d3b1/c2a3−d3,y^=A−d2w^b1+w^/a3 exists if c2a3>d3 and y^>0.(iii)
Predator free:E5r−d1/a1,A/d2,0 exists if r>d1.


### 4.1.4. Interior (Coexistence) Equilibrium Point

Solving the following system of equations,(21)r1+ky−d1−a1x−a21−my/a+yyb1+w+1−my/a+yx=0,A−d2w−a3wyb1+w+1−my/a+yx=0,c1a21−my/a+yxb1+w+1−my/a+yx+c2a3wb1+w+1−my/a+yx−d3=0,we can obtain E∗x∗,w∗,y∗ using the software MATHEMATICA with the following existence conditions:(1)
d3<minc1a2,c2a3(2)
0<w∗<minA/d2,d3b1/c2a3−d3(3)
r>d1+a1x∗+b1d3−c2a3−d3w∗A−d2w∗/a3x∗w∗c1a2−d3 (otherwise, predator population goes into extinction)


## 4.1.1. Trivial Equilibrium Point

Extinction:E00,0,0.


## 4.1.2. Axial Equilibrium Points

(i)
Subsidy only:E10,A/d2,0(ii)
Prey only:E2r−d1/a1,0,0 exists if r>d1 and A=0


## 4.1.3. Planer Equilibrium Points

(i)
Subsidy free:E3x¯,0,y¯ exists if A=0, x¯>0, and y¯>0 where x¯ and y¯ can be obtained by solving the equations:(17)r1+ky−d1−a1x−a21−my/a+yyb1+1−my/a+yx=0,c1a21−my/a+yxb1+1−my/a+yx−d3=0,and we get(18)y¯=ad3b1−x¯c1a2−d31−mc1a2−d3x¯−d3b1,m≠1,wherex¯ is a positive root of the equation:(19)B0x4+B1x3+B2x2+B3x+B4=0.Here,(20)B0=a11−m1−m−akc1a2−d32,B1=a1b1d31−mak−1c1a2−d3−a1b1d31−m−akc1a2−d3+d11−m1−m−akc1a2−d32−r1−m2c1a2−d32,B2=b1d1d31−mak−1c1a2−d3−a1b12d32ak−1−b1d1d31−m−akc1a2−d3−d3c11−m−akc1a2−d32+2b1d3r1−mc1a2−d3,B3=−b12d32r+b1d32c11−m−akc1a2−d3−b1d32c1ak−1c1a2−d3−b12d1d32ak−1,B4=b12d33c1ak−1,andy¯=ad3b1−x¯c1a2−d3/1−mc1a2−d3x¯−d3d1 exists if 1−mc1a2−d3x¯<d3b1<c1a2−d3x¯ and c1a2>d3.(ii)
Prey free:E40,w^=d3b1/c2a3−d3,y^=A−d2w^b1+w^/a3 exists if c2a3>d3 and y^>0.(iii)
Predator free:E5r−d1/a1,A/d2,0 exists if r>d1.


## 4.1.4. Interior (Coexistence) Equilibrium Point

Solving the following system of equations,(21)r1+ky−d1−a1x−a21−my/a+yyb1+w+1−my/a+yx=0,A−d2w−a3wyb1+w+1−my/a+yx=0,c1a21−my/a+yxb1+w+1−my/a+yx+c2a3wb1+w+1−my/a+yx−d3=0,we can obtain E∗x∗,w∗,y∗ using the software MATHEMATICA with the following existence conditions:(1)
d3<minc1a2,c2a3(2)
0<w∗<minA/d2,d3b1/c2a3−d3(3)
r>d1+a1x∗+b1d3−c2a3−d3w∗A−d2w∗/a3x∗w∗c1a2−d3 (otherwise, predator population goes into extinction)


## 4.2. Stability Analysis

Now, we will study the stability conditions of all equilibria for the proposed system (2).The Jacobian matrixJ00,0,0 at E00,0,0 is given by(22)J00,0,0=r−d1000−d2000−d3.The eigenvalues ofJ00,0,0 are r−d1, −d2<0, and −d3<0. Then. we have stated the following theorem.Theorem 3.
Trivial equilibrium pointE00,0,0 is locally asymptotically stable (LAS) if r<d1 and unstable if r>d1.The Jacobian matrixJ10,A/d2,0 at E10,A/d2,0 is as follows:(23)J10,Ad2,0=r−d1000−d2−a3A/d2b1+A/d200−d3+c2a3A/d2b1+A/d2.We see that(24)r−d1,−d2<0,−d3+c2a3A/d2b1+A/d2,are the eigenvalues of the matrix J10,A/d2,0. Thus, we have the following theorem.Theorem 4.
Axial equilibrium point (subsidy only)E10,A/d2,0 is locally asymptotically stable if r<d1 and c2a3A/d2/b1+A/d2<d3 and unstable if either r>d1,c2a3A/d2/b1+A/d2<d3 or r<d1,c2a3A/d2/b1+A/d2>d3 or r>d1,c2a3A/d2/b1+A/d2>d3.The Jacobian matrixJ2r−d1/a1,0,0 at E2r−d1/a1,0,0 is given by(25)J2r−d1a1,0,0=−r+d10r−d1a1−rk−a2b1+r−d1/a10−d2000−d3+c1a2r−d1/a1b1+r−d1/a1.The eigenvalues ofJ2r−d1/a1,0,0 are −r+d1<0, −d2<0 or −d3+c1a2r−d1/a1/b1+r−d1/a1. Thus. we have stated the following theorem.Theorem 5.
Axial equilibrium point (prey only)E2r−d1/a1,0,0 is locally asymptotically stable if(26)c1a2r−d1/a1b1+r−d1/a1<d3⟹0<r<d1+a1b1d3c1a2−d3withc1a2>d3.
Ifr>d1+a1b1d3/c1a2−d3 with c1a2>d3, then E2r−d1/a1,0,0 is unstable.The Jacobian matrixJ3x¯,0,y¯ at E3x¯,0,y¯ is as follows:(27)J3x¯,0,y¯=b11b12b13b21b22b23b31b32b33,where(28)b11=x¯−a1+a2y¯a+1−my¯2b1a+y¯+a+1−my¯x¯2,b12=a2x¯y¯a+1−my¯a+y¯b1a+y¯+a+1−my¯x¯2,b13=x¯−rk1+ky¯2−a2a+21−my¯b1a+y¯+a+1−my¯x¯+a2y¯a+1−my¯b1+1−mx¯b1a+y¯+a+1−my¯x¯2,b21=0,b22=−d2−a3y¯b1+1−my¯/a+y¯x¯,b23=0,b31=y¯c1a2a+1−my¯b1a+y¯+a+1−my¯x¯−c1a2a+1−my¯2x¯b1a+y¯+a+1−my¯x¯2,b32=c2a3y¯a+y¯b1a+y∗+a+1−my¯x¯,b33=y¯c1a21−mx¯b1a+y¯+a+1−my¯x¯−c1a2a+1−my¯x¯b1+1−mx¯b1a+y¯+a+1−my¯x¯2.The characteristic equation corresponding toJ3x¯,0,y¯ is expressed as(29)λ3+B1λ2+B2λ+B3=0,where B1=−b11+b22+b33, B2=b22b33+b11b33−b13b31+b11b22, and B3=−b22b11b33−b13b31.Theorem 6.
Subsidy-free equilibrium pointE3x¯,0,y¯ is locally asymptotically stable if b11<0, b33<0, and b13<0.The Jacobian matrixJ40,w^,y^ corresponding to E40,w^,y^ is given by(30)J40,w^,y^=c11c12c13c21c22c23c31c32c33,where(31)c11=r1+ky^−d1−a2a+1−my^y^b1+w^a+y^,c12=0,c13=0,c21=a31−ma+y^w^y^2b1+w^a+y^2,c22=−d2−a3y^a+y^b1+w^a+y^,c23=−a3w^a+2y^b1+w^a+y^+a3y^w^a+y^b1+w^b1+w^a+y^2,c31=c1a2y^a+1−my^b1+w^a+y^,c32=c2a3y^a+y^b1+w^a+y^−c2a3w^y^a+y^2b1+w^a+y^2=c2a3y^b1b1+w^a+y^2,c33=c2a3w^y^b1+w^a+y^−c2a3w^y^a+y^b1+w^b1+w^a+y^2=0.The characteristic equation corresponding toJ40,w^,y^ is expressed as(32)λ3+C1λ2+C2λ+C3=0,where(33)C1=−c11+c22,C2=−c23c32+c11c22,C3=+c11c23c32.Theorem 7.
Prey-free equilibrium pointE40,w^,y^ is locally asymptotically stable if c11<0 and c23<0.The Jacobian matrixJ5r−d1/a1,A/d2,0 corresponding to E5r−d1/a1,A/d2,0 is given by(34)J5r−d1a1,Ad2,0=d11d12d13d21d22d23d31d32d33,where(35)d11=−r+d1,d12=0,d13=r−d1a1−rk−a2b1+A/d2+r−d1/a1,d21=0,d22=−d2,d23=−a3A/d2b1+A/d2+r−d1/a1,d31=d32=0,d33=c1a2r−d1/a1b1+A/d2+r−d1/a1+c2a3A/d2b1+A/d2+r−d1/a1−d3.The characteristic equation corresponding toJ5r−d1/a1,A/d2,0 is expressed as(36)λ3+D1λ2+D2λ+D3=0,where(37)D1=−d11+d22+d33,D2=d22d33+d11d33+d11d22,D3=−d11d22d33.Theorem 8.
Predator free equilibrium pointE5r−d1/a1,A/d2,0 is locally asymptotically stable if(38)c1a2r−d1/a1b1+A/d2+r−d1/a1+c2a3A/d2b1+A/d2+r−d1/a1<d3⟹0<r<d1+a1b1d3−c2a3−d3A/d2c1a2−d3,provided d3<minc1a2,c2a3 and 0<A/d2<b1d3/c2a3−d3.The Jacobian matrixJ∗x∗,w∗,y∗ corresponding to E∗x∗,w∗,y∗ is as follows:(39)J∗x∗,w∗,y∗=a11a12a13a21a22a23a31a32a33,where(40)a11=x∗−a1+a2a+1−my∗2y∗b1+w∗a+y∗+a+1−my∗x∗2,a12=a2a+1−my∗a+y∗x∗y∗b1+w∗a+y∗+a+1−my∗x∗2,a13=−rkx∗1+ky∗2−a2a+21−my∗x∗b1+w∗a+y∗+a+1−my∗x∗+a2a+1−my∗x∗y∗b1+w∗+1−mx∗b1+w∗a+y∗+a+1−my∗x∗2,a21=a3a+y∗a+1−my∗w∗y∗b1+w∗a+y∗+a+1−my∗x∗2,a22=−d2−a3a+y∗y∗b1+w∗a+y∗+a+1−my∗x∗+a3a+y∗2w∗y∗b1+w∗a+y∗+a+1−my∗x∗2,a23=−a3a+2y∗w∗b1+w∗a+y∗+a+1−my∗x∗+a3a+y∗b1+w∗+1−mx∗b1+w∗a+y∗+a+1−my∗x∗2,a31=c1a2a+1−my∗y∗b1+w∗a+y∗+a+1−my∗x∗−c1a2a+1−my∗2x∗y∗b1+w∗a+y∗+a+1−my∗x∗2−c2a3w∗y∗a+y∗a+1−my∗b1+w∗a+y∗+a+1−my∗x∗2,a32=−c1a2a+1−my∗x∗y∗a+y∗b1+w∗a+y∗+a+1−my∗x∗2+c2a3a+y∗y∗b1+w∗a+y∗+a+1−my∗x∗−c2a3w∗y∗a+y∗2b1+w∗a+y∗+a+1−my∗x∗2,a33=c1a21−mx∗y∗+c2a3w∗y∗b1+w∗a+y∗+a+1−my∗x∗−c1a2a+1−my∗x∗y∗+c2a3w∗y∗a+y∗b1+w∗+1−mx∗b1+w∗a+y∗+a+1−my∗x∗2.The characteristic equation corresponding toJ∗x∗,w∗,y∗ is expressed as(41)λ3+A1λ2+A2λ+A3=0,where(42)A1=−a11+a22+a33,A2=a22a33−a23a32+a11a33−a13a31+a11a22−a12a21,A3=−a11a22a33−a23a32+a12a23a31−a21a33+a13a21a32−a22a31.Theorem 9.
The coexistence equilibriumE∗x∗,w∗,y∗ is locally asymptotically stable if A1>0, A3>0, and A1A2>A3, where A1,A2,andA3 are stated in (41).


## 4.3. Local Bifurcations of Codimension 1

### 4.3.1. Transcritical Bifurcation

Theorem 10.
System (2) undergoes a transcritical bifurcation around E10,A/d2,0 if d1TC=r and c2a3A/d2/b1+A/d2<d3 (TC stands for transcritical bifurcation).Proof.
We apply Sotomayor’s theorem [43] to prove the occurrence of a transcritical bifurcation around E1 with d1 as bifurcation parameter. For applicability of Sotomayor’s theorem, exactly one of the eigenvalues of the Jacobian matrix at E1 must be zero and other eigenvalues must have negative real parts. So, we need to fulfill the condition c2a3A/d2/b1+A/d2<d3.
The eigenvectors ofJE1=vij and JE1T corresponding to the zero eigenvalue of E10,A/d2,0 are obtained as V=v1,v2,v3T and W=1,0,0T, respectively, where v1=1, v2=v3=0, v11=0, v12=0, v13=0, v21=0, v22=−d2, v23=0, v31=0, v32=0, v33=−d3+c2a3A/d2/b1+A/d2.
ComputeΔ1, Δ2, and Δ3 as follows:(43)Δ1=WT⋅Fd10,Ad2,0;d1TC=1,0,0⋅∂F1∂d1∂F2∂d1∂F3∂d1E1=1,0,0⋅−x00E1=0,where F=F1,F2,F3T and F1, F2, and F3 are given by(44)F1=rx1+ky−d1x−a1x2−a21−my/a+yxyb1+w+1−my/a+yx,F2=A−d2w−a3wyb1+w+1−my/a+yx,F3=c1a21−my/a+yxyb1+w+1−my/a+yx+c2a3wyb1+w+1−my/a+yx−d3y,Δ2=WT⋅DFd10,Ad2,0;d1TCV=1,0,0⋅∂2F1∂x∂d1∂2F1∂w∂d1∂2F1∂y∂d1∂2F2∂x∂d1∂2F2∂w∂d1∂2F2∂y∂d1∂2F3∂x∂d1∂2F3∂w∂d1∂2F3∂y∂d1E1⋅100=1,0,0⋅−100000000E1⋅100=−1≠0,Δ3=WT⋅D2F0,Ad2,0;d1TCV,V=1,0,0⋅D∂F1∂xv1+∂F1∂wv2+∂F1∂yv3∂F2∂xv1+∂F2∂wv2+∂F2∂yv3∂F3∂xv1+∂F3∂wv2+∂F3∂yv3E1⋅v1v2v3=1,0,0⋅∂2F1∂x2v12+∂2F1∂w2v22+∂2F1∂y2v32+2∂2F1∂x∂wv1v2+2∂2F1∂x∂yv1v3+2∂2F1∂w∂yv2v3∂2F2∂x2v12+∂2F2∂w2v22+∂2F2∂y2v32+2∂2F2∂x∂wv1v2+2∂2F2∂x∂yv1v3+2∂2F2∂w∂yv2v3∂2F3∂x2v12+∂2F3∂w2v22+∂2F3∂y2v32+2∂2F3∂x∂wv1v2+2∂2F3∂x∂yv1v3+2∂2F3∂w∂yv2v3E1=−2a1≠0.
Therefore, by Sotomayor’s theorem [43], system (2) undergoes a transcritical bifurcation at d1=d1TC around the axial equilibrium point E1.Theorem 11.
System (2) exhibits a transcritical bifurcation around E2r−d1/a1,0,0 if(45)d3TC=c1a2r−d1/a1b1+r−d1/a1.Proof.
Let us apply Sotomayor’s theorem [43] to prove the occurrence of a transcritical bifurcation around E2 with d3 as bifurcation parameter. For applicability of Sotomayor’s theorem, exactly one of the eigenvalues of the Jacobian matrix at E2 must be zero and other eigenvalues must have negative real parts.
The eigenvectors ofJE2=tij and JE2T corresponding to the zero eigenvalue of E2r−d1/a1,0,0 are obtained as V=v1,v2,v3T and W=0,0,1T, respectively, where v1=1/a1−rk−a2/b1+r−d1/a1<0, v2=0, v3=1, t11=−r+d1, t12=0, t13=r−d1/a1−rk−a2/b1+r−d1/a1, t21=0, t22=−d2, t23=0, t31=0, t32=0, t33=0.
ComputeΔ1, Δ2, and Δ3 as follows:(46)Δ1=WT⋅Fd1r−d1a1,0,0;d3TC=0,0,1⋅∂F1∂d3∂F2∂d3∂F3∂d3E2=0,0,1⋅00−yE2=0,where F=F1,F2,F3T and F1, F2, and F3 are given by(47)F1=rx1+ky−d1x−a1x2−a21−my/a+yxyb1+w+1−my/a+yx,F2=A−d2w−a3wyb1+w+1−my/a+yx,F3=c1a21−my/a+yxyb1+w+1−my/a+yx+c2a3wyb1+w+1−my/a+yx−d3y,Δ2=WT⋅DFd3r−d1a1,0,0;d3TCV=0,0,1⋅∂2F1∂x∂d3∂2F1∂w∂d3∂2F1∂y∂d3∂2F2∂x∂d3∂2F2∂w∂d3∂2F2∂y∂d3∂2F3∂x∂d3∂2F3∂w∂d3∂2F3∂w∂d3E2⋅v1v2v3=0,0,1⋅00000000−1E2⋅v1v2v3=−v3=−1≠0,Δ3=WT⋅D2Fr−d1a1,0,0;d3TCV,V=0,0,1⋅D∂F1∂xv1+∂F1∂wv2+∂F1∂yv3∂F2∂xv1+∂F2∂wv2+∂F2∂yv3∂F3∂xv1+∂F3∂wv2+∂F3∂yv3E2⋅v1v2v3=0,0,1⋅∂2F1∂x2v12+∂2F1∂w2v22+∂2F1∂y2v32+2∂2F1∂x∂wv1v2+2∂2F1∂x∂yv1v3+2∂2F1∂w∂yv2v3∂2F2∂x2v12+∂2F2∂w2v22+∂2F2∂y2v32+2∂2F2∂x∂wv1v2+2∂2F2∂x∂yv1v3+2∂2F2∂w∂yv2v3∂2F3∂x2v12+∂2F3∂w2v22+∂2F3∂y2v32+2∂2F3∂x∂wv1v2+2∂2F3∂x∂yv1v3+2∂2F3∂w∂yv2v3E2=2a2b1c1v1−m/ab1+r−d1/a12≠0.
Therefore, by Sotomayor’s theorem [43], system (2) exhibits a transcritical bifurcation at d3=d3TC around the axial equilibrium point E2.Theorem 12.
System (2) undergoes a transcritical bifurcation around E40,w^,y^ if(48)d1TC=r1+ky^−a2a+1−my^y^b1+y^a+y^,c2a3w^y^b1+w^a+y^<c2a3w^y^a+y^b1+w^b1+w^a+y^2.Proof.
Proof is the same as in Theorem10.Theorem 13.
System (2) undergoes a transcritical bifurcation around E5r−d1/a1,A/d2,0 if(49)d3TC=c1a2r−d1/a1b1+A/d2+r−d1/a1+c2a3A/d2b1+A/d2+r−d1/a1.Proof.
Proof is the same as in Theorem11.


### 4.3.2. Hopf Bifurcation aroundE∗x∗,w∗,y∗

Let us considerk as a bifurcation parameter of system (2) where the characteristic equation at E∗ is(50)λ3+A1kλ2+A2kλ+A3k=0.Then, Hopf bifurcation theorem is stated as follows.Theorem 14.
(Hopf bifurcation theorem [44]). If A1k, A2k, and A3k are the smooth functions of k in NεkH, ε>0, kH∈ℝ for which the characteristic equation (50) has the following:(i)
A pair of imaginary eigenvaluesλ=q1k±iq2k with q1k and q2k∈ℝ so that they become purely complex at k=kH and dq1/dk|k=kH≠0(ii)
The other eigenvalue is negative atk=kH; then, a Hopf bifurcation appears around E∗ at k=kHTheorem 15.
IfE∗x∗,w∗,y∗ is locally asymptotically stable, then a Hopf bifurcation is exhibited around E∗x∗,w∗,y∗ when k passes through its critical value kH provided A1kH>0, A3kH>0, and A1kHA2kH=A3kH (kH is a positive root of equation A1kA2k−A3k=0).Proof.
Atk=kH, we can write equation (50) as(51)λ2+A2λ+A1=0.
The roots of equation (51) are λ1=iA2, λ2=−iA2, and λ3=−A1. Also A1, A2, and A3 are the smooth functions of k. So, the roots of equation (59) have the form λ1=p1k+ip2k, λ2=p1k−ip2k, and λ3=p3k where pik are real functions of k in an open neighborhood of kH for i=1,2,3. Next, we verify the transversality condition:(52)ddkReλik|k=kH≠0,i=1,2.
Puttingλk=p1k+ip2k in (59), we get(53)p1k+ip2k3+A1kp1k+ip2k2+A2kp1k+ip2k+A3k=0.
Differentiating both sides with respect tok, we have(54)3p1k+ip2k2p˙1k+ip˙2k+2A1kp1k+ip2kp˙1k+ip˙2k+A˙1kp1k+ip2k2+A2p˙1k+ip˙2k+A˙2kp1k+ip2k+A˙3k=0.
Comparing real and imaginary parts from both sides, we obtain(55)X1p˙1−X2p˙2+X3=0,(56)X2p˙1+X1p˙2+X4=0,where(57)X1=3p12−p22+2A1p1+A2,X2=6p1p2+2A1p2,X3=A˙1p12−p22+A˙2p1+A˙3≠A˙1p12−p22+A˙2p1+A1A˙2+A2A˙1,sinceA3≠A1A2in a deleted neighbour hood ofkH,and X4=2A˙1p1p2+A˙2p2.
Multiplying (55) by X1 and (56) by X2 and then adding, we get(58)X12+X22p˙1+X1X3+X2X4=0⟹p˙1=−X1X3+X2X4X12+X22.
Atk=kH,Case 1:p1=0, p2=A2. Then, X1=−2A2, X2=2A1A2, X3≠A1A˙2, and X4=A˙2A2. ∴X1X3+X2X4≠0.Case 2:p1=0, p2=−A2. Then, X1=−2A2, X2=−2A1A2, X3≠A1A˙2, and X4=−A˙2A2. ∴X1X3+X2X4≠0.
Also,λ3=−A1kH<0.
Hence, this theorem is proved by virtue of Theorem14.


### 4.3.3. Hopf Bifurcation aroundE3x¯,0,y¯

Let us considerk as a bifurcation parameter of system (2) where the characteristic equation of E3 is(59)λ3+B1kλ2+B2kλ+B3k=0,and then Hopf bifurcation theorem is stated as follows.Theorem 16.
(Hopf bifurcation theorem [44]). If B1k, B2k, and B3k are the smooth functions of k in Nεk∗, ε>0, k∗∈ℝ for which the characteristic equation (59) has the following:(i)
A pair of imaginary eigenvaluesλ=p1′k±ip2′k with p1′k and p2′k∈ℝ so that they become purely imaginary at k=k∗ and dp1′/dk|k=k∗≠0(ii)
The other eigenvalue is negative atk=k∗; then a Hopf bifurcation occurs around E3x¯,0,y¯ at k=k∗Theorem 17.
IfE3x¯,0,y¯ is locally asymptotically stable, then a Hopf bifurcation appears around subsidy-free equilibrium E3x¯,0,y¯ when k passes through its critical value k∗ provided B1k∗>0, B3k∗>0, and B1k∗B2k∗=B3k∗ (k∗ is a positive root of equation B1kB2k−B3k=0).Proof.
Proof is the same as in Theorem15.


## 4.3.1. Transcritical Bifurcation

Theorem 10.
System (2) undergoes a transcritical bifurcation around E10,A/d2,0 if d1TC=r and c2a3A/d2/b1+A/d2<d3 (TC stands for transcritical bifurcation).Proof.
We apply Sotomayor’s theorem [43] to prove the occurrence of a transcritical bifurcation around E1 with d1 as bifurcation parameter. For applicability of Sotomayor’s theorem, exactly one of the eigenvalues of the Jacobian matrix at E1 must be zero and other eigenvalues must have negative real parts. So, we need to fulfill the condition c2a3A/d2/b1+A/d2<d3.
The eigenvectors ofJE1=vij and JE1T corresponding to the zero eigenvalue of E10,A/d2,0 are obtained as V=v1,v2,v3T and W=1,0,0T, respectively, where v1=1, v2=v3=0, v11=0, v12=0, v13=0, v21=0, v22=−d2, v23=0, v31=0, v32=0, v33=−d3+c2a3A/d2/b1+A/d2.
ComputeΔ1, Δ2, and Δ3 as follows:(43)Δ1=WT⋅Fd10,Ad2,0;d1TC=1,0,0⋅∂F1∂d1∂F2∂d1∂F3∂d1E1=1,0,0⋅−x00E1=0,where F=F1,F2,F3T and F1, F2, and F3 are given by(44)F1=rx1+ky−d1x−a1x2−a21−my/a+yxyb1+w+1−my/a+yx,F2=A−d2w−a3wyb1+w+1−my/a+yx,F3=c1a21−my/a+yxyb1+w+1−my/a+yx+c2a3wyb1+w+1−my/a+yx−d3y,Δ2=WT⋅DFd10,Ad2,0;d1TCV=1,0,0⋅∂2F1∂x∂d1∂2F1∂w∂d1∂2F1∂y∂d1∂2F2∂x∂d1∂2F2∂w∂d1∂2F2∂y∂d1∂2F3∂x∂d1∂2F3∂w∂d1∂2F3∂y∂d1E1⋅100=1,0,0⋅−100000000E1⋅100=−1≠0,Δ3=WT⋅D2F0,Ad2,0;d1TCV,V=1,0,0⋅D∂F1∂xv1+∂F1∂wv2+∂F1∂yv3∂F2∂xv1+∂F2∂wv2+∂F2∂yv3∂F3∂xv1+∂F3∂wv2+∂F3∂yv3E1⋅v1v2v3=1,0,0⋅∂2F1∂x2v12+∂2F1∂w2v22+∂2F1∂y2v32+2∂2F1∂x∂wv1v2+2∂2F1∂x∂yv1v3+2∂2F1∂w∂yv2v3∂2F2∂x2v12+∂2F2∂w2v22+∂2F2∂y2v32+2∂2F2∂x∂wv1v2+2∂2F2∂x∂yv1v3+2∂2F2∂w∂yv2v3∂2F3∂x2v12+∂2F3∂w2v22+∂2F3∂y2v32+2∂2F3∂x∂wv1v2+2∂2F3∂x∂yv1v3+2∂2F3∂w∂yv2v3E1=−2a1≠0.
Therefore, by Sotomayor’s theorem [43], system (2) undergoes a transcritical bifurcation at d1=d1TC around the axial equilibrium point E1.Theorem 11.
System (2) exhibits a transcritical bifurcation around E2r−d1/a1,0,0 if(45)d3TC=c1a2r−d1/a1b1+r−d1/a1.Proof.
Let us apply Sotomayor’s theorem [43] to prove the occurrence of a transcritical bifurcation around E2 with d3 as bifurcation parameter. For applicability of Sotomayor’s theorem, exactly one of the eigenvalues of the Jacobian matrix at E2 must be zero and other eigenvalues must have negative real parts.
The eigenvectors ofJE2=tij and JE2T corresponding to the zero eigenvalue of E2r−d1/a1,0,0 are obtained as V=v1,v2,v3T and W=0,0,1T, respectively, where v1=1/a1−rk−a2/b1+r−d1/a1<0, v2=0, v3=1, t11=−r+d1, t12=0, t13=r−d1/a1−rk−a2/b1+r−d1/a1, t21=0, t22=−d2, t23=0, t31=0, t32=0, t33=0.
ComputeΔ1, Δ2, and Δ3 as follows:(46)Δ1=WT⋅Fd1r−d1a1,0,0;d3TC=0,0,1⋅∂F1∂d3∂F2∂d3∂F3∂d3E2=0,0,1⋅00−yE2=0,where F=F1,F2,F3T and F1, F2, and F3 are given by(47)F1=rx1+ky−d1x−a1x2−a21−my/a+yxyb1+w+1−my/a+yx,F2=A−d2w−a3wyb1+w+1−my/a+yx,F3=c1a21−my/a+yxyb1+w+1−my/a+yx+c2a3wyb1+w+1−my/a+yx−d3y,Δ2=WT⋅DFd3r−d1a1,0,0;d3TCV=0,0,1⋅∂2F1∂x∂d3∂2F1∂w∂d3∂2F1∂y∂d3∂2F2∂x∂d3∂2F2∂w∂d3∂2F2∂y∂d3∂2F3∂x∂d3∂2F3∂w∂d3∂2F3∂w∂d3E2⋅v1v2v3=0,0,1⋅00000000−1E2⋅v1v2v3=−v3=−1≠0,Δ3=WT⋅D2Fr−d1a1,0,0;d3TCV,V=0,0,1⋅D∂F1∂xv1+∂F1∂wv2+∂F1∂yv3∂F2∂xv1+∂F2∂wv2+∂F2∂yv3∂F3∂xv1+∂F3∂wv2+∂F3∂yv3E2⋅v1v2v3=0,0,1⋅∂2F1∂x2v12+∂2F1∂w2v22+∂2F1∂y2v32+2∂2F1∂x∂wv1v2+2∂2F1∂x∂yv1v3+2∂2F1∂w∂yv2v3∂2F2∂x2v12+∂2F2∂w2v22+∂2F2∂y2v32+2∂2F2∂x∂wv1v2+2∂2F2∂x∂yv1v3+2∂2F2∂w∂yv2v3∂2F3∂x2v12+∂2F3∂w2v22+∂2F3∂y2v32+2∂2F3∂x∂wv1v2+2∂2F3∂x∂yv1v3+2∂2F3∂w∂yv2v3E2=2a2b1c1v1−m/ab1+r−d1/a12≠0.
Therefore, by Sotomayor’s theorem [43], system (2) exhibits a transcritical bifurcation at d3=d3TC around the axial equilibrium point E2.Theorem 12.
System (2) undergoes a transcritical bifurcation around E40,w^,y^ if(48)d1TC=r1+ky^−a2a+1−my^y^b1+y^a+y^,c2a3w^y^b1+w^a+y^<c2a3w^y^a+y^b1+w^b1+w^a+y^2.Proof.
Proof is the same as in Theorem10.Theorem 13.
System (2) undergoes a transcritical bifurcation around E5r−d1/a1,A/d2,0 if(49)d3TC=c1a2r−d1/a1b1+A/d2+r−d1/a1+c2a3A/d2b1+A/d2+r−d1/a1.Proof.
Proof is the same as in Theorem11.


## 4.3.2. Hopf Bifurcation aroundE∗x∗,w∗,y∗

Let us considerk as a bifurcation parameter of system (2) where the characteristic equation at E∗ is(50)λ3+A1kλ2+A2kλ+A3k=0.Then, Hopf bifurcation theorem is stated as follows.Theorem 14.
(Hopf bifurcation theorem [44]). If A1k, A2k, and A3k are the smooth functions of k in NεkH, ε>0, kH∈ℝ for which the characteristic equation (50) has the following:(i)
A pair of imaginary eigenvaluesλ=q1k±iq2k with q1k and q2k∈ℝ so that they become purely complex at k=kH and dq1/dk|k=kH≠0(ii)
The other eigenvalue is negative atk=kH; then, a Hopf bifurcation appears around E∗ at k=kHTheorem 15.
IfE∗x∗,w∗,y∗ is locally asymptotically stable, then a Hopf bifurcation is exhibited around E∗x∗,w∗,y∗ when k passes through its critical value kH provided A1kH>0, A3kH>0, and A1kHA2kH=A3kH (kH is a positive root of equation A1kA2k−A3k=0).Proof.
Atk=kH, we can write equation (50) as(51)λ2+A2λ+A1=0.
The roots of equation (51) are λ1=iA2, λ2=−iA2, and λ3=−A1. Also A1, A2, and A3 are the smooth functions of k. So, the roots of equation (59) have the form λ1=p1k+ip2k, λ2=p1k−ip2k, and λ3=p3k where pik are real functions of k in an open neighborhood of kH for i=1,2,3. Next, we verify the transversality condition:(52)ddkReλik|k=kH≠0,i=1,2.
Puttingλk=p1k+ip2k in (59), we get(53)p1k+ip2k3+A1kp1k+ip2k2+A2kp1k+ip2k+A3k=0.
Differentiating both sides with respect tok, we have(54)3p1k+ip2k2p˙1k+ip˙2k+2A1kp1k+ip2kp˙1k+ip˙2k+A˙1kp1k+ip2k2+A2p˙1k+ip˙2k+A˙2kp1k+ip2k+A˙3k=0.
Comparing real and imaginary parts from both sides, we obtain(55)X1p˙1−X2p˙2+X3=0,(56)X2p˙1+X1p˙2+X4=0,where(57)X1=3p12−p22+2A1p1+A2,X2=6p1p2+2A1p2,X3=A˙1p12−p22+A˙2p1+A˙3≠A˙1p12−p22+A˙2p1+A1A˙2+A2A˙1,sinceA3≠A1A2in a deleted neighbour hood ofkH,and X4=2A˙1p1p2+A˙2p2.
Multiplying (55) by X1 and (56) by X2 and then adding, we get(58)X12+X22p˙1+X1X3+X2X4=0⟹p˙1=−X1X3+X2X4X12+X22.
Atk=kH,Case 1:p1=0, p2=A2. Then, X1=−2A2, X2=2A1A2, X3≠A1A˙2, and X4=A˙2A2. ∴X1X3+X2X4≠0.Case 2:p1=0, p2=−A2. Then, X1=−2A2, X2=−2A1A2, X3≠A1A˙2, and X4=−A˙2A2. ∴X1X3+X2X4≠0.
Also,λ3=−A1kH<0.
Hence, this theorem is proved by virtue of Theorem14.


## 4.3.3. Hopf Bifurcation aroundE3x¯,0,y¯

Let us considerk as a bifurcation parameter of system (2) where the characteristic equation of E3 is(59)λ3+B1kλ2+B2kλ+B3k=0,and then Hopf bifurcation theorem is stated as follows.Theorem 16.
(Hopf bifurcation theorem [44]). If B1k, B2k, and B3k are the smooth functions of k in Nεk∗, ε>0, k∗∈ℝ for which the characteristic equation (59) has the following:(i)
A pair of imaginary eigenvaluesλ=p1′k±ip2′k with p1′k and p2′k∈ℝ so that they become purely imaginary at k=k∗ and dp1′/dk|k=k∗≠0(ii)
The other eigenvalue is negative atk=k∗; then a Hopf bifurcation occurs around E3x¯,0,y¯ at k=k∗Theorem 17.
IfE3x¯,0,y¯ is locally asymptotically stable, then a Hopf bifurcation appears around subsidy-free equilibrium E3x¯,0,y¯ when k passes through its critical value k∗ provided B1k∗>0, B3k∗>0, and B1k∗B2k∗=B3k∗ (k∗ is a positive root of equation B1kB2k−B3k=0).Proof.
Proof is the same as in Theorem15.


## 5. Delayed Dynamical System

In biological point of view, many processes, both natural and man-made, include time-delay. The study of delay factor makes our system much more realistic than non-delayed system. Also, a delay differential equation reveals much more complicated dynamics than an ordinary differential equation (for details, see [10–13, 23, 45–49]).In reality, after sensing the vocal cue, individuals of prey species take some time for assessing the predation risk. So, the effect of fear (felt by prey) of predator does not respond spontaneously on the birth rate of prey population; some time lag must be needed. In view of this fact, the predator-prey-subsidy interactions (2) can be modified as follows:(60)dxdt=rx1+kyt−τ−d1x−a1x2−a21−my/a+yxyb1+w+1−my/a+yx,dwdt=A−d2w−a3wyb1+w+1−my/a+yx,dydt=c1a21−my/a+yxyb1+w+1−my/a+yx+c2a3wyb1+w+1−my/a+yx−d3y.The initial conditions are assumed as (i=1,2,3)(60a)ψiϕ>0,ϕ∈−τ,0,wherexϕ=ψ1ϕ,wϕ=ψ2ϕ,yϕ=ψ3ϕ,For biological feasibility:ψ10>0,ψ20>0,ψ30>0.Let us linearize (60) using the following transformations:(61)X=x−x¯,W=w−0,Y=y−y¯.It leads to(62)dUdt=B1′Ut+B2′Ut−τ,where U=X,W,YT,(63)B1′=b11b12b13′b21b22b23b31b32b33,B2′=00b13′′000000,b11=x¯−a1+a2y¯a+1−my¯2b1a+y¯+a+1−my¯x¯2,b12=a2x¯y¯a+1−my¯a+y¯b1a+y¯+a+1−my¯x¯2,b13′=x¯−a2a+21−my¯b1a+y¯+a+1−my¯x¯+a2y¯a+1−my¯b1+1−mx¯b1a+y¯+a+1−my¯x¯2,b21=0,b22=−d2−a3y¯b1+1−my¯/a+y¯x¯,b23=0,b31=y¯c1a2a+1−my¯b1a+y¯+a+1−my¯x¯−c1a2a+1−my¯2x¯b1a+y¯+a+1−my¯x¯2,b32=c2a3y¯a+y¯b1a+y∗+a+1−my¯x¯,b33=y¯c1a21−mx¯b1a+y¯+a+1−my¯x¯−c1a2a+1−my¯x¯b1+1−mx¯b1a+y¯+a+1−my¯x¯2,b13″=−rkx¯1+ky¯2.The characteristic equation corresponding to (62) is(64)λ3+L1λ2+L2λ+L3+M1λ+M2e−λτ=0,where(65)L1=−b11+b22+b33,L2=b22b33+b11b33−b13′b31+b11b22,L3=−b11b22b33−b13′b22b31,M1=−b13″b31,M2=b13″b31b22.Ifτ≠0, E3 of system (60) is LAS provided equation (64) has no purely imaginary roots and it is also LAS for τ=0. Further, it has been shown that stability nature of E3 switches at τ=τ′∗. Already, it has been derived that E3 is LAS provided B1>0, B3>0, and B1B2>B3 for τ=0 (non-delayed system). Let us discuss if the real part of the roots of equation (64) gradually increases to reach zero and eventually turns to a positive value when τ increases.Substitutingλ=q1′+iq2′ in equation (64), we have(66)q1′+iq2′3+L1q1′+iq2′2+L2q1′+iq2′+L3+M1q1′+iq2′+M2e−q1′τcosq2′τ−isinq2′τ=0.Equating respective real and complex parts from both sides, we obtain(67)q1′3−3q1′q2′2+L1q1′2−q2′2+L2q1′+L3+M1q1′e−q1′τcosq2′τ+M2e−q1′τcosq2′τ+M1q2′e−q1′τsinq2′τ=0,(68)3q1′2q2′−q2′3+2L1q1′q2′+L2q2′+M1q2′e−q1′τcosq2′τ−M1q1′e−q1′τsinq2′τ−M2e−q1′τsinq2′τ=0.Now, let us examine whether equation (64) has purely imaginary roots or not. For this purpose, let us take q1′=0. Then, equations (67) and (68) become(69)M1q2′sinq2′τ+M2cosq2′τ=L1q2′2−L3,(70)M1q2′cosq2′τ−M2sinq2′τ=q2′3−L2q2′.Eliminatingτ from (69) and (70) (squaring and adding), we get(71)q2′6+q2′4L12−2L2+q2′2L22−2L1L3−M12+L32−M22=0.Puttingq2′2=β, we have(72)Lβ≡β3+L12−2L2β2+L22−2L1L3−M12β+L32−M22=0.This is a cubic equation ofβ. It is noticed that L∞=∞. So, equation (72) has exactly one positive real root if L0<0, i.e., if L32<M22.Letβ=β+ be a positive root of (77); then, q2′=β+.Lemma 1.
(see [50]). Consider the exponential polynomial:(73)Pλ≡Pλ,τ1,τ2,…,τm≡λn+p10λn−1+⋯+pn−10λ+pn0+p11λn−1+⋯+pn−11λ+pn1e−λτ1+⋯+p1mλn−1+⋯+pn−1mλ+pnme−λτm,where τi≥0i=1,2,…,m and pjii=0,1,…,m;j=1,2,…,n are constants. As τ1,τ2,…,τm vary, the sum of the orders of zero of Pλ in the open half plane can change only if a zero appears on or crosses the imaginary axis.Now, let us discuss the existence of Hopf bifurcation aroundE3 with τ as a bifurcation parameter.Theorem 18.
SupposeE3 exists and is locally asymptotically stable for system (2) when τ=0. If L32<M22, then there exists a critical value τ′∗ such that E3 of system (60) is LAS when τ∈0,τ′∗ and unstable when τ>τ′∗, where(74)τ+′j=cos−1M2L1β+−L3+M1β+2−L2β+/M22+M12β+β++2πjβ+,j=0,1,2,3,…,and τ′∗=τ+′0 (minimum value). Also, system (60) exhibits Hopf bifurcation around E3 at τ=τ′∗ provided K′M′−L′N′≠0, where(75)K′=−3β++L2+M1cosβ+τ′∗−M2τ′∗cosβ+τ′∗−M1β+τ′∗sinβ+τ′∗,L′=−2L1β++M1sinβ+τ′∗−M2τ′∗sinβ+τ′∗+M1β+τ′∗cosβ+τ′∗,M′=M2β+sinβ+τ′∗−M1β+cosβ+τ′∗,N′=M2β+cosβ+τ′∗+M1β+sinβ+τ′∗.Proof.
IfL32<M22, then (72) has exactly one positive root β+, i.e., from (69) and (70), τ+′j,j=0,1,2,…, are obtained as functions of β+:(76)τ+′j=cos−1M2L1β+−L3+M1β+2−L2β+/M22+M12β+β++2πjβ+,j=0,1,2,3,….
IfE3 is locally asymptotically stable, the stability behaviour of E3 will remain unaltered for τ<τ′∗ (using Butler’s Lemma [51]).
To check the transversality condition,d/dτReλττ=τ′∗≠0, let us differentiate (67) and (68) with respect to τ and set q1′=0 and τ=τ′∗. The following equations are obtained:(77)K′ddτReλττ=τ′∗+L′ddτImλττ=τ′∗=M′,(78)−L′ddτReλττ=τ′∗+K′ddτImλττ=τ′∗=N′,where(79)K′=−3β++L2+M1cosβ+τ′∗−M2τ′∗cosβ+τ′∗−M1β+τ′∗sinβ+τ′∗,L′=−2L1β++M1sinβ+τ′∗−M2τ′∗sinβ+τ′∗+M1β+τ′∗cosβ+τ′∗,M′=M2β+sinβ+τ′∗−M1β+cosβ+τ′∗,N′=M2β+cosβ+τ′∗+M1β+sinβ+τ′∗.
Solving (77) and (78),(80)dReλτdττ=τ′∗=K′M′−L′N′K′2+L′2.
Now, we havedReλτ/dττ=τ′∗≠0, if K′M′−L′N′≠0. Hence, the transversality condition is satisfied and a Hopf bifurcation occurs around E∗ when τ passes through its critical value τ′∗.Now, linearize system (60) using the transformations X=x−x∗, W=w−w∗, and Y=y−y∗:(81)dUdt=A1′Ut+A2′Ut−τ,where U=X,W,YT,(82)A1′=a11a12a13′a21a22a23a31a32a33,A2′=00a13″000000,a11=x∗−a1+a2a+1−my∗2y∗b1+w∗a+y∗+a+1−my∗x∗2,a12=a2a+1−my∗a+y∗x∗y∗b1+w∗a+y∗+a+1−my∗x∗2,a13′=−a2a+21−my∗x∗b1+w∗a+y∗+a+1−my∗x∗+a2a+1−my∗x∗y∗b1+w∗+1−mx∗b1+w∗a+y∗+a+1−my∗x∗2,a21=a3a+y∗a+1−my∗w∗y∗b1+w∗a+y∗+a+1−my∗x∗2,a22=−d2−a3a+y∗y∗b1+w∗a+y∗+a+1−my∗x∗+a3a+y∗2w∗y∗b1+w∗a+y∗+a+1−my∗x∗2,a23=−a3a+2y∗w∗b1+w∗a+y∗+a+1−my∗x∗+a3a+y∗b1+w∗+1−mx∗b1+w∗a+y∗+a+1−my∗x∗2,a31=c1a2a+1−my∗y∗b1+w∗a+y∗+a+1−my∗x∗−c1a2a+1−my∗2x∗y∗b1+w∗a+y∗+a+1−my∗x∗2−c2a3w∗y∗a+y∗a+1−my∗b1+w∗a+y∗+a+1−my∗x∗2a32=−c1a2a+1−my∗x∗y∗a+y∗b1+w∗a+y∗+a+1−my∗x∗2+c2a3a+y∗y∗b1+w∗a+y∗+a+1−my∗x∗−c2a3w∗y∗a+y∗2b1+w∗a+y∗+a+1−my∗x∗2,a33=c1a21−mx∗y∗+c2a3w∗y∗b1+w∗a+y∗+a+1−my∗x∗−c1a2a+1−my∗x∗y∗+c2a3w∗y∗a+y∗b1+w∗+1−mx∗b1+w∗a+y∗+a+1−my∗x∗2,a13′′=−rkx∗1+ky∗2.The characteristic equation corresponding to (81) is(83)λ3+R1λ2+R2λ+R3+S1λ+S2e−λτ=0,where(84)R1=−a11+a22+a33,R2=a22a33−a23a32+a11a33−a13′a31+a11a22−a12a21,R3=−a11a22a33−a23a32+a12a23a31−a21a33+a13′a21a32−a22a31,S1=−a13″a31,S2=−a13″a21a32−a31a22.Ifτ≠0, E∗ of system (60) is LAS provided equation (83) has no purely imaginary roots and it is LAS for τ=0. Furthermore, it has to be noted that changes of stability occur at τ=τ∗. Already, it has been discussed that E∗ is LAS when τ=0 provided A1>0, A3>0, and A1A2>A3. Here, equation (83) is a transcendental equation, so it contains infinitely many eigenvalues. In this situation, we cannot apply the Routh–Hurwitz criteria to determine the stability of system (60). To understand the stability behaviour, our necessity is to check the sign of the real parts of the eigenvalues of equation (83).Now, puttingλ=q1+iq2 in equation (83), we have(85)q1+iq23+R1q1+iq22+R2q1+iq2+R3+S1q1+iq2+S2e−q1τcosq2τ−isinq2τ=0.Equating respective real and complex parts from both sides, we get(86)q13−3q1q22+R1q12−q22+R2q1+R3+S1q1e−q1τcosq2τ+S2e−q1τcosq2τ+S1q2e−q1τsinq2τ=0,(87)3q12q2−q23+2R1q1q2+R2q2+S1q2e−q1τcosq2τ−S1q1e−q1τsinq2τ−S2e−q1τsinq2τ=0.To check whether (83) has purely imaginary roots or not, set q1=0; then, (86) and (87) become(88)S1q2sinq2τ+S2cosq2τ=R1q22−R3,(89)S1q2cosq2τ−S2sinq2τ=q23−R2q2.Eliminatingτ from (88) and (89) (squaring and adding), we get(90)q26+q24R12−2R2+q22R22−2R1R3−S12+R32−S22=0.Puttingq22=σ, we have(91)R′σ≡σ3+R12−2R2σ2+R22−2R1R3−S12σ+R32−S22=0.This is a cubic equation ofσ. It is noted that R′∞=∞. So, equation (91) has exactly one positive root if R′0<0, i.e., if R32<S22.Letσ=σ+ be a positive root of (91); then, q2=σ+.Let us study the existence of Hopf bifurcation aroundE∗ with τ as bifurcation parameter.Theorem 19.
SupposeE∗ exists and is locally asymptotically stable for system (2) when τ=0. If R32<S22, then there exists a critical value τ∗ such that E∗ of system (60) is LAS when τ∈0,τ∗ and unstable when τ>τ∗, where(92)τ+j=cos−1S2R1σ+−R3+S1σ+2−R2σ+/S22+S12σ+σ++2πjσ+,j=0,1,2,3,…,and τ∗=τ+0 (minimum value). Also, a supercritical Hopf bifurcation is exhibited around E∗ at τ=τ∗ provided K″M″−L″N″≠0, where(93)K″=−3σ++R2+S1cosσ+τ∗−S2τ∗cosσ+τ∗−S1σ+τ∗sinσ+τ∗,L″=−2R1σ++S1sinσ+τ∗−S2τ∗sinσ+τ∗+S1σ+τ∗cosσ+τ∗,M″=S2σ+sinσ+τ∗−S1σ+cosσ+τ∗,N″=S2σ+cosσ+τ∗+S1σ+sinσ+τ∗.Proof.
Proof is similar to that in Theorem18.Theorem 20.
Suppose interior (coexistence) equilibrium pointE∗ exists and is locally asymptotically stable for system (2) when τ=0. Let equation (91) have exactly two positive roots σi,i=1,2σ1>σ2 when R32>S22 and R22−2R1R3−S12<0 irrespective of sign of R12−2R2. Moreover, let(94)τji=cos−1S2R1σi−R3+S1σi2−R2σi/S22+S12σiσi+2πjσi,i=1,2;j=0,1,2,3…,τk+=minτki:i=1,2,τk−=maxτki:i=1,2,k=0,1,2,3…,then there is a positive integer k such that(95)0<τ0+<τ0−<τ1+<τ1−<τ2+<⋯<τk−1−<τk+,and there are k switches from stability to instability to stability; that is, when(96)τ∈0,τ0+,τ0−,τ1+,τ1−,τ2+,…,τk−1−,τk+,then E∗ is locally asymptotically stable and when(97)τ∈τ0+,τ0−,τ1+,τ1−,τ2+,τ2−,…,τk−1+,τk−1−,τ>τk+,then E∗ is unstable. Further, at τ=τk±,k=0,1,2,…, system (60) experiences Hopf bifurcation provided(98)Pk+Rk+−Qk+Sk+≠0,fork=0,1,2,…,Pk−Rk−−Qk−Sk−≠0,fork=0,1,2,…,where(99)Pk+=−3σ1+R2+S1cosσ1τk+−S2τk+cosσ1τk+−S1σ1τk+sinσ1τk+k=0,1,2,…,Pk−=−3σ2+R2+S1cosσ2τk−−S2τk−cosσ2τk−−S1σ2τk−sinσ2τk−k=0,1,2,…,Qk+=−2R1σ1+S1sinσ1τk+−S2τk+sinσ1τk++S1σ1τk+cosσ1τk+k=0,1,2,…,Qk−=−2R1σ2+S1sinσ2τk−−S2τk−sinσ2τk−+S1σ2τk−cosσ2τk−k=0,1,2,…,Rk+=S2σ1sinσ1τk+−S1σ1cosσ1τk+k=0,1,2,…,Rk−=S2σ2sinσ2τk−−S1σ2cosσ2τk−k=0,1,2,…,Sk+=S2σ1cosσ1τk++S1σ1sinσ1τk+k=0,1,2,…,Sk−=S2σ2cosσ2τk−+S1σ2sinσ2τk−k=0,1,2,….Proof.
Proof is similar to that in Theorem18.


## 6. Numerical Computations

Here, we have illustrated numerical simulations to verify the analytical findings of the proposed system (2). We select a parameter set: {r=5.5, d1=6.5, d2=4, k=0.2, d3=5, a1=2, a2=0.3, a3=0.25, b1=1.5, c1=0.7, c2=0.5, A=0, a=1.2, m=0.01}. Under this set of parametric values, the stable nature of E00,0,0 is shown in Figure 1. If we take subsidy input rate A=2 and other parametric values are chosen from the data set of Figure 1, then the subsidy only equilibrium E10,A/d2,0≡E10,0.5,0 exists and stable nature of E10,0.5,0 with time t is depicted in Figure 2. Now, we choose another parameter set: {r=5.5, d1=1.5, d2=0.52, k=1.9, d3=2, a1=0.5, a2=0.3, a3=0.25, b1=1.5, c1=0.7, c2=0.5, A=0, a=1.2, m=0.01}. Under this set of parametric values, the prey only equilibrium E2r−d1/a1,0,0≡E28,0,0 exists and stable behaviour of E28,0,0 is presented in Figure 3. Let us choose the parameters as follows:(100)r=5.5,d1=0.4,d2=0.3,k=0.2,d3=0.2,a1=0.6,a2=0.98,a3=0.8,b1=2.5,c1=0.85,c2=0.7,a=1.1,m=0.4.Figure 1
Stable behaviour ofE00,0,0 with respect to time t corresponding to the data set {r=5.5, d1=6.5, d2=4, k=0.2, d3=5, a1=2, a2=0.3, a3=0.25, b1=1.5, c1=0.7, c2=0.5, A=0, a=1.2, m=0.01}.Figure 2
Stable behaviour ofE10,0.5,0 with respect to time t taking the parameters as {r=5.5, d1=6.5, d2=4, k=0.2, d3=5, a1=2, a2=0.3, a3=0.25, b1=1.5, c1=0.7, c2=0.5, A=2, a=1.2, m=0.01}.Figure 3
Stable nature ofE28,0,0 with respect to time t regarding the parameters as {r=5.5, d1=1.5, d2=0.52, k=1.9, d3=2, a1=0.5, a2=0.3, a3=0.25, b1=1.5, c1=0.7, c2=0.5, A=0, a=1.2, m=0.01}.If we take subsidy input rateA=0 and other parameters are taken from set (100), then subsidy-free equilibrium point E3x¯,0,y¯≡E31.20053,0,6.49267 exists and is locally asymptotically stable. Stable time series and stable phase diagram are represented in Figure 4. In the same manner, if we change the value of the parameter k=0.01 and others are the same as in the data set of Figure 4, then it is observed that E31.20053,0,6.49267 is unstable accompanied with a limit cycle (see Figure 5). From Figures 4 and 5, it can be easily noted that there must exist a threshold value of k, say k∗=0.019 for which unstable behaviour of E3 changes to stable spiral. Since the vector fields for k<k∗ and k>k∗ are qualitatively different, a Hopf bifurcation is created around E3 taking k as bifurcation parameter (see Figure 6). For the set of parameter values {r=5.5, d1=6.5, d2=0.3, k=0.2, d3=0.2, a1=0.6, a2=0.98, a3=0.8, b1=2.5, c1=0.85, c2=0.7, A=2, a=1.1, m=0.01}, prey free equilibrium point E40,w^,y^≡E40,1.3889,5.5417 exists and is stable (see Figure 7). Next, let us take a different set of parameters of system (2): {r=5.5, d1=0.5, d2=0.3, k=0.2, d3=2.2, a1=0.6, a2=0.98, a3=0.8, b1=2.5, c1=0.85, c2=0.7, A=2, a=1.1, m=0.01}. Then, predator free equilibrium point E5r−d1/a1,A/d2,0≡E58.3333,6.6667,0 is locally asymptotically stable. The stable behaviour with time t is shown in Figure 8.Figure 4
(a) Stable nature ofE31.20053,0,6.49267 with t and (b) stable phase portrait of E31.20053,0,6.49267 when k=0.2>k∗ (threshold value) =0.019 and other are taken as {r=5.5, d1=0.4, d2=0.3, d3=0.2, a1=0.6, a2=0.98, a3=0.8, b1=2.5, c1=0.85, c2=0.7, A=0, a=1.1, m=0.4}.
(a)(b)Figure 5
(a) Oscillatory behaviour ofE3 with time t and (b) phase diagram (isolated closed orbit) when k=0.01<k∗=0.019 and all other parameters are fixed as in Figure 4.
(a)(b)Figure 6
Bifurcation diagrams for the Hopf bifurcation aroundE31.20053,0,6.49267 regarding k as bifurcation parameter and others are the same as in Figure 4. (a) Bifurcation diagram of x. (b) Bifurcation diagram of y.
(a)(b)Figure 7
Stability nature ofE40,1.3889,5.5417 with time t regarding the parameters as {r=5.5, d1=6.5, d2=0.3, k=0.2, d3=0.2, a1=0.6, a2=0.98, a3=0.8, b1=2.5, c1=0.85, c2=0.7, A=2, a=1.1, m=0.01}.Figure 8
Stability nature ofE58.3333,6.6667,0 with time t regarding the parameters as {r=5.5, d1=0.5, d2=0.3, k=0.2, d3=2.2, a1=0.6, a2=0.98, a3=0.8, b1=2.5, c1=0.85, c2=0.7, A=2, a=1.1, m=0.01}.If we take subsidy input rateA=2 and others are fixed as in the data set of Figure 4, then E∗x∗,w∗,y∗≡E∗0.385717,0.950363,8.34509 exists and is locally asymptotically stable. Figure 9 depicts the stable behaviour of E∗. Comparing Figures 4 and 9, it is observed that subsidy input rate A enhances the value of y component of E∗ and decreases the value of x component of E∗. Also, from Figure 10, it is noticed that the prey population is leading towards extinction and the predator population has enormous growth (due to huge supply of food source) at high subsidy input rate (when A∈2,10) in the presence of fear felt by prey population. So, it can be concluded that it is not possible to control prey population from extinction in presence of nonlinear prey refuge because they cannot get enough time to protect themselves from predation risk. After extinction of prey species, predator can easily survive with the help of resource subsidy. Thus, the parameter A has great importance in the proposed population dynamics.Figure 9
(a) Stable nature ofE∗0.385717,0.950363,8.34509 with time t and (b) stable phase diagram of E∗0.385717,0.950363,8.34509 when k=0.2>kH=0.025 and others are chosen as {r=5.5, d1=0.4, d2=0.3, d3=0.2, a1=0.6, a2=0.98, a3=0.8, b1=2.5, c1=0.85, c2=0.7, A=2, a=1.1, m=0.4}.
(a)(b)Figure 10
Nature of steady stateE∗ when subsidy input rate A varies from 1 to 10 and other parameters are fixed as in Figure 9.Moreover, Figure11 represents the unstable nature of E∗ when k=0.02 and other parametric values are the same as in Figure 9. So, the parameter k has an interesting nature because there exists a threshold value kH=0.025 of k for which unstable nature (limit cycle) of E∗ switches to stable behaviour (stable spiral) when k passes through its critical value kH; i.e., the vector fields for k>kH and k<kH are topologically different. Hence, a Hopf bifurcation occurs around E∗ and Figure 12 depicts the corresponding bifurcation diagram taking k as bifurcation parameter. Also, it has to be noted from Figure 13 that, in the absence of fear effect, the oscillatory behaviour of E∗ changes to stable nature when subsidy input rate A crosses its critical value A∗=7.9 (approximately) and since predator population has huge growth rate at very large value of subsidy input rate A, the prey population cannot persist in ecosystem in presence nonlinear prey refuge. This phenomenon is very interesting because the prey refuge cannot control the prey population from extinction due to enormous growth of predator when subsidy input rate is very high. In this manner, system (2) is not persistent.Figure 11
(a), (b) Oscillatory behaviour ofE∗ with time and (c) isolated closed orbit around E∗ when k=0.02<kH=0.025 and all others are fixed as in Figure 9.
(a)(b)(c)Figure 12
Bifurcation diagrams for the Hopf bifurcation aroundE∗0.385717,0.950363,8.34509 regarding k as bifurcation parameter and others are the same as in Figure 9. (a) Bifurcation diagram of x. (b) Bifurcation diagram of w. (c) Bifurcation diagram of y.
(a)(b)(c)Figure 13
Bifurcation diagram with respect to subsidy input rateA when k=0 and the remaining parameters are fixed as in Figure 9. (a) Bifurcation diagram of x. (b) Bifurcation diagram of w. (c) Bifurcation diagram of y.
(a)(b)(c)Further, Figure14 depicts that the nature of steady states E3 and E∗ when m∈0,1. Here, the predator population cannot go extinct for large value of coefficient of prey refuge parameter. Also, Figure 15 shows the changes of predator’s growth at the steady states E3 and E∗ for three different fear levels k when m varies from 0 to 1. Here, also the predator can persist for large m. From here, it may be concluded that system (2) is always persistent for small subsidy input rate in the presence of nonlinear prey refuge function. Again, Figures 16 and 17, respectively, show that, in the absence of fear effect k=0, the equilibria E3 and E∗ are approaching towards stable state by excluding the existence of oscillatory behaviour taking m as the bifurcating parameter. In this manner, predator population also survives in ecosystem for large coefficient of prey refuge parameter m.Figure 14
Impact ofm on the nature of steady states E3 and E∗ when other parameters are fixed as in Figure 9. (a) Nature of steady state E3 (A = 0). (b) Nature of steady state E∗ (A ≠ 0).
(a)(b)Figure 15
Impact ofm on the predator’s growth for different values of k when others are fixed as in Figure 9. ((i) red color: k=0.2, (ii) blue color: k=0.5, and (iii) green color: k=0.8). (a) Nature at E3 (A = 0). (b) Nature at E∗ (A ≠ 0).
(a)(b)Figure 16
Bifurcation diagram aroundE3 with respect to coefficient of refuge parameter m when k=0 and the remaining parameters are fixed as in Figure 4. (a) Bifurcation diagram of x. (b) Bifurcation diagram of y.
(a)(b)Figure 17
Bifurcation diagram aroundE∗ when k=0 and the remaining parameters are fixed as in Figure 9. (a) Bifurcation diagram of x. (b) Bifurcation diagram of w. (c) Bifurcation diagram of y.
(a)(b)(c)A “transcritical bifurcation” (BP) occurs atd1TC=5.5 around E1. At this point, exactly one eigenvalue of the Jacobian matrix is zero and others have negative real parts. Bifurcation diagram 32 depicts that when d1<d1TC, then E1 is unstable, and when d1>d1TC, then E1 is stable. Also, Figure 18 exhibits that when d3<d3TC=0.1768, then E2 is unstable, and when d3>d3TC=0.1768, then E2 is stable. So, a transcritical bifurcation is exhibited at d3TC=0.1768 around E2. Similarly, Figures 19 and 20 depict the transcritical bifurcation diagrams around equilibrium E4 and E5 taking d1d1TC=1.2238 and d3d3TC=0.61 as bifurcation parameter, respectively.Figure 18
Transcritical bifurcation diagram aroundE28.0,0,0 considering d3 as bifurcation parameter and others parameters are the same as in Figure 3. Here d3TC=0.1768.Figure 19
Bifurcation diagram for the transcritical bifurcation aroundE40,1.9889,5.5417 considering d1 as the independent variable and others parameters are the same as in Figure 7. Here, d1TC=1.2238.Figure 20
Transcritical bifurcation diagram aroundE58.3333,6.6667,0 considering d3 as independent variable and others are the same as in Figure 8. Here, d3TC=0.61.


### 6.1. Effect of Time-Delay on Population Dynamics

Now, let us perform the numerical computations to validate the analytical results of the delayed model (60). For the parameter set {r=5.5, d1=0.4, d2=0.3, d3=0.2, a1=0.6, a2=0.98, a3=0.8, b1=2.5, c1=0.85, c2=0.7, A=0, a=1.1, m=0.4}, equation (72) has exactly one positive root 0.2246 (correct up to four decimal places). So, from Theorem 18, the planer equilibrium point E31.20053,0,6.49267 is stable when τ∈0,τ′∗=2.9272 and unstable when τ>τ′∗=2.9272. Stable time series and stable phase trajectory of E31.20053,0,6.49267 are shown in Figure 21 when τ=2<τ′∗=2.9272. Also, Figure 22 depicts the corresponding unstable behaviour of E3 when τ=3.5>τ′∗=2.9272. Moreover, Figure 23 presents the supercritical Hopf bifurcation diagram around E3 taking τ as bifurcation parameter.Figure 21
(a) Stable nature with time and (b) stable spiral ofE31.20053,0,6.49267 when τ=2<τ′∗=2.9272 corresponding to the data set of Figure 4.
(a)(b)Figure 22
(a, b) Oscillatory nature ofx, y, respectively, and (c) stable limit cycle around E31.20053,0,6.49267 when τ=3.5>τ′∗=2.9272 corresponding to the data set of Figure 4. (a) Time series of prey. (b) Time series of predator. (c) Phase diagram.
(a)(b)(c)Figure 23
Supercritical Hopf bifurcation diagram aroundE31.20053,0,6.49267 considering τ as bifurcation parameter and other parameters are chosen from the data set of Figure 4. (a) Bifurcation diagram of prey. (b) Bifurcation diagram of predator.
(a)(b)For the parameter set {r=5.5, d1=0.4, d2=0.3, d3=0.2, a1=0.6, a2=0.98, a3=0.8, b1=2.5, c1=0.85, c2=0.7, A=2, a=1.1, m=0.4}, equation (91) has exactly two positive roots, 0.1127 and 0.0526 (correct up to four decimal places). Then from Theorem 20, we have calculated τ0+=4.0652, τ0−=10.3836, τ1+=22.7845, τ1−=37.7710, and τ2+=41.5037. The interior equilibrium E∗0.385717,0.950363,8.34509 is locally asymptotically stable when τ∈0,4.0652,10.3836,22.7845,37.7710,41.5037 and unstable when τ∈4.0652,10.3836,22.7845,37.7710 and τ>41.5037. At τ=τk±,k=0,1 and τ2+, Hopf bifurcation appears around E∗. Figures 24–26 depict the stable nature of E∗ for τ∈0,4.0652,10.3836,22.7845,37.7710,41.5037 respectively. Also, unstable behaviour of E∗ is presented in Figures 27–29 for τ∈4.0652,10.3836,22.7845,37.7710 and τ>τ2+=41.5037, respectively. The corresponding bifurcation diagrams are depicted in Figures 30 and 31 .Figure 24
(a) Stable behaviour with time and (b) stable trajectory ofE∗0.385717,0.950363,8.34509 when τ=2.5∈0,τ0+=4.0652 corresponding to the data set of Figure 9.
(a)(b)Figure 25
(a) Stable nature with time and (b) stable phase portrait ofE∗0.385717,0.950363,8.34509 when τ=14.5∈τ0−=10.3836,τ1+=22.7845 corresponding to the data set of Figure 9.
(a)(b)Figure 26
(a) Stable behaviour with time and (b) stable trajectory ofE∗0.385717,0.950363,8.34509 when τ=39∈τ1−=37.7710,τ2+=41.5037 corresponding to the data set of Figure 9.
(a)(b)Figure 27
(a, b, c) Oscillatory nature of prey, subsidy, and predator, respectively, with time and (d) stable limit cycle aroundE∗0.385717,0.950363,8.34509 when τ=5∈τ0+=4.0652,τ0−=10.3836 corresponding to the data set of Figure 9. (a) Oscillation of prey. (b). Oscillation of subsidy. (c) Oscillation of predator. (d) Phase portrait (closed orbit).
(a)(b)(c)(d)Figure 28
(a–c) oscillatory nature of prey, subsidy, and predator, respectively, with time and (d) stable limit cycle aroundE∗0.385717,0.950363,8.34509 when τ=25∈τ1+=22.7845,τ1−=37.7710 corresponding to the data set of Figure 9. (a) Oscillation of prey. (b). Oscillation of subsidy. (c) Oscillation of predator. (d) Phase portrait (closed orbit).
(a)(b)(c)(d)Figure 29
(a–c) Oscillatory nature of prey, subsidy, and predator, respectively, with time aroundE∗0.385717,0.950363,8.34509 when τ=50>τ2+=41.5037 corresponding to the data set of Figure 9. (a) Oscillation of prey. (b). Oscillation of subsidy. (c) Oscillation of predator.
(a)(b)(c)Figure 30
Bifurcation diagram aroundE∗ considering τ as independent variable and the remaining parameters are fixed as in dataset of Figure 9. Here, τ runs from 0 to 30. (a) Bifurcation diagram of x. (b) Bifurcation diagram of w. (c) Bifurcation diagram of y.
(a)(b)(c)Figure 31
Bifurcation diagram aroundE∗ regarding τ as bifurcation parameter when τ runs from 30 to 70 and parameters are chosen from the data set of Figure 9. (a) Bifurcation diagram of x. (b) Bifurcation diagram of w. (c) Bifurcation diagram of y.
(a)(b)(c)


## 6.1. Effect of Time-Delay on Population Dynamics

Now, let us perform the numerical computations to validate the analytical results of the delayed model (60). For the parameter set {r=5.5, d1=0.4, d2=0.3, d3=0.2, a1=0.6, a2=0.98, a3=0.8, b1=2.5, c1=0.85, c2=0.7, A=0, a=1.1, m=0.4}, equation (72) has exactly one positive root 0.2246 (correct up to four decimal places). So, from Theorem 18, the planer equilibrium point E31.20053,0,6.49267 is stable when τ∈0,τ′∗=2.9272 and unstable when τ>τ′∗=2.9272. Stable time series and stable phase trajectory of E31.20053,0,6.49267 are shown in Figure 21 when τ=2<τ′∗=2.9272. Also, Figure 22 depicts the corresponding unstable behaviour of E3 when τ=3.5>τ′∗=2.9272. Moreover, Figure 23 presents the supercritical Hopf bifurcation diagram around E3 taking τ as bifurcation parameter.Figure 21
(a) Stable nature with time and (b) stable spiral ofE31.20053,0,6.49267 when τ=2<τ′∗=2.9272 corresponding to the data set of Figure 4.
(a)(b)Figure 22
(a, b) Oscillatory nature ofx, y, respectively, and (c) stable limit cycle around E31.20053,0,6.49267 when τ=3.5>τ′∗=2.9272 corresponding to the data set of Figure 4. (a) Time series of prey. (b) Time series of predator. (c) Phase diagram.
(a)(b)(c)Figure 23
Supercritical Hopf bifurcation diagram aroundE31.20053,0,6.49267 considering τ as bifurcation parameter and other parameters are chosen from the data set of Figure 4. (a) Bifurcation diagram of prey. (b) Bifurcation diagram of predator.
(a)(b)For the parameter set {r=5.5, d1=0.4, d2=0.3, d3=0.2, a1=0.6, a2=0.98, a3=0.8, b1=2.5, c1=0.85, c2=0.7, A=2, a=1.1, m=0.4}, equation (91) has exactly two positive roots, 0.1127 and 0.0526 (correct up to four decimal places). Then from Theorem 20, we have calculated τ0+=4.0652, τ0−=10.3836, τ1+=22.7845, τ1−=37.7710, and τ2+=41.5037. The interior equilibrium E∗0.385717,0.950363,8.34509 is locally asymptotically stable when τ∈0,4.0652,10.3836,22.7845,37.7710,41.5037 and unstable when τ∈4.0652,10.3836,22.7845,37.7710 and τ>41.5037. At τ=τk±,k=0,1 and τ2+, Hopf bifurcation appears around E∗. Figures 24–26 depict the stable nature of E∗ for τ∈0,4.0652,10.3836,22.7845,37.7710,41.5037 respectively. Also, unstable behaviour of E∗ is presented in Figures 27–29 for τ∈4.0652,10.3836,22.7845,37.7710 and τ>τ2+=41.5037, respectively. The corresponding bifurcation diagrams are depicted in Figures 30 and 31 .Figure 24
(a) Stable behaviour with time and (b) stable trajectory ofE∗0.385717,0.950363,8.34509 when τ=2.5∈0,τ0+=4.0652 corresponding to the data set of Figure 9.
(a)(b)Figure 25
(a) Stable nature with time and (b) stable phase portrait ofE∗0.385717,0.950363,8.34509 when τ=14.5∈τ0−=10.3836,τ1+=22.7845 corresponding to the data set of Figure 9.
(a)(b)Figure 26
(a) Stable behaviour with time and (b) stable trajectory ofE∗0.385717,0.950363,8.34509 when τ=39∈τ1−=37.7710,τ2+=41.5037 corresponding to the data set of Figure 9.
(a)(b)Figure 27
(a, b, c) Oscillatory nature of prey, subsidy, and predator, respectively, with time and (d) stable limit cycle aroundE∗0.385717,0.950363,8.34509 when τ=5∈τ0+=4.0652,τ0−=10.3836 corresponding to the data set of Figure 9. (a) Oscillation of prey. (b). Oscillation of subsidy. (c) Oscillation of predator. (d) Phase portrait (closed orbit).
(a)(b)(c)(d)Figure 28
(a–c) oscillatory nature of prey, subsidy, and predator, respectively, with time and (d) stable limit cycle aroundE∗0.385717,0.950363,8.34509 when τ=25∈τ1+=22.7845,τ1−=37.7710 corresponding to the data set of Figure 9. (a) Oscillation of prey. (b). Oscillation of subsidy. (c) Oscillation of predator. (d) Phase portrait (closed orbit).
(a)(b)(c)(d)Figure 29
(a–c) Oscillatory nature of prey, subsidy, and predator, respectively, with time aroundE∗0.385717,0.950363,8.34509 when τ=50>τ2+=41.5037 corresponding to the data set of Figure 9. (a) Oscillation of prey. (b). Oscillation of subsidy. (c) Oscillation of predator.
(a)(b)(c)Figure 30
Bifurcation diagram aroundE∗ considering τ as independent variable and the remaining parameters are fixed as in dataset of Figure 9. Here, τ runs from 0 to 30. (a) Bifurcation diagram of x. (b) Bifurcation diagram of w. (c) Bifurcation diagram of y.
(a)(b)(c)Figure 31
Bifurcation diagram aroundE∗ regarding τ as bifurcation parameter when τ runs from 30 to 70 and parameters are chosen from the data set of Figure 9. (a) Bifurcation diagram of x. (b) Bifurcation diagram of w. (c) Bifurcation diagram of y.
(a)(b)(c)


## 7. Conclusion

We have analyzed a system for generalist predator which utilizes more than one food source: predator-prey-subsidy model of non-Kolmogorov form introducing nonlinear prey refuge function and the effect of fear felt by prey population. Our main interest is to find the situations such that dynamical stability and instability appear so as to make out more fully how subsidy may influence the predator and their prey. It has been shown that the solutions of system (2) remain positively invariant always and they are asymptotically uniformly bounded. These, in turn, imply that system (2) is biologically well-behaved. Existence criteria and stable behaviour of all the biologically meaningful equilibria have been discussed. It has to be noted that Hopf bifurcations are exhibited around E3 (subsidy free) and E∗ (interior) of system (2) considering k as a bifurcating parameter (see Figures 4–6, 9, 11, and 12). Also, observing Figures 6 and 12, it can be concluded that high levels of fear can stabilize system (2) by excluding the existence of periodic solutions. These phenomena are biologically significant because prey species are aware after a certain level of fear; i.e., after a certain level of fear, they are not affected as they are aware and show signs of habituation.Moreover, this work derives transcritical bifurcations (local bifurcation of codimension 1) at the various equilibrium pointsE1, E2, E4, and E5, respectively (see Figures 18–20 and 32).Figure 32
Transcritical bifurcation diagram aroundE10,0.5,0 considering d1 as bifurcation parameter and others parameters are the same as in Figure 2. Here d1TC=5.5.Also, we have discussed numerically the influences of coefficient of prey refuge parameterm on the nature of the equilibrium points E3 (zero subsidy input rate) and E∗ (fixed small subsidy input rate) irrespective of fear level k. Noting Figures 14–17, it is observed that both the prey and predator species always persist in ecosystem due to continuous increment of coefficient of prey refuge. But Figures 10 and 13 depict that, irrespective of fear level, a highly subsidized predator should indeed drive the prey population towards extinction regardless of whether the prey and subsidy arise in the same habitat. This phenomenon is ecologically meaningful because the prey population cannot get enough time to protect themselves from predation risk for enormous growth of predator at high subsidy input rate. So, the prey population is leading towards extinction, but the predator species can easily survive in ecosystem with the help of resource subsidy. Thus, the study of system (2) is ecologically very significant.In reality, fear effect does not instantaneously reduce the birth of a prey population, but some time lag should be needed to create an impact on the birth rate of the prey population. We have considered that there is a time-delay on the impact of fear to the birth rate of prey, from the instance it perceives the fear of predator through any means. So, the incorporation of time-delay makes system (60) more realistic. It is noted that delay parameter τ has a significant role because there exists a threshold value τ′∗ such that stable behaviour of planer equilibrium point E3 (in the absence of subsidy input rate) switches to oscillatory nature when τ passes through its threshold value τ′∗; i.e., the vector fields for τ<τ′∗ and τ>τ′∗ are qualitatively different. So, system (60) exhibits a supercritical Hopf bifurcation around E3 considering τ as bifurcation parameter (see Figures 21–23). Also, a rigorous study of the stability and bifurcation of interior equilibrium point E∗ has been performed. Our analysis describes that the delay within a certain specified range could maintain the stable behaviour of E∗. On the other hand, the delay could drive the system into an unstable state. Hence, the study of the time-delay parameter has a regulatory impact on the whole system.



---

*Source: 9963031-2021-07-05.xml*