# Dynamics on Effect of Prey Refuge Proportional to Predator in Discrete-Time Prey-Predator Model

Authors: G. S. Mahapatra; P. K. Santra; Ebenezer Bonyah
Journal: Complexity
 (2021)

---

## Abstract

Prey-predator models with refuge effect have great importance in the context of ecology. Constant refuge and refuge proportional to prey are the most popular concepts of refuge in the existing literature. Now, there are new different types of refuge concepts attracting researchers. This study considers a refuge concept proportional to the predator due to the fear induced by predators. When predators increase, fears also increase and that is why prey refuges also increase. Here, we examine the influence of prey refuge proportional to predator effect in a discrete prey-predator interaction with the Holling type II functional response model. Is this refuge stabilizing or destabilizing the system? That is the central question of this study. The existence and stability of fixed points, Period-Doubling Bifurcation, Neimark–Sacker Bifurcation, the influence of prey refuge, and chaos are analyzed. This work provides the bifurcation diagrams and Lyapunov exponents to analyze the refuge parameter of the model. The proposed discrete model indicates rich dynamics as the effect of prey refuge through numerical simulations.


---

## Body

## 1. Introduction

In the current biomathematical literature, prey-predator interaction has become an exciting subject matter due to its influence on the environment. Many researchers [1–4] have devoted considerable time to explore several perspectives of the dynamical behaviour of this subject matter in ecology and the associated growth of population models [5–7]. Nature offers some amount of protection to some prey populations in the form of refuge dimensions. In some instances, refugia give a prolonged prey-predator interaction [8, 9] to reduce the probability of extinction due to predation activity. This phenomenon has been investigated extensively by many researchers in the discrete domain of refuse concepts [10–12].In the context of populations characterized by overlapping generations, the nature of the birth process is a continuous matter; therefore, the predator-prey relationship is frequently developed through different models utilizing a deterministic approach such as ordinary differential equations. Several species, including monocarpic plants and semelparous animals, possess a discrete nonoverlapping generation character and have birth and regular breeding seasons. Their interactions are presented in the form of difference equations or other forms such as discrete-time mappings [13–17]. A discrete-time prey-predator model is usually characterized by more complicated dynamics behaviour than the associated continuous-time models [18–21]. Several researchers have focused on this concept to present a comprehensive rich dynamics of this phenomenon, including stability analysis of equilibria [22–27], Period-Doubling Bifurcation [28], Neimark–Sacker bifurcation [29], and chaos control [30]. Liu and Cai [31] explored the presence of bifurcation and chaos in a discrete-time prey-predator system. Elettreby et al. [32] used mixed functional response to investigate the dynamics of a discrete prey-predator model. AlSharawi et al. [33] studied a discrete prey-predator model with a nonmonotonic functional response and strong Allee effect in prey. Wang and Fečkan [34] studied the dynamics of a discrete nonlinear prey-predator model. Rech [35] investigated two discrete-time counterparts of a continuous-time prey-predator model. Khan and Khalique [36] discussed Neimark–Sacker bifurcation and hybrid control in a discrete-time Lotka–Volterra model. For some more dynamical investigations related to different versions of prey-predator models, we refer to Elsadany et al. [37]; Baydemir et al. [38]; Rozikov and Shoyimardonov [39]; Singh and Deolia [40]; and references therein.In this paper, we consider a refuge term proportional to predator density due to fear induced by predators. The earlier literature has pointed out that the use of refuge by a fixed number or a small percentage of prey exerts a stabilizing effect on the dynamics of the interacting populations. This article is concerned with the above statement assuming that the quantity of prey in refugia is proportional to that of predators due to the fear induced by the predator. We analyze the dynamic properties of such a prey-predator model with prey self-limitation. The rest of the paper is partitioned as follows. A discrete-time prey-predator model with refuge is formulated in the second section. Section3 is devoted to the local stability analysis of the model. Bifurcation analysis of the proposed model is presented in Section 4. Section 5 deals with the influence of prey refuge on the proposed model. Section 6 gives a chaos control procedure. Section 7 presents numerical simulations with support of the presented dynamical analysis of the proposed model. Section 8 presents the conclusion that is the last aspect of the findings from this study.


## 2. Development of a Discrete Prey-Predator System

The densities of prey and predator populations vary with time, have a uniform distribution over space, and have no stage structure for neither prey nor predators. The proposed model takes into consideration a generalised prey-predator model incorporating the logistic growth of prey with functional responseφx of the predator population as dx/dt=rx1−x/k−cφxy, and the predator as a population is epitomised by the expression dy/dt=dφxy−fy, where x and y are the density of prey and predator populations correspondingly at any time t. The parameters r,k,c,d,andf are all taken to be positive constants and have its appropriate biological interpretations, respectively. r denotes the intrinsic per capita growth rate of prey population, k represents the environmental carrying capacity of prey population, c denotes the maximal per capita consumption rate of predators, d considers the efficiency with which predators transform consumed prey into new predators, and f represents the per capita death rate of predators. φx deals with the functional response of the predator population and fulfils the assumption φ0=0andφ′x>0 for x>0. There exists a quantity xr of prey population that incorporates refuges, and therefore we have modified the functional response to incorporate prey refuges, and we consider φx=x/e+x, and xr=by. The mathematical form of the proposed predator-prey model is as follows:(1)dxdt=rx1−xk−cx−byye+x−by,dydt=dx−byye+x−by−fy.To arrive at a discrete-time model from (1), let dx/dt=xt+h−xt/h and dy/dt=yt+h−yt/h, where xt and yt are the densities of the prey and predator populations, respectively, in discrete time t. Let h⟶1 and f=1, we obtain the equations for the n+1th generation of the prey and predator populations substituting t by n as follows:(2)xn+1=r+1xn1−rkr+1xn−cxn−bynyne+xn−byn,yn+1=dxn−bynyne+xn−byn.Letr/kr+1=1 and r+1=a, we obtain the discrete-time predator-prey system from (2) as follows:(3)xn+1=axn1−xn−cxn−bynyne+xn−byn,yn+1=dxn−bynyne+xn−byn,where a,b,c,d, and e are all positive constants. By the biological meaning of the model variables, we only consider the system in the region Ω=x,y:x≥0,y≥0 in the x,y−plane.


## 3. Fixed Points and Stability Analysis

Fixed points of the discrete system (3) are determined by solving the following nonlinear system of equations:(4)x=ax1−x−cx−byye+x−by,y=dx−byye+x−by.We get three nonnegative fixed points by solving above equations: (i)P0=0,0, (ii) P1=a−1/a,0, a>1, and (iii) P2=x2,y2; here, x2 is positive root of the equation x2+Ax+B=0, where A=1−a/a+c/abd, B=−ce/bdd−1, and y2=x2/b−e/bd−1.


### 3.1. Dynamic Behaviour of the Discrete System

In this part, we examine the local behaviour of the discrete system (3) for every steady point of the model. The stability of system (3) is accomplished by deriving a Jacobian matrix connected to each equilibrium point. The Jacobian Matrix J for system (3) is(5)J=a11a12a21a22,where a11=a1−2x−cey/e+x−by2, a12=−cex−2by+x−by2/e+x−by2, a21=dey/e+x−by2, and a22=dex−2by+x−by2/e+x−by2.The characteristic equation of the matrixJ is λ2−TrJλ+DetJ=0, where TrJ = trace of matrix = a1−2x+dex−2by+x−by2−cey/e+x−by2 and DetJ = determinant of matrix = a1−2xdex−2by+x−by2/e+x−by2.Therefore, the discrete-time system (3) possesses the following:(i)
A dissipative dynamical system ifa1−2xdex−2by+x−by2/e+x−by2<1(ii)
A conservative dynamical system if and only ifa1−2xdex−2by+x−by2/e+x−by2=1(iii)
An undissipated dynamical system otherwise


#### 3.1.1. Stability and Dynamic Behaviour of the Model atP0

The dynamical behaviour of the system is studied utilizing the Jacobian matrix at the fixed pointP00,0. The Jacobian matrix at P00,0 is J=a000. The equilibrium point P00,0 is called (i) sink if a<1; (ii) saddle if a>1; and (iii) nonhyperbolic if a=1.


#### 3.1.2. Stability and Dynamic Behaviour atP1

The Jacobian matrix at the fixed pointP1a−1/a,0 is J=2−a−ca−1/ae+a−10da−1/ae+a−1. The equilibrium point P1 is said to be(i)
Sink if2−a<1 and da−1/ae+a−1<1(ii)
Source if2−a>1 and da−1/ae+a−1>1(iii)
Saddle if2−a>1 and da−1/ae+a−1<1 or 2−a<1 and da−1/ae+a−1>1(iv)
Nonhyperbolic if2−a=1 or da−1/ae+a−1=1


#### 3.1.3. Dynamic Behaviour around the Interior Fixed Point

Following the Jacobian matrix at the interior fixed pointP2x2,y2, we get(6)1−TrL+DetL=1−a1−2x2−dex2−2by2+x2−by22−cey2e+x2−by22+a1−2x2dex2−2by2+x2−by22e+x2−by22,1+TrL+DetL=1+a1−2x2+dex2−2by2+x2−by22−cey2e+x2−by22+a1−2x2dex2−2by2+x2−by22e+x2−by22.If1−TrL+DetL>0, then interior equilibrium point P2x2,y2 is as follows:(i)
Sink if1+TrL+DetL>0 and DetL<1(ii)
Source if1+TrL+DetL>0 and DetL>1(iii)
Saddle if1+TrL+DetL<0(iv)
Nonhyperbolic if1+TrL+DetL=0 and TrL≠0,2 or TrL2−4DetL<0 and DetL=1If1−TrL+DetL>0, 1+TrL+DetL=0, and TrL≠0,2, then P2x2,y2 can undergo Period-Doubling Bifurcation (PDB).If1−TrL+DetL>0, TrL2−4DetL<0, and DetL=1, then P2x2,y2 can undergo Neimark–Sacker Bifurcation (NSB).


## 3.1. Dynamic Behaviour of the Discrete System

In this part, we examine the local behaviour of the discrete system (3) for every steady point of the model. The stability of system (3) is accomplished by deriving a Jacobian matrix connected to each equilibrium point. The Jacobian Matrix J for system (3) is(5)J=a11a12a21a22,where a11=a1−2x−cey/e+x−by2, a12=−cex−2by+x−by2/e+x−by2, a21=dey/e+x−by2, and a22=dex−2by+x−by2/e+x−by2.The characteristic equation of the matrixJ is λ2−TrJλ+DetJ=0, where TrJ = trace of matrix = a1−2x+dex−2by+x−by2−cey/e+x−by2 and DetJ = determinant of matrix = a1−2xdex−2by+x−by2/e+x−by2.Therefore, the discrete-time system (3) possesses the following:(i)
A dissipative dynamical system ifa1−2xdex−2by+x−by2/e+x−by2<1(ii)
A conservative dynamical system if and only ifa1−2xdex−2by+x−by2/e+x−by2=1(iii)
An undissipated dynamical system otherwise


### 3.1.1. Stability and Dynamic Behaviour of the Model atP0

The dynamical behaviour of the system is studied utilizing the Jacobian matrix at the fixed pointP00,0. The Jacobian matrix at P00,0 is J=a000. The equilibrium point P00,0 is called (i) sink if a<1; (ii) saddle if a>1; and (iii) nonhyperbolic if a=1.


### 3.1.2. Stability and Dynamic Behaviour atP1

The Jacobian matrix at the fixed pointP1a−1/a,0 is J=2−a−ca−1/ae+a−10da−1/ae+a−1. The equilibrium point P1 is said to be(i)
Sink if2−a<1 and da−1/ae+a−1<1(ii)
Source if2−a>1 and da−1/ae+a−1>1(iii)
Saddle if2−a>1 and da−1/ae+a−1<1 or 2−a<1 and da−1/ae+a−1>1(iv)
Nonhyperbolic if2−a=1 or da−1/ae+a−1=1


### 3.1.3. Dynamic Behaviour around the Interior Fixed Point

Following the Jacobian matrix at the interior fixed pointP2x2,y2, we get(6)1−TrL+DetL=1−a1−2x2−dex2−2by2+x2−by22−cey2e+x2−by22+a1−2x2dex2−2by2+x2−by22e+x2−by22,1+TrL+DetL=1+a1−2x2+dex2−2by2+x2−by22−cey2e+x2−by22+a1−2x2dex2−2by2+x2−by22e+x2−by22.If1−TrL+DetL>0, then interior equilibrium point P2x2,y2 is as follows:(i)
Sink if1+TrL+DetL>0 and DetL<1(ii)
Source if1+TrL+DetL>0 and DetL>1(iii)
Saddle if1+TrL+DetL<0(iv)
Nonhyperbolic if1+TrL+DetL=0 and TrL≠0,2 or TrL2−4DetL<0 and DetL=1If1−TrL+DetL>0, 1+TrL+DetL=0, and TrL≠0,2, then P2x2,y2 can undergo Period-Doubling Bifurcation (PDB).If1−TrL+DetL>0, TrL2−4DetL<0, and DetL=1, then P2x2,y2 can undergo Neimark–Sacker Bifurcation (NSB).


## 3.1.1. Stability and Dynamic Behaviour of the Model atP0

The dynamical behaviour of the system is studied utilizing the Jacobian matrix at the fixed pointP00,0. The Jacobian matrix at P00,0 is J=a000. The equilibrium point P00,0 is called (i) sink if a<1; (ii) saddle if a>1; and (iii) nonhyperbolic if a=1.


## 3.1.2. Stability and Dynamic Behaviour atP1

The Jacobian matrix at the fixed pointP1a−1/a,0 is J=2−a−ca−1/ae+a−10da−1/ae+a−1. The equilibrium point P1 is said to be(i)
Sink if2−a<1 and da−1/ae+a−1<1(ii)
Source if2−a>1 and da−1/ae+a−1>1(iii)
Saddle if2−a>1 and da−1/ae+a−1<1 or 2−a<1 and da−1/ae+a−1>1(iv)
Nonhyperbolic if2−a=1 or da−1/ae+a−1=1


## 3.1.3. Dynamic Behaviour around the Interior Fixed Point

Following the Jacobian matrix at the interior fixed pointP2x2,y2, we get(6)1−TrL+DetL=1−a1−2x2−dex2−2by2+x2−by22−cey2e+x2−by22+a1−2x2dex2−2by2+x2−by22e+x2−by22,1+TrL+DetL=1+a1−2x2+dex2−2by2+x2−by22−cey2e+x2−by22+a1−2x2dex2−2by2+x2−by22e+x2−by22.If1−TrL+DetL>0, then interior equilibrium point P2x2,y2 is as follows:(i)
Sink if1+TrL+DetL>0 and DetL<1(ii)
Source if1+TrL+DetL>0 and DetL>1(iii)
Saddle if1+TrL+DetL<0(iv)
Nonhyperbolic if1+TrL+DetL=0 and TrL≠0,2 or TrL2−4DetL<0 and DetL=1If1−TrL+DetL>0, 1+TrL+DetL=0, and TrL≠0,2, then P2x2,y2 can undergo Period-Doubling Bifurcation (PDB).If1−TrL+DetL>0, TrL2−4DetL<0, and DetL=1, then P2x2,y2 can undergo Neimark–Sacker Bifurcation (NSB).


## 4. Bifurcation Analysis

In this section, we investigate the existence of NSB and PDB at the positive fixed pointP2x2,y2 of the discrete-time of (3). Countless numbers of dynamical properties of the system can be analyzed with regard to the emergence of NSB and PDB. Bifurcation frequently occurs when there is a change in the stability of a fixed point. Thus, the qualitative properties of a dynamical system vary. The NSB and PDB are explored in the positive fixed point direction P2x2,y2 of system (3) considering refuge parameter b as a bifurcation parameter. Utilizing the emergence of NSB, closed invariant circles are obtained. The bifurcation can be classified as supercritical or subcritical, leading to a stable or unstable closed invariant curve, in that order. A PDB in a discrete dynamical system is a bifurcation where a small perturbation in a parameter value brings about the system moving to a new behaviour with the double period of the original system.


### 4.1. Neimark–Sacker Bifurcation

Here, we examine the NSB of the discrete prey-predator model (3) at P2x2,y2 for the parameters lying in the following set ANSB = {a,b=bNSB,c,d,e:1−TrL+DetL>0, TrL2−4DetL<0, DetL=1}.In exploring the NSB,b is utilised as the bifurcation parameter. Further, b∗ is the perturbation of b where b∗⋘1, and perturbation of the model is considered as follows:(7)xn+1=axn1−xn−cxn−b+b∗ynyne+xn−b+b∗yn≡fxn,yn,b∗,yn+1=dxn−b+b∗ynyne+xn−b+b∗yn≡gxn,yn,b∗.Letun=xn−x2andvn=yn−y2, then equilibrium P2x2,y2 is converted into the origin, and further making expansion f and g as a Taylor series at un,vn=0,0 to the third order, model (7) turns to(8)un+1=α1un+α2vn+α11un2+α12unvn+α22vn2+α111un3+α112un2vn+α122unvn2+α222vn3+Oun+vn4,vn+1=β1un+β2vn+β11un2+β12unvn+β22vn2+β111un3+β112un2vn+β122unvn2+β222vn3+Oun+vn4,where α1=fxx2,y2,0, α2=fyx2,y2,0, α11=fxxx2,y2,0, α12=fxyx2,y2,0, α22=fyyx2,y2,0, α111=fxxxx2,y2,0, α112=fxxyx2,y2,0, α122=fxyyx2,y2,0, α222=fyyyx2,y2,0, β1=gxx2,y2,0, β2=gyx2,y2,0, β11=gxxx2,y2,0, β12=gxyx2,y2,0, β22=gyyx2,y2,0, β111=gxxxx2,y2,0, β112=gxxyx2,y2,0, β122=gxyyx2,y2,0, and β222=gyyyx2,y2,0.Note that the characteristic equation associated with the linearization of model (8) at un,vn=0,0 is given by λ2−TrL1b∗λ+DetL1b∗=0. The roots of the characteristic equation are λ1,2b∗=TrL1b∗±i4DetL1b∗−TrL1b∗2/2.Fromλ1,2b∗=1 and b∗=0, we have λ1,2b∗=DetL1b∗1/2 and l=dλ1,2b∗/db∗h∗=0≠0.In addition, it is required that whenb∗=0, λ1,2i≠1, i=1,2,3,4, which is equivalent to TrL10≠−2,−1,1,2.To study the normal form, letγ=Imλ1,2 and δ=Reλ1,2. We define T=01γδ, and using the transformation unvn=Tx¯ny¯n, model (8) becomes(9)x¯n+1=δx¯n−γy¯n+f1x¯n,y¯n,y¯n+1=γx¯n+δy¯n+g1x¯n,y¯n,where the functions f1 and g1 denote the terms in model (9) in variables x¯n,y¯n with the order at least two.To undergo NSB, we require that the following discriminatory quantityΩ be nonzero:(10)Ω=−Re1−2λ¯λ¯21−λξ11ξ20−12ξ112−ξ022+Reλ¯ξ21,where(11)ξ20=δ82β22−δα22−α12+4γα22+i4γα22−2α22−2δα22+γ4α12+i184γβ22+2γ2α22−2α11+β128+δα11−2β11+δ3α22−δ2β22−δ2α12+δβ124γ,ξ11=γ2β22−δα22+i12γ2α22+α11+δα12+δ2α22+β11−δα11+δβ12−δ2α12−2δ2β22+2δ3α222γ,ξ02=14γ2δα22+α12+β22+i14β12+2δβ22−2δα12−α11−β11−δα11+δβ12−δ2α12+δ2β22−δ3α224γ+14α22iγ2−3δ2,ξ21=38β222γ2+δ2+β1128+δ4α112+δ4β122+α122γ28+3δ28−δ4+38α111+i38α222γ2+2δ2+i3γδ8α122−18β122γi−i3γδ8β222−i3β111−3δα1118γ−i3δβ112−3δ2α1128γ−i3δ2β122−3δ3α1228γ−i3δ3β222−3δ4α2228γ.Finally, the above analysis leads to the following associated results.Theorem 1.
IfΩ≠0, then model (3) undergoes NSB at P2x2,y2 when the parameter b changes in a small neighbourhood of bNSB. If Ω<0Ω>0, then an attracting (repelling) invariant closed curve bifurcates from P2x2,y2.


### 4.2. Period-Doubling Bifurcation

Here, we consider one of the eigenvalues of the positive fixed point thatP2x2,y2 is λ1=−1, and the other λ2 is neither 1 nor −1, whenever parameters of the model are located in the following set APDB = {a,b=bPDB,c,d,e: 1−TrJ+DetJ>0, 1+TrJ+DetJ=0,TrJ≠0,2}. Here, we discuss PDB of model (3) at P2x2,y2 when parameters change in a small neighbourhood of APDB. In studying the flip bifurcation, b is made use as the bifurcation parameter. Furthermore, b∗ is the perturbation of b where b∗⋘1, and perturbation of the model is considered as follows:(12)xn+1=axn1−xn−cxn−b+b∗ynyne+xn−b+b∗yn≡fxn,yn,b∗,yn+1=dxn−b+b∗ynyne+xn−b+b∗yn≡gxn,yn,b∗.Letun=xn−x2 and vn=yn−y2; then, the equilibrium P2x2,y2 is transformed into the origin, and further expanding f and g as a Taylor series at un,vn,b∗=0,0,0 up to the third order, system (12) becomes(13)un+1=α1un+α2vn+α11un2+α12unvn+α13unb∗+α23vnb∗+α111un3+α112un2vn+α113un2b∗+α123unvnb∗+Oun+vn+b∗4,vn+1=β1un+β2vn+β11un2+β12unvn+β22vn2+β13unb∗+β23vnb∗+β111un3+β112un2vn+β113un2b∗+β123unvnb∗+β223vn2b∗+Oun+vn+b∗4,where α1=fxx2,y2,0, α2=fyx2,y2,0, α11=fxxx2,y2,0, α12=fxyx2,y2,0, α13=fxb∗x2,y2,0, α23=fyb∗x2,y2,0, α111=fxxxx2,y2,0, α112=fxxyx2,y2,0, α113=fxxb∗x2,y2,0, α123=fxyb∗x2,y2,0, β1=gxx2,y2,0, β2=gyx2,y2,0, β11=gxxx2,y2,0, β12=gxyx2,y2,0, β22=gyyx2,y2,0, β13=gxb∗x2,y2,0, β23=gyb∗x2,y2,0, β111=gxxxx2,y2,0, β112=gxxyx2,y2,0, β113=gxxb∗x2,y2,0, β123=gxyb∗x2,y2,0, and β223=gyyb∗x2,y2,0.We defineT=α2α2−1−α1λ2−α1, and it is obvious that T is invertible. Using the transformation unvn=Tx¯ny¯n, system (13) becomes(14)x¯n+1=−x¯n+f1un,vn,b∗,y¯n+1=λ2y¯n+g1un,vn,b∗,where the functions f1 and g1 denote the terms in model (14) in variables un,vn,b∗ with order at least two.Using the centre manifold theorem, we know that there exists a centre manifoldWc0,0,0 of system (14) at 0,0 in a small neighbourhood of b∗=0, which can be approximately described as follows:(15)Wc0,0,0=x¯n,y¯n,b∗εR3:y¯n+1=α¯1x¯n2+α¯2x¯nb∗+Ox¯n+b∗3,where(16)α¯1=α21+α1α11+α2β111−λ22+β221+α121−λ22−1+α1α121+α1+α2β121−λ22,α¯2=1+α1α231+α1+α2β23α21+λ22−1+α1α13+α2β131+λ22.System (14) which is restricted to the centre manifold Wc0,0,0 has the following form:(17)x¯n+1=−x¯n+h1x¯n2+h2x¯nb∗+h3x¯n2b∗+h4x¯nb∗2+h5x¯n3+Ox¯n+b∗3≡Fx¯n,b∗,where(18)h1=α¯2λ2−α¯1α11−α¯2β111+λ2−β221+α¯121+λ2−1+α¯1λ2−α¯1α12−α¯2β121+λ2,h2=λ2−α¯1α13−α¯2β131+λ2−1+α¯1λ2−α¯1α23−α¯2β23α¯21+λ2,h3=λ2−α1α¯1α13−α2β131+λ2+λ2−α1α23−α2β23λ2−α1α¯1α21+λ2−1+α1λ2−α1α123−α2β1231+λ2+α2λ2−α1α113−α2β1131+λ2−β2231+α121+λ2+2α2α¯2λ2−α1α11−α2β111+λ2−2β22α¯21+α1λ2−α11+λ2+α¯2λ2−α1α12−α2β12λ2−1−2α11+λ2,h4=α¯2λ2−α1α13−α2β131+λ2+λ2−α1α23−α2β23λ2−α1α¯2α21+λ2+2α2α¯2λ2−α1α11−α2β111+λ2+2β22α¯21+α1λ2−α11+λ2+α¯2λ2−α1α12−α2β12λ2−1−2α11+λ2,h5=2α2α¯1λ2−α1α11−α2β111+λ2+2β22α¯1λ2−α11+α11+λ2+λ2−α1α11−α2β11λ2−1−2α1α¯11+λ2+α¯22λ2−α1α111−α2β1111+λ2−α¯21+α1λ2−α1α112−α2β1121+λ2.For PDB, the two discriminatory quantitiesξ1 and ξ2 be nonzero,(19)ξ1=∂2F∂x¯∂b∗+12∂F∂b∗∂2F∂x¯2|0,0,ξ2=16∂3F∂x¯3+12∂2F∂x¯22|0,0.Finally, from the above analysis, we have the following result.Theorem 2.
Ifξ1≠0 and ξ2≠0, then model (3) undergoes PDB at P2x2,y2 when the parameter b varies in a small neighbourhood of bPDB. If ξ2>0ξ2<0, then the period-two orbits that bifurcate from P2x2,y2 are stable (unstable).


## 4.1. Neimark–Sacker Bifurcation

Here, we examine the NSB of the discrete prey-predator model (3) at P2x2,y2 for the parameters lying in the following set ANSB = {a,b=bNSB,c,d,e:1−TrL+DetL>0, TrL2−4DetL<0, DetL=1}.In exploring the NSB,b is utilised as the bifurcation parameter. Further, b∗ is the perturbation of b where b∗⋘1, and perturbation of the model is considered as follows:(7)xn+1=axn1−xn−cxn−b+b∗ynyne+xn−b+b∗yn≡fxn,yn,b∗,yn+1=dxn−b+b∗ynyne+xn−b+b∗yn≡gxn,yn,b∗.Letun=xn−x2andvn=yn−y2, then equilibrium P2x2,y2 is converted into the origin, and further making expansion f and g as a Taylor series at un,vn=0,0 to the third order, model (7) turns to(8)un+1=α1un+α2vn+α11un2+α12unvn+α22vn2+α111un3+α112un2vn+α122unvn2+α222vn3+Oun+vn4,vn+1=β1un+β2vn+β11un2+β12unvn+β22vn2+β111un3+β112un2vn+β122unvn2+β222vn3+Oun+vn4,where α1=fxx2,y2,0, α2=fyx2,y2,0, α11=fxxx2,y2,0, α12=fxyx2,y2,0, α22=fyyx2,y2,0, α111=fxxxx2,y2,0, α112=fxxyx2,y2,0, α122=fxyyx2,y2,0, α222=fyyyx2,y2,0, β1=gxx2,y2,0, β2=gyx2,y2,0, β11=gxxx2,y2,0, β12=gxyx2,y2,0, β22=gyyx2,y2,0, β111=gxxxx2,y2,0, β112=gxxyx2,y2,0, β122=gxyyx2,y2,0, and β222=gyyyx2,y2,0.Note that the characteristic equation associated with the linearization of model (8) at un,vn=0,0 is given by λ2−TrL1b∗λ+DetL1b∗=0. The roots of the characteristic equation are λ1,2b∗=TrL1b∗±i4DetL1b∗−TrL1b∗2/2.Fromλ1,2b∗=1 and b∗=0, we have λ1,2b∗=DetL1b∗1/2 and l=dλ1,2b∗/db∗h∗=0≠0.In addition, it is required that whenb∗=0, λ1,2i≠1, i=1,2,3,4, which is equivalent to TrL10≠−2,−1,1,2.To study the normal form, letγ=Imλ1,2 and δ=Reλ1,2. We define T=01γδ, and using the transformation unvn=Tx¯ny¯n, model (8) becomes(9)x¯n+1=δx¯n−γy¯n+f1x¯n,y¯n,y¯n+1=γx¯n+δy¯n+g1x¯n,y¯n,where the functions f1 and g1 denote the terms in model (9) in variables x¯n,y¯n with the order at least two.To undergo NSB, we require that the following discriminatory quantityΩ be nonzero:(10)Ω=−Re1−2λ¯λ¯21−λξ11ξ20−12ξ112−ξ022+Reλ¯ξ21,where(11)ξ20=δ82β22−δα22−α12+4γα22+i4γα22−2α22−2δα22+γ4α12+i184γβ22+2γ2α22−2α11+β128+δα11−2β11+δ3α22−δ2β22−δ2α12+δβ124γ,ξ11=γ2β22−δα22+i12γ2α22+α11+δα12+δ2α22+β11−δα11+δβ12−δ2α12−2δ2β22+2δ3α222γ,ξ02=14γ2δα22+α12+β22+i14β12+2δβ22−2δα12−α11−β11−δα11+δβ12−δ2α12+δ2β22−δ3α224γ+14α22iγ2−3δ2,ξ21=38β222γ2+δ2+β1128+δ4α112+δ4β122+α122γ28+3δ28−δ4+38α111+i38α222γ2+2δ2+i3γδ8α122−18β122γi−i3γδ8β222−i3β111−3δα1118γ−i3δβ112−3δ2α1128γ−i3δ2β122−3δ3α1228γ−i3δ3β222−3δ4α2228γ.Finally, the above analysis leads to the following associated results.Theorem 1.
IfΩ≠0, then model (3) undergoes NSB at P2x2,y2 when the parameter b changes in a small neighbourhood of bNSB. If Ω<0Ω>0, then an attracting (repelling) invariant closed curve bifurcates from P2x2,y2.


## 4.2. Period-Doubling Bifurcation

Here, we consider one of the eigenvalues of the positive fixed point thatP2x2,y2 is λ1=−1, and the other λ2 is neither 1 nor −1, whenever parameters of the model are located in the following set APDB = {a,b=bPDB,c,d,e: 1−TrJ+DetJ>0, 1+TrJ+DetJ=0,TrJ≠0,2}. Here, we discuss PDB of model (3) at P2x2,y2 when parameters change in a small neighbourhood of APDB. In studying the flip bifurcation, b is made use as the bifurcation parameter. Furthermore, b∗ is the perturbation of b where b∗⋘1, and perturbation of the model is considered as follows:(12)xn+1=axn1−xn−cxn−b+b∗ynyne+xn−b+b∗yn≡fxn,yn,b∗,yn+1=dxn−b+b∗ynyne+xn−b+b∗yn≡gxn,yn,b∗.Letun=xn−x2 and vn=yn−y2; then, the equilibrium P2x2,y2 is transformed into the origin, and further expanding f and g as a Taylor series at un,vn,b∗=0,0,0 up to the third order, system (12) becomes(13)un+1=α1un+α2vn+α11un2+α12unvn+α13unb∗+α23vnb∗+α111un3+α112un2vn+α113un2b∗+α123unvnb∗+Oun+vn+b∗4,vn+1=β1un+β2vn+β11un2+β12unvn+β22vn2+β13unb∗+β23vnb∗+β111un3+β112un2vn+β113un2b∗+β123unvnb∗+β223vn2b∗+Oun+vn+b∗4,where α1=fxx2,y2,0, α2=fyx2,y2,0, α11=fxxx2,y2,0, α12=fxyx2,y2,0, α13=fxb∗x2,y2,0, α23=fyb∗x2,y2,0, α111=fxxxx2,y2,0, α112=fxxyx2,y2,0, α113=fxxb∗x2,y2,0, α123=fxyb∗x2,y2,0, β1=gxx2,y2,0, β2=gyx2,y2,0, β11=gxxx2,y2,0, β12=gxyx2,y2,0, β22=gyyx2,y2,0, β13=gxb∗x2,y2,0, β23=gyb∗x2,y2,0, β111=gxxxx2,y2,0, β112=gxxyx2,y2,0, β113=gxxb∗x2,y2,0, β123=gxyb∗x2,y2,0, and β223=gyyb∗x2,y2,0.We defineT=α2α2−1−α1λ2−α1, and it is obvious that T is invertible. Using the transformation unvn=Tx¯ny¯n, system (13) becomes(14)x¯n+1=−x¯n+f1un,vn,b∗,y¯n+1=λ2y¯n+g1un,vn,b∗,where the functions f1 and g1 denote the terms in model (14) in variables un,vn,b∗ with order at least two.Using the centre manifold theorem, we know that there exists a centre manifoldWc0,0,0 of system (14) at 0,0 in a small neighbourhood of b∗=0, which can be approximately described as follows:(15)Wc0,0,0=x¯n,y¯n,b∗εR3:y¯n+1=α¯1x¯n2+α¯2x¯nb∗+Ox¯n+b∗3,where(16)α¯1=α21+α1α11+α2β111−λ22+β221+α121−λ22−1+α1α121+α1+α2β121−λ22,α¯2=1+α1α231+α1+α2β23α21+λ22−1+α1α13+α2β131+λ22.System (14) which is restricted to the centre manifold Wc0,0,0 has the following form:(17)x¯n+1=−x¯n+h1x¯n2+h2x¯nb∗+h3x¯n2b∗+h4x¯nb∗2+h5x¯n3+Ox¯n+b∗3≡Fx¯n,b∗,where(18)h1=α¯2λ2−α¯1α11−α¯2β111+λ2−β221+α¯121+λ2−1+α¯1λ2−α¯1α12−α¯2β121+λ2,h2=λ2−α¯1α13−α¯2β131+λ2−1+α¯1λ2−α¯1α23−α¯2β23α¯21+λ2,h3=λ2−α1α¯1α13−α2β131+λ2+λ2−α1α23−α2β23λ2−α1α¯1α21+λ2−1+α1λ2−α1α123−α2β1231+λ2+α2λ2−α1α113−α2β1131+λ2−β2231+α121+λ2+2α2α¯2λ2−α1α11−α2β111+λ2−2β22α¯21+α1λ2−α11+λ2+α¯2λ2−α1α12−α2β12λ2−1−2α11+λ2,h4=α¯2λ2−α1α13−α2β131+λ2+λ2−α1α23−α2β23λ2−α1α¯2α21+λ2+2α2α¯2λ2−α1α11−α2β111+λ2+2β22α¯21+α1λ2−α11+λ2+α¯2λ2−α1α12−α2β12λ2−1−2α11+λ2,h5=2α2α¯1λ2−α1α11−α2β111+λ2+2β22α¯1λ2−α11+α11+λ2+λ2−α1α11−α2β11λ2−1−2α1α¯11+λ2+α¯22λ2−α1α111−α2β1111+λ2−α¯21+α1λ2−α1α112−α2β1121+λ2.For PDB, the two discriminatory quantitiesξ1 and ξ2 be nonzero,(19)ξ1=∂2F∂x¯∂b∗+12∂F∂b∗∂2F∂x¯2|0,0,ξ2=16∂3F∂x¯3+12∂2F∂x¯22|0,0.Finally, from the above analysis, we have the following result.Theorem 2.
Ifξ1≠0 and ξ2≠0, then model (3) undergoes PDB at P2x2,y2 when the parameter b varies in a small neighbourhood of bPDB. If ξ2>0ξ2<0, then the period-two orbits that bifurcate from P2x2,y2 are stable (unstable).


## 5. Influence of Refuge

This section presents the impact of prey refuge on the model when the coexistence steady-state pointP2x2,y2 exists and is stable. It is simple to observe that x and y are continuous differential functions of parameter b.Here,x2 is the positive root of the equation x2+Ax+B=0 where A=1−a/a+c/abd, B=−ce/bdd−1, and y2=x2/b−e/bd−1.Ifdx2/db>y2, then dy2/db>0; hence, y2 is strictly increasing function of the parameter b.Ifdx2/db<y2, then dy2/db<0; hence, y2 is strictly decreasing function of the parameter b.Clearly, if0<dx2/db<y2, then x2 is strictly increasing function of the parameter b and y2 is strictly decreasing function of the parameter b.Ifdx2/db<0, then x2 and y2 are strictly decreasing function of the parameter b.Ifdx2/db>y2, then dy2/db>0; hence, x2 and y2 are strictly increasing function of the parameter b.


## 6. Chaos Control

In discrete-time models, several approaches can be used in chaos control. These include the following: hybrid control method, pole-placement technique, and the state feedback method [9, 14, 18, 30]. In this regard, a feedback control method is employed to stabilize the chaotic orbits at an unstable positive fixed point of model (3) as follows:(20)xn+1=axn1−xn−cxn−bynyne+xn−byn+S,yn+1=dxn−bynyne+xn−byn.With the feedback control law as the control force,S=−q1xn−x2−q2yn−y2, where q1 and q2 are the feedback gain and x2,y2 is a positive fixed point of the model.The Jacobian matrixJ for system (20) at x2,y2 is J=a11−q1a12−q2a21a22, where a11=a1−2x−cey/e+x−by2, a12=−cex−2by+x−by2/e+x−by2, a21=dey/e+x−by2, and a22=dex−2by+x−by2/e+x−by2.The characteristic equation of the matrixJ is given by λ2−a11+a22−q1λ+a22a11−q1−a21a12−q2=0.Letλ1 and λ2 be the eigenvalues, then(21)λ1+λ2=a11+a22−q1,λ1λ2=a22a11−q1−a21a12−q2.The lines of marginal stability are determined by solving the equationλ1=±1 and λ1λ2=1. These conditions guarantee that the eigenvalues λ1 and λ2 have modulus less than 1.Supposeλ1λ2=1, we have a line l1 from (21) as a22q1−a21q2=a22a11−a21a12−1.Supposeλ1=±1, we have the line l2 and l3 from (21) as 1−a22q1+a21q2=a11+a22−1−a22a11+a21a12 and 1+a22q1−a21q2=a11+a22+1+a22a11−a21a12.The stable eigenvalues lie within a triangular region by the linesl1,l2, and l3.


## 7. Numerical Simulation

To support the analytical results, this section presents the numerical simulation based on system (3) through the following examples.Example 1.
The equilibrium points and their nature fora=4, b=0.4, c=3, d=3, and e=0.7 of system (3) are derived.
The discrete prey-predator system (3) has three equilibrium points. One is the trivial equilibrium point 0,0, another one is axial equilibrium point 0.75,0, and the last one is unique positive equilibrium point 0.53,0.46. The characteristic equation of the Jacobian matrix of system (3) evaluated at the trivial equilibrium point 0,0 is given by(22)Fλ=λ2−4λ=0.
Now,F1=−3<0, F−1=5>0, and F0=0<1. The roots of (22) are given by λ1=0 and λ2=4. Therefore, the trivial equilibrium point 0,0 is saddle.
The characteristic equation of the Jacobian matrix of system (3) evaluated at the axial equilibrium point 0.75,0 is given by(23)Fλ=λ2+0.4483λ−3.1034=0.
Now,F1=−1.6551<0, F−1=−2.5517<0, and F0=−3.1034<1. The roots of (23) are given by λ1=−2 and λ2=1.5517. Therefore, the axial equilibrium point 0.75,0 is source. The characteristic equation of the Jacobian matrix of system (3) evaluated at the positive equilibrium point 0.53,0.46 is given by(24)Fλ=λ2+0.4837λ−0.1534=0.
Now,F1=1.3303>0, F−1=0.3629>0, and F0=−0.1534<1, and then λ1<1 and λ2<1. The roots of (24) are given by λ1=−0.7022 and λ2=0.2185. Therefore, the positive equilibrium point 0.53,0.46 is sink.Example 2.
We consider the system parametersa=4, b∈0,1, c=3, d=3, e=0.7, and initial population 0.6,0.3. Then, system (3) undergoes the Neimark–Sacker bifurcation as b varies in a small neighbourhood of b=0.0355. The corresponding bifurcation diagrams and Lyapunov exponents are shown in Figures 1 and 2, respectively. For a=4, b=0.0355, c=3, d=3, and e=0.7, system (3) has three equilibrium points. One is the trivial equilibrium point 0,0, another one is axial equilibrium point 0.75,0, and the last one is unique positive equilibrium point 0.37,0.56. Now, we will study the dynamics of the interior equilibrium point 0.37,0.56. Furthermore, the characteristic equation of the Jacobian matrix of system (3) for a=4, b=0.0355, c=3, d=3, and e=0.7 evaluated at positive steady-state 0.37,0.56 is given by(25)λ2−0.9359λ+1.0009=0.
The roots of (25) are given by λ1,2=0.4680±0.8842i with λ1,2=1. Thus, ANSB=a=4,b=0.0355,c=3,d=3,e=0.7.Figure 1
Bifurcation diagrams fora=4,b∈0,1,c=3,d=3,e=0.7, and initial population 0.6,0.3. (a) Bifurcation diagram for prey x. (b) Bifurcation diagram for predator y.
(a)(b)Figure 2
Lyapunov exponents with respect to refuge parameterb.Example 3.
We consider parametersa=4, b∈0,1, c=3, d=3, e=0.7, and initial population 0.6,0.3. Then, system (3) undergoes the Period-Doubling Bifurcation as b varies in a small neighbourhood of b=0.52. The corresponding bifurcation diagram and Lyapunov exponents are shown in Figures 1 and 2, respectively. For a=4, b=0.52, c=3, d=3,ande=0.7, system (3) has three equilibrium points. One is the trivial equilibrium point 0,0, another one is axial equilibrium point 0.75,0, and the last one is unique positive equilibrium point 0.56,0.42. Now, we will study the dynamics of the interior equilibrium point 0.56,0.42. Furthermore, the characteristic equation of the Jacobian matrix of system (3) for a=4, b=0.52, c=3, d=3, and e=0.7 evaluated at positive steady-state 0.56,0.42 is given by(26)λ2+0.7318λ−0.2693=0.
The roots of (26) are given by λ1=−1 and λ2=0.269 with λ2≠1. Thus, APDB=a=4,b=0.52,c=3,d=3,e=0.7.
With the comparison of bifurcation analysis through the pictorial representation in Figure3 with Figure 1 for the same parameter value with only the difference in functional response parameter e, we observed that Neimark–Sacker bifurcation vanishes from Figure 3, but Period-Doubling Bifurcation occurs. Figure 4 is the corresponding Lyapunov exponent diagram of Figure 3. Furthermore, the bifurcation diagram by changing the parameter value a is shown in Figure 5, and compared with Figures 1 and 3, here both bifurcations are vanished.Figure 3
Bifurcation diagrams fora=4,b∈0,1,c=3,d=3,e=0.76, and initial population 0.6,0.3. (a) Bifurcation diagram for prey x. (b) Bifurcation diagram for predator y.
(a)(b)Figure 4
Lyapunov exponents with respect to refuge parameterb for a=4,b∈0,1,c=3,d=3,e=0.76, and initial population 0.6,0.3.Figure 5
Bifurcation diagrams fora=3,b∈0,1,c=3,d=3,e=0.7, and initial population 0.6,0.3. (a) Bifurcation diagram for prey x. (b) Bifurcation diagram for predator y.
(a)(b)Example 4.
For the parameter seta=4, b=0.95, c=3, d=3, e=0.75, and the initial population 0.6,0.3, we observe the chaotic behaviour of prey predator shown in Figure 6(a). In this case, the fixed point 0.64,0.28 is unstable. For feedback gain q1=−0.5 and q2=0.5, we observe that fixed point 0.64,0.28 is stable as shown in Figure 6(b). Model (3) with the feedback control law as the control force S=0.5xn−0.64−0.5yn−0.28 is as follows:(27)xn+1=axn1−xn−cxn−bynyne+xn−byn+S,yn+1=dxn−bynyne+xn−byn.
The Jacobian matrixJ for system (27) at 0.64,0.28 is J=−1.1187−1.02450.49870.5245. The characteristic equation of the matrix J is given by(28)λ2+0.5942λ−0.0759=0.
The roots of (28) are given by λ1=−0.7023 and λ2=0.1081. Now, line l1 is 0.5245q1−0.4987q2=−1.5874, line l2 is 0.4755q1+0.4987q2=−1.5068, and line l3 is 1.5245q1−0.4987q2=−0.6816.
The stable eigenvalues lie within a triangular region by the linesl1,l2, and l3 shown in Figure 7. The system is stable for the triangular region ABC bounded by the marginal lines l1, l2, and l3.Figure 6
Chaos control of the system for the parameter seta=4, b=0.95, c=3, d=3, and e=0.75. (a) Chaotic behaviour of system (3). (b) Stable behaviour of system (27) for feedback gain q1=−0.5 and q2=0.5.
(a)(b)Figure 7
Stability region of the controlled system (27) for the parameter set a=4, b=0.95, c=3, d=3, e=0.75, and feedback gain q1=−0.5 and q2=0.5.


## 8. Conclusion

In this paper, the impacts of refuges used by prey on a prey-predator interaction are studied through the analytical approach and graphical representation. This study considered the prey refuge proportional to predator density due to fear induced by a predator on prey. The detailed effects of the local stability of the interior steady-state point and the long-term dynamics of the interacting populations on the model have been investigated. The results show that the impact of refuges can stabilize and destabilize the interior equilibrium point of the proposed prey-predator model for different parameter sets.Different observations have been made on the types of refuge and their impacts on the proposed system. Generally, in the existing literature, we always found the stabilizing effect of refuge. However, here, we observed an unstable to stable to unstable behaviour in Figure1. For a different parameter set in Figure 3, we have also found a stable to unstable behaviour of the system. In Figure 5, another interesting observation is that refuge is profitable for prey and adverse effect on the predator. We have studied the bifurcations of the proposed discrete prey-predator model with refuge. It is shown that the model exhibits various bifurcations of codimension one, including PDB and NSB, as the values of the parameters vary.



---

*Source: 6209908-2021-10-25.xml*