# Dynamical Complexities of a Discrete Ivlev-Type Predator-Prey System

Authors: Sarker Md. Sohel Rana
Journal: Discrete Dynamics in Nature and Society
 (2023)

---

## Abstract

In this study, we examine a discrete predator-prey system from the following two perspectives: (i) the functional response is of the Ivlev type and (ii) the prey growth rate is of the Gompertz type. We define the stability requirement for feasible fixed points. We demonstrate algebraically that if the bifurcation (control) parameter rises over its threshold value, the system encounters flip and Neimark–Sacker (NS) bifurcations in the vicinity of the interior fixed point. We explicitly establish the existence requirements and direction of bifurcations via the center manifold theory. Analytical findings are validated by numerical simulations, which are used to highlight the occurrence of instability and chaotic dynamics in the system. In order to regulate the chaotic trajectories that exist in the system, we adopt a feedback control approach.


---

## Body

## 1. Introduction

In population dynamics, the central goal is to understand and analyze the interaction between predator and prey species due to their universal existence and importance. Several models have been developed to describe the dynamic relationship between predators and their prey using different ecological conditions. Prey growth rate and prey loss rate owing to predation are both included in the mathematical structure of predator-prey systems, which has been progressively studied through field and laboratory studies. Prey loss rate as a result of predation is defined as the functional response of a predator, which is the change in the prey density attached per unit time per predator as the prey density fluctuates. A traditional Gauss-type predator-prey system is represented by the following equation:(1)x˙=xGx,K−fxy,y˙=βfx−dy,where x is the prey (victim) density and y is the predator density, respectively? in the absence of a predator, Gx,K stands for a specific rate of prey growth; The positive constants K,β,d represent the prey’s carrying capacity, the rate at which a prey is converted into a predator, and the rate at which a predator dies, respectively. The rate at which prey are captured per predator, or the predator’s functional reaction, is fx. The majority of examples in the literature assume that(2)H:fx∈C10,∞,f0=0,f′x>0for allx>0andlimx⟶∞fx=m<∞.In the ecological literature, the most frequently investigated mathematical model of predator-prey interaction is system (1) with conventional logistic type prey growth rate Gx,K=r1−x/K,r>0 stands for prey intrinsic growth rate and Holling type II predator functional response fx=x/a+x,a>0 is the half-saturation constant [1]. For various ecological models, various functional responses have been employed. But according to research, nonmonotonic responses take place when microbial dynamics are inhibited and population dynamics are defended by groups. For several other categories of functional reactions, we refer to [1–3].Gompertz [4] created an alternate equation for the prey’s growth that is comparable in effect to logistic growth, gx,K=rxlnK/x, to study the dynamics of a community made up of populations of several interacting species. On the other hand, Ivlev [5] proposed a different functional response based on empirical evidence, known as the Ivlev functional response: px=α1−e−ax, where α and a are positive constants and stand for the maximum rate of predation and the predator’s effectiveness in capturing prey, respectively. Additionally, it supports hypothesis H. Figure 1 demonstrates the comparison: growth curves and functional responses. We find that the Gompertz curve grows faster than the logistic curve. Also, the Gompertz curve achieves carrying capacity a little bit earlier than the logistic curve, and its point of inflection happens sooner than that of the logistic curve. Moreover, the predation rate achieves its peak substantially sooner in the Ivlev-type functional response than in the Holling-type II functional response. In terms of biology, this implies that the predation rate saturates to a constant α when the prey population is large and is proportionate to the prey population when the prey population is low. Wide-ranging applications of Ivlev-type predator-prey interactions may be found in ecology. Examples include dynamics in the predator-prey system [6–13], host-parasitoid systems [14, 15], and animal coat patterns [16]. The authors investigated the presence and uniqueness of limit cycles as well as the numerical computation of phase portraits in these empirical works.Although there are numerous and comprehensive results on the dynamics of continuous predator-prey systems describing the interactions between populations with generations overlapping life cycles, many explanatory research works have suggested that discrete-time models based on difference equations [17] are better suited to characterize the living situations where the birth occurs during a regular mating season in populations with nonoverlapping generations, such as insects, birds, fish, the phytoplankton–zooplankton reproduction cycle, and plants herbivores. The analysis demonstrates that the discrete-time model has distinct features and a different structure from the continuous one. The discretized system exhibits richer and more complex dynamics including periodicity, quasiperiodicity, and chaos with respect to different parameters. Additionally, we point out that the algorithm for the numerical simulations is automatically provided by the discrete version of the continuous system. Therefore, the discrete form is a logical link between the simulations and the real model. On the other hand, mathematical models are frequently dependent on biological and chemical facts, which are frequently seen and gathered at discrete times [18]. In addition, the low-dimensional differential equation models (1.1) do not capture the dynamics of predators and prey, such as chaos. In reality, chaos does exist in nature, not just in mathematics. Examples include ecosystems and chemical reactions [19–21]. More significantly, discrete-time models are comparably simpler and less computationally intensive [22]. Because of this, compared to continuous forms, discrete-time model formulation and simulation are typically more straightforward, practical, and accurate. In recent years, the predator-prey relationship expressed in discrete-time form has attracted much attention [23–34]. These studies examined both monotonic and nonmonotonic functional responses as well as the conventional logistic-type of prey growth rate. In these research efforts, the possibilities of bifurcations, periodic orbits, and chaotic sets, all of which had been computed numerically or by means of center manifold theory, were all investigated.Taking into consideration the above research works, we express our interest in looking at system (1) with the Ivlev-type functional response of the predator and Gompertz growth on the prey population [6].(3)x˙=rxlnKx−α1−e−axy,y˙=β1−e−axy−dy,where x and y stand densities of prey and predator, respectively; r,K,a,d are all positive constants. The prey’s intrinsic growth rate is represented by parameter r, while its carrying capacity is represented by parameter K. Applying the forward Euler scheme, the discrete form of (3) is as follows:(4)xy↦x+δrxlnKx−α1−e−axyy+δβ1−e−axy−dy,where δ is the integral step size. From a biological standpoint, assume that the system is defined in the region R+2=x,y:x≥0,y≥0 in the x,y-plane. A recent study in the literature [35, 36] looked at a discrete Ivlev-type predator-prey system. For instance, in [35, 36], the authors examined the predator-prey system using the traditional logistic type of prey individual growth, whereas in [37], Gompertz growth of prey individuals with Michaelis–Menten functional response was taken into account. Through the use of the center manifold theory, all of these works focused on determining the stability conditions of flip and NS bifurcations. The task is to identify how parameters influence the dynamics of system (4). We specifically use center manifold theory [38] to systematically examine the criterion for the existence and direction of bifurcations, namely, flip and NS bifurcations in R+2. Flip bifurcation and NS bifurcation are both significant mechanisms for the formation of complex dynamics in the discrete system, and both bifurcations open a path to chaos through periodic as well as quasiperiodic orbits.The structure of this study is as follows: Section2 provides the stability analysis for feasible fixed points of system (4) in the interior of R+2. In Section 3, we identify the stability conditions for the bifurcation of system (4) under specific parametric spaces. By altering the values of the control parameters, the system dynamics are numerically displayed in Section 4. A feedback control method has been employed in Section 5 to stabilize unstable trajectories. Section 6 provides a succinct conclusion.Figure 1
Comparison: growth curves and functional responses: (a) growth curves; (b) functional responses forr=0.5,K=1.65,a=1.2,α=0.5,β=0.65,d=0.45,P0=0.05.
(a)(b)Figure 2
Flip bifurcation diagrams of system (4) for δ: (a) in prey, (b) in predator, (c) MLEs, and (d) FD.
(a)(b)(c)(d)


## 2. Analysis of the Stability and Existence Criteria for Fixed Points

### 2.1. Study of Fixed Points

According to a straightforward algebraic calculation, system (4) has the following two fixed points for all possible parameter values:(1)
E1K,0, the predator-free fixed point(2)
E2x∗,y∗, the singular fixed point if d<1−e−aKβ where x∗=−1/aln1−d/β and y∗=βr/αdx∗lnK/x∗Biologically,E1K,0 interprets that whenever no predators in the system, the prey population reaches its carrying capacity whereas E2x∗,y∗ stands for the coexistence of both populations.


### 2.2. Dynamical Behavior: Stability Analysis

When examining the stability of system (4) at fixed points Ex,y, it is essential to bear in mind that the absolute value of the eigenvalues of the Jacobian matrix derived at the point influences the local stability of that fixed point. Given at Ex,y, the Jacobian matrix of system (4) takes the form(5)Jx,y=j11j12j21j22,where(6)j11=1+δ−r−ae−axyα+rlnKx,j12=α−1+ae−axδ,j21=aβe−axyδ,j22=1+−d+1−e−axβδ.The characteristic equation of a matrixJ is as follows:(7)λ2−trJx,yλ+detJx,y=0.Now, using Jury’s criteria [39], the topological categorization of stability around fixed points is represented as follows:Proposition 1.
For the predator-free fixed pointE1K,0, the topological categorization listed below is accurate as follows:(i)
Ifd−1−e−aKβ>0 then E1 is a(1)
Sink if0<δ<min2/r,2/d−1−e−aKβ(2)
Source ifδ>max2/r,2/d−1−e−aKβ(3)
Nonhyperbolic ifδ=2/ror2/d−1−e−aKβ(ii)
Ifd−1−e−aKβ<0 then E1 is a(1)
Source ifδ>2/r(2)
Saddle ifδ<2/r(3)
Nonhyperbolic ifδ=2/r(iii)
Ifd−1−e−aKβ=0 then E1 is always nonhyperbolic.It is obvious that whenδ=2/ror2/d−1−e−aKβ, then one eigenvalue of JE1 is −1 and the other is not equal to ±1. Thus, system (4) experiences a flip bifurcation if parameters change in a small vicinity of FBE11 or FBE12.(8)FBE11=r,K,a,d,α,β,δ∈R+2:δ=2r,δ≠2d−1−e−aKβ,d−1−e−aKβ>0.or(9)FBE12=r,K,a,d,α,β,δ∈R+2:δ=2d−1−e−aKβ,δ≠2r,d−1−e−aKβ>0.AtE2x∗,y∗, we write (7) as(10)Fλ≔=λ2−2+Mδλ+1+Mδ+Nδ2=0,where(11)M=−d−r−ae−ax∗y∗α+1−e−ax∗β+rlnKx∗,N=−ae−ax∗−1+e−ax∗y∗αβ+−d+1−e−ax∗β−r−ae−ax∗y∗α+rlnKx∗.Therefore,F1=Nδ2>0 and F−1=4+2Mδ+Nδ2. Let L=M2−4N. For the topological classification of E2, we state the following Proposition.Proposition 2.
Suppose thatd−1−e−aKβ<0. Then, a fixed point E2x∗,y∗ of system (4) exists and the topological categorization listed as follows is accurate:(i)
If any of the following applies,E2 is a sink(1)
L≥0andδ<−M−L/N(2)
L<0andδ<−M/N(ii)
If any of the following applies,E2 is a source(1)
L≥0andδ>−M−L/N(2)
L<0andδ>−M/N(iii)
If any of the following applies,E2 is a nonhyperbolic(1)
L≥0andδ=−M−L/N;δ≠−2/M,−4/M(2)
L<0andδ=−M/N(iv)
If otherwise,E2 is a saddleThe eigenvalues ofJE2 are λ1=−1 and λ2≠±1 if condition (iii.1) of Proposition 2 is valid and if (iii.2) is true, two eigenvalues λ1,2 are complex having magnitude one.Let(12)FBE21=r,K,a,d,α,β,δ∈R+2:δ=−M−L/N,L≥0,δ≠−2/M,−4/M,or(13)FBE22=r,K,a,d,α,β,δ∈R+2:δ=−M−L/N,L≥0,δ≠−2/M,−4/M.Then, system (4) undergoes a flip bifurcation at E2 if parameters r,K,a,d,α,β,δ∈FBE21 or FBE22.Also, let(14)NSBE2=r,K,a,d,α,β,δ∈R+2:δ=−M/N,L<0.Therefore, an NS bifurcation occurs atE2 in system (4) if the parameters r,K,a,d,α,β,δ∈NSBE2.


## 2.1. Study of Fixed Points

According to a straightforward algebraic calculation, system (4) has the following two fixed points for all possible parameter values:(1)
E1K,0, the predator-free fixed point(2)
E2x∗,y∗, the singular fixed point if d<1−e−aKβ where x∗=−1/aln1−d/β and y∗=βr/αdx∗lnK/x∗Biologically,E1K,0 interprets that whenever no predators in the system, the prey population reaches its carrying capacity whereas E2x∗,y∗ stands for the coexistence of both populations.


## 2.2. Dynamical Behavior: Stability Analysis

When examining the stability of system (4) at fixed points Ex,y, it is essential to bear in mind that the absolute value of the eigenvalues of the Jacobian matrix derived at the point influences the local stability of that fixed point. Given at Ex,y, the Jacobian matrix of system (4) takes the form(5)Jx,y=j11j12j21j22,where(6)j11=1+δ−r−ae−axyα+rlnKx,j12=α−1+ae−axδ,j21=aβe−axyδ,j22=1+−d+1−e−axβδ.The characteristic equation of a matrixJ is as follows:(7)λ2−trJx,yλ+detJx,y=0.Now, using Jury’s criteria [39], the topological categorization of stability around fixed points is represented as follows:Proposition 1.
For the predator-free fixed pointE1K,0, the topological categorization listed below is accurate as follows:(i)
Ifd−1−e−aKβ>0 then E1 is a(1)
Sink if0<δ<min2/r,2/d−1−e−aKβ(2)
Source ifδ>max2/r,2/d−1−e−aKβ(3)
Nonhyperbolic ifδ=2/ror2/d−1−e−aKβ(ii)
Ifd−1−e−aKβ<0 then E1 is a(1)
Source ifδ>2/r(2)
Saddle ifδ<2/r(3)
Nonhyperbolic ifδ=2/r(iii)
Ifd−1−e−aKβ=0 then E1 is always nonhyperbolic.It is obvious that whenδ=2/ror2/d−1−e−aKβ, then one eigenvalue of JE1 is −1 and the other is not equal to ±1. Thus, system (4) experiences a flip bifurcation if parameters change in a small vicinity of FBE11 or FBE12.(8)FBE11=r,K,a,d,α,β,δ∈R+2:δ=2r,δ≠2d−1−e−aKβ,d−1−e−aKβ>0.or(9)FBE12=r,K,a,d,α,β,δ∈R+2:δ=2d−1−e−aKβ,δ≠2r,d−1−e−aKβ>0.AtE2x∗,y∗, we write (7) as(10)Fλ≔=λ2−2+Mδλ+1+Mδ+Nδ2=0,where(11)M=−d−r−ae−ax∗y∗α+1−e−ax∗β+rlnKx∗,N=−ae−ax∗−1+e−ax∗y∗αβ+−d+1−e−ax∗β−r−ae−ax∗y∗α+rlnKx∗.Therefore,F1=Nδ2>0 and F−1=4+2Mδ+Nδ2. Let L=M2−4N. For the topological classification of E2, we state the following Proposition.Proposition 2.
Suppose thatd−1−e−aKβ<0. Then, a fixed point E2x∗,y∗ of system (4) exists and the topological categorization listed as follows is accurate:(i)
If any of the following applies,E2 is a sink(1)
L≥0andδ<−M−L/N(2)
L<0andδ<−M/N(ii)
If any of the following applies,E2 is a source(1)
L≥0andδ>−M−L/N(2)
L<0andδ>−M/N(iii)
If any of the following applies,E2 is a nonhyperbolic(1)
L≥0andδ=−M−L/N;δ≠−2/M,−4/M(2)
L<0andδ=−M/N(iv)
If otherwise,E2 is a saddleThe eigenvalues ofJE2 are λ1=−1 and λ2≠±1 if condition (iii.1) of Proposition 2 is valid and if (iii.2) is true, two eigenvalues λ1,2 are complex having magnitude one.Let(12)FBE21=r,K,a,d,α,β,δ∈R+2:δ=−M−L/N,L≥0,δ≠−2/M,−4/M,or(13)FBE22=r,K,a,d,α,β,δ∈R+2:δ=−M−L/N,L≥0,δ≠−2/M,−4/M.Then, system (4) undergoes a flip bifurcation at E2 if parameters r,K,a,d,α,β,δ∈FBE21 or FBE22.Also, let(14)NSBE2=r,K,a,d,α,β,δ∈R+2:δ=−M/N,L<0.Therefore, an NS bifurcation occurs atE2 in system (4) if the parameters r,K,a,d,α,β,δ∈NSBE2.


## 3. Analysis of Bifurcations: Direction and Stability

This part will focus to recapitulate the conditions for the stability and direction of flip and NS bifurcations of system (4) around E2 by application of the bifurcation theory [38]. We take δ as a bifurcation parameter, otherwise stated.


### 3.1. Flip Bifurcation: Stability and Direction

Choose parametersr,K,a,d,α,β,δ∈FBE21 for system (4) around E2x∗,y∗. One can choose parameters in FBE22 to do the analysis. Let δ=δF=−M−L/N, then JE2 returns two eigenvalues λ1δF=−1 and λ2δF=3+MδF. The condition λ2δF≠1 leads to the following equation:(15)δF≠−2M,−4M.We setx˜=x−x∗,y˜=y−y∗. Then, system (4) reduces to(16)x˜y˜⟶Aδx˜y˜+F1x˜,y˜,δF2x˜,y˜,δ,where Aδ=Jx∗,y∗ and(17)F1x˜,y˜,δ=−16a2e−ax∗x˜2βδ−3y˜+ax˜y∗+rx˜3δ6x∗2−rx˜2δ2x∗+12ae−ax∗x˜βδ−2y˜+ax˜y∗+OX4,F2x˜,y˜,δ=16a2e−ax∗x˜2βδ−3y˜+ax˜y∗−12ae−ax∗x˜βδ−2y˜+ax˜y∗+OX4.It is possible to express system (16) as follows:(18)Xn+1=AXn+12BXn,Xn+16CXn,Xn,Xn+OXn4,X=x˜,y˜T,where Bx,y=B1x,yB2x,y and x,y,u=C1x,y,uC2x,y,u are symmetric multilinear vector functions of x,y,u∈R2 and defined as follows:(19)B1x,y=∑j,k=12δ2F1ξ,δδξjδξkξ=0xjyk=−rx1y1δx∗+ae−ax∗αδ−x2y1−x1y2+ax1y1y∗,B2x,y=∑j,k=12δ2F2ξ,δδξjδξkξ=0xjyk=−ae−ax∗αδ−x2y1−x1y2+ax1y1y∗,C1x,y,u=∑j,k,l=12δ2F1ξ,δδξjδξkδξlξ=0xjykul=−a2e−ax∗αδ−u2x1y1−u1x2y1+x1y2−ay1y∗+ru1x1y1δx∗2,C2x,y,u=∑j,k,l=12δ2F2ξ,δδξjδξkδξlξ=0xjykul=a2e−ax∗δ−u2x1y1−u1x2y1+x1y2−ay1y∗,and δ=δF.The direct calculation yields the following two eigenvectorsp,q∈R2 of A for eigenvalue λ1δF=−1 satisfying AδFq=−q and ATδFp=−p.(20)q∼2−dδF+β1−e−ax∗δF,−ae−ax∗βδFy∗T,p∼2−dδF+β1−e−ax∗δF,1−e−ax∗αδFT.InR2, we employ the conventional scalar product defined by p,q=p1q1+p2q2, to normalize p,q such that p,q=1. To do this, we set p=γF2−dδF+β1−e−ax∗δF,1−e−ax∗αδFT, where(21)γF=12−dδF+β1−e−ax∗δF2−ae−ax∗αβδF2y∗1−e−ax∗.Then, the coefficient of the normal form is(22)l1δF=16p,Cq,q,q−12p,Bq,A−I−1Bq,q.According to the study above, we present the following conclusion regarding the stability and direction of the flip bifurcation:Theorem 1.
Assume that (15) is true for the fixed point E2x∗,y∗. Then, system (4) encounters a flip bifurcation at E2 if l1δF≠0 and δ fluctuate its value in a limited proximity of FBE21. Moreover, stable (resp., unstable) period-2 orbits split off from E2 if l1δF is positive (resp., negative).


### 3.2. NS Bifurcation: Stability and Direction

Next, we take into account system (4) at E2x∗,y∗ with r,K,a,d,α,β,δ∈NSBE2. From equation (7), the complex roots (eigenvalues) are provided by the following equation:(23)λ,λ¯=1+Mδ2±iδ24N−M2.Let(24)δ=δNS=−MN.Therefore, it follows thatλ=1+Mδ+Nδ2 with 1+Mδ+Nδ2δ=δNS=1. The condition of the transversality condition yields(25)dλδdδδ=δNS=−M2≠0.Moreover, the nondegenerate conditionpδNS≠0,1, obviously satisfies(26)M2N≠2,3,and we have(27)λkδNS≠1fork=1,2,3,4.By direct calculation, we obtain the following two eigenvectorsq,p∈C2 of AδNS and ATδNS for eigenvalues λδNS and λ¯δNS such that AδNSq=λδNSq,ATδNSp=λ¯δNSp.(28)q∼1−dδNS+β1−e−ax∗δNS−λ,−ae−ax∗βδNSy∗T,p∼1−dδNS+β1−e−ax∗δNS−λ¯,1−e−ax∗αδNST.Forp,q=1, we set p=γNS1−dδNS+β1−e−ax∗δNS−λ¯,1−e−ax∗αδNST, where γNS=1/1−dδNS+β1−e−ax∗δNS−λ¯2−ae−ax∗αβδNS2y∗1−e−ax∗.Now, decomposing the vectorX∈R2 as X=wq+w¯q¯, for δ close to δNS and w∈C. Obviously, w=<p,X>. As a result, for δ near δNS, we derive the transformed form of system (16) as follows:(29)w↦λδw+gw,w¯,δ,where λδ=1+ϕδeiθδ with ϕδNS=0 and gw,w¯,δ is a smooth complex-valued function. Then, we obtain(30)gw,w¯,δ=∑k+l≥21k!l!gklδwkw¯l,withgkl∈C,k,l=0,1,The coefficientsgkl are determined via multilinear symmetric vector functions.(31)g20δNS=<p,Bq,q>,g11δNS=<p,Bq,q¯>,g02δNS=<p,Bq¯,q¯>,g21δNS=<p,Cq,q,q¯>.The coefficient of the critical normal form is as follows:(32)l2δNS=Ree−iθδNSg212−Re1−2eiθδNSe−2iθδNS21−eiθδNSg20g11−12g112−14g022,where eiθδNS=λδNS. As a conclusion to the aforementioned study, we provide the following theorem on the stability and direction of the NS bifurcation.Theorem 2.
Assume that (26) is valid and l2δNS≠0. System (4) encounters a Neimark–Sacker bifurcation at E2x∗,y∗ if δ fluctuates its value in a limited vicinity of NSBE2. Moreover, a singular invariant closed curve bifurcates from E2 that is attracting (resp., repelling) if l2δNS is negative (resp., positive), and the NS bifurcation is supercritical (resp., subcritical).


## 3.1. Flip Bifurcation: Stability and Direction

Choose parametersr,K,a,d,α,β,δ∈FBE21 for system (4) around E2x∗,y∗. One can choose parameters in FBE22 to do the analysis. Let δ=δF=−M−L/N, then JE2 returns two eigenvalues λ1δF=−1 and λ2δF=3+MδF. The condition λ2δF≠1 leads to the following equation:(15)δF≠−2M,−4M.We setx˜=x−x∗,y˜=y−y∗. Then, system (4) reduces to(16)x˜y˜⟶Aδx˜y˜+F1x˜,y˜,δF2x˜,y˜,δ,where Aδ=Jx∗,y∗ and(17)F1x˜,y˜,δ=−16a2e−ax∗x˜2βδ−3y˜+ax˜y∗+rx˜3δ6x∗2−rx˜2δ2x∗+12ae−ax∗x˜βδ−2y˜+ax˜y∗+OX4,F2x˜,y˜,δ=16a2e−ax∗x˜2βδ−3y˜+ax˜y∗−12ae−ax∗x˜βδ−2y˜+ax˜y∗+OX4.It is possible to express system (16) as follows:(18)Xn+1=AXn+12BXn,Xn+16CXn,Xn,Xn+OXn4,X=x˜,y˜T,where Bx,y=B1x,yB2x,y and x,y,u=C1x,y,uC2x,y,u are symmetric multilinear vector functions of x,y,u∈R2 and defined as follows:(19)B1x,y=∑j,k=12δ2F1ξ,δδξjδξkξ=0xjyk=−rx1y1δx∗+ae−ax∗αδ−x2y1−x1y2+ax1y1y∗,B2x,y=∑j,k=12δ2F2ξ,δδξjδξkξ=0xjyk=−ae−ax∗αδ−x2y1−x1y2+ax1y1y∗,C1x,y,u=∑j,k,l=12δ2F1ξ,δδξjδξkδξlξ=0xjykul=−a2e−ax∗αδ−u2x1y1−u1x2y1+x1y2−ay1y∗+ru1x1y1δx∗2,C2x,y,u=∑j,k,l=12δ2F2ξ,δδξjδξkδξlξ=0xjykul=a2e−ax∗δ−u2x1y1−u1x2y1+x1y2−ay1y∗,and δ=δF.The direct calculation yields the following two eigenvectorsp,q∈R2 of A for eigenvalue λ1δF=−1 satisfying AδFq=−q and ATδFp=−p.(20)q∼2−dδF+β1−e−ax∗δF,−ae−ax∗βδFy∗T,p∼2−dδF+β1−e−ax∗δF,1−e−ax∗αδFT.InR2, we employ the conventional scalar product defined by p,q=p1q1+p2q2, to normalize p,q such that p,q=1. To do this, we set p=γF2−dδF+β1−e−ax∗δF,1−e−ax∗αδFT, where(21)γF=12−dδF+β1−e−ax∗δF2−ae−ax∗αβδF2y∗1−e−ax∗.Then, the coefficient of the normal form is(22)l1δF=16p,Cq,q,q−12p,Bq,A−I−1Bq,q.According to the study above, we present the following conclusion regarding the stability and direction of the flip bifurcation:Theorem 1.
Assume that (15) is true for the fixed point E2x∗,y∗. Then, system (4) encounters a flip bifurcation at E2 if l1δF≠0 and δ fluctuate its value in a limited proximity of FBE21. Moreover, stable (resp., unstable) period-2 orbits split off from E2 if l1δF is positive (resp., negative).


## 3.2. NS Bifurcation: Stability and Direction

Next, we take into account system (4) at E2x∗,y∗ with r,K,a,d,α,β,δ∈NSBE2. From equation (7), the complex roots (eigenvalues) are provided by the following equation:(23)λ,λ¯=1+Mδ2±iδ24N−M2.Let(24)δ=δNS=−MN.Therefore, it follows thatλ=1+Mδ+Nδ2 with 1+Mδ+Nδ2δ=δNS=1. The condition of the transversality condition yields(25)dλδdδδ=δNS=−M2≠0.Moreover, the nondegenerate conditionpδNS≠0,1, obviously satisfies(26)M2N≠2,3,and we have(27)λkδNS≠1fork=1,2,3,4.By direct calculation, we obtain the following two eigenvectorsq,p∈C2 of AδNS and ATδNS for eigenvalues λδNS and λ¯δNS such that AδNSq=λδNSq,ATδNSp=λ¯δNSp.(28)q∼1−dδNS+β1−e−ax∗δNS−λ,−ae−ax∗βδNSy∗T,p∼1−dδNS+β1−e−ax∗δNS−λ¯,1−e−ax∗αδNST.Forp,q=1, we set p=γNS1−dδNS+β1−e−ax∗δNS−λ¯,1−e−ax∗αδNST, where γNS=1/1−dδNS+β1−e−ax∗δNS−λ¯2−ae−ax∗αβδNS2y∗1−e−ax∗.Now, decomposing the vectorX∈R2 as X=wq+w¯q¯, for δ close to δNS and w∈C. Obviously, w=<p,X>. As a result, for δ near δNS, we derive the transformed form of system (16) as follows:(29)w↦λδw+gw,w¯,δ,where λδ=1+ϕδeiθδ with ϕδNS=0 and gw,w¯,δ is a smooth complex-valued function. Then, we obtain(30)gw,w¯,δ=∑k+l≥21k!l!gklδwkw¯l,withgkl∈C,k,l=0,1,The coefficientsgkl are determined via multilinear symmetric vector functions.(31)g20δNS=<p,Bq,q>,g11δNS=<p,Bq,q¯>,g02δNS=<p,Bq¯,q¯>,g21δNS=<p,Cq,q,q¯>.The coefficient of the critical normal form is as follows:(32)l2δNS=Ree−iθδNSg212−Re1−2eiθδNSe−2iθδNS21−eiθδNSg20g11−12g112−14g022,where eiθδNS=λδNS. As a conclusion to the aforementioned study, we provide the following theorem on the stability and direction of the NS bifurcation.Theorem 2.
Assume that (26) is valid and l2δNS≠0. System (4) encounters a Neimark–Sacker bifurcation at E2x∗,y∗ if δ fluctuates its value in a limited vicinity of NSBE2. Moreover, a singular invariant closed curve bifurcates from E2 that is attracting (resp., repelling) if l2δNS is negative (resp., positive), and the NS bifurcation is supercritical (resp., subcritical).


## 4. Numerical Simulations

In order to verify the validity of our theoretical findings, the dynamics of system (4) around the fixed point E2 are investigated numerically. We consider several sets of the parameter values listed in Table 1 for the bifurcation analysis.Example 1.
The numerical simulation for flip bifurcation with respect to parameterδ.
Taking parameter values as given in case (i) withx0,y0=0.915,5.53, system (4) experiences a flip bifurcation (see Figures 2(a) and 2(b)) around E20.928774,5.56516 when δ crosses its threshold value δF∼0.918401. Also, the eigenvalues are λ1=−1,λ2=0.961275, l1δF=4.22783 and r,K,a,d,α,β,δ∈FBE21. This validates the results in Theorem 1. Figure 2 reveals that there exists a period-doubling phenomenon of periods 2, 4, 8, and 16 orbits at δ∼0.925,1.05,1.078,and1.082 in the range δ∈0.75,1.084 trigger chaotic set for δ∈1.084,1.135 and the period 6 orbit at δ∼1.113 which is in chaotic range. The MLEs and FD associated with Figures 2(a) and 2(b) are plotted in Figures 2(c) and 2(d), which justify the periodic and chaotic states that exist in system (4). The phase portrait diagrams are displayed in Figure 3.Table 1
Parameter values.
CasesVarying parameterFixed parametersDynamics of system (4)(i)0.75≤δ≤1.135r=2.5,K=1.5,a=0.55,α=0.45,β=0.25,d=0.1Flip bifurcation(ii)0.1≤δ≤0.785r=0.5,K=1.65,a=4.6,α=0.5,β=0.65,d=0.45NS bifurcation(iii)0.45≤β≤0.65r=0.5,K=1.65,a=4.6,α=0.5,d=0.45,δ=0.785NS bifurcationFigure 3
Phase portraits (inxy-plane) for different values of δ associated with Figures 2(a) and 2(b).Example 2.
The numerical simulation for NS bifurcation with respect to parameterδ
By taking parameters as given in case (ii) withx0,y0=0.253,0.683 and after calculation, we observe that system (4) encounters a NS bifurcation at E20.256229,0.689312 when δ crosses its threshold value through δNS∼0.257801. The eigenvalues are λ,λ¯=0.992705±0.120567i, dλδ/dδδ=δNS=0.0282961, and M2/N=0.0145895. Also, we have g20=−0.00317817−0.0334566i,g11=0.0352861+0.00423371i,g02=−0.0599838+0.0318261i,g21=−0.0140495+0.011453i, l2δNS=−0.00305513<0 and r,K,a,d,α,β,δ∈NSBE2. The correctness of Theorem 2 is ascertained. The bifurcation diagrams are depicted in Figures 4(a) and 4(b), which illustrate that there exists an attracting invariant closed cycle for δ>δNS and as δ increases, the disappearance of the closed curve occurs suddenly through periods 27, 80, 53, 26, and 52 orbits leading to chaos in system (4). The MLEs and FD connected to Figures 4(a) and 4(b) are disposed of in Figures 4(c) and 4(d), which illustrates the presence of periodic orbits and chaos as parameter δ rises. The phase portraits of system (4) associated with Figures 4(a) and 4(b) are plotted in Figure 5 illustrating the act of the smooth invariant curve and how it bifurcates from the stable fixed point and increases its radius. System (4) possesses a sub-NS bifurcation and a flip bifurcation with a rise in the δ value.Figure 4
NS bifurcation diagrams of system (4) for δ: (a) in prey, (b) local amplification in prey, (c) MLEs, and (d) FD.
(a)(b)(c)(d)Example 3.
The numerical simulation for the NS bifurcation of system (4) with varying parameter β
Fixingδ=0.785 with parameters as given in case (ii) and considering β as the bifurcation parameter, we discover that system (4) experiences a NS bifurcation at E20.491437,0.664528 when β crosses its threshold value through βNS∼0.502393. Figure 6 depicts the bifurcation diagrams, MLEs, and FD, respectively, which illustrate that an attracting invariant closed cycle appears for β>βNS and on the route to chaos, period 26 and period 52 orbits (at β∼0.618and0.644) exist (see Figure 7).Figure 5
Phase portraits (inxy-plane) for different values of δ connected to Figures 4(a) and 4(b).Figure 6
NS bifurcation diagrams of system (4) for β: (a) in prey, (b) local amplification in prey, (c) MLEs, and (d) FD x0,y0=0.2462,0.6793.
(a)(b)(c)(d)Figure 7
Phase portraits (inxy-plane) associated with Figures 6(a) and 6(b).Example 4.
Dynamics of system (4) with respect to two control parameters
When two additional parameters go over their critical thresholds, the dynamics of system (4) might become more complex. The plots of the 2D projection of 3D MLEs for two parameters are shown in Figure 8. This graphic serves as a tool to see how the system dynamics alter qualitatively as parameter values rise. Finding control parameter values for which system (4) dynamics are nonchaotic, periodic, or chaotic is straightforward. We vary two parameters in a certain range and fix the rest parameters as in case (ii). Therefore, we first plot (Figure 8(a)) the 2D projection of 3D MLEs for δ,β∈0.2,0.785×0.45,0.65. Figure 8(b) is the plot in δ,a plane with δ,a∈0.2,0.785×3.0,4.8. Figure 8(c) is the plot in β,a plane with β,a∈0.45,0.65×4.0,4.8. From these graphics, bifurcation parameter values may now be easily identified to observe how does the dynamics of system (4) switch from static to periodic or chaotic and we observe the following:(i)
The dynamics of system (4) changes from being nonperiodic to having an attractive fixed point or stable periodic cycle when the control parameters δ and β values grow(ii)
As control parameter valuesδ and a rise, the dynamics of system (4) shift from nonchaotic to chaotic states(iii)
The system dynamics remarkably shift from static to chaotic when the values of the control parametersβ and a grow
For example, the nonchaotic dynamics exist forδ=0.22,a=4.6 and the chaotic dynamics exist for δ=0.785,a=4.6 (see Figure 5), which are both consistent with the signs of the MLEs shown in Figure 8(c).Figure 8
Dynamics of system (4) in 2D parameter space. The 2D projection of 3D MLEs onto (a)δ,β plane, (b)δ,a plane, and (c)δ,β plane.
(a)(b)(c)


### 4.1. Fractal Dimension of System (4)

Strange attractors in a system are characterized by their fractal dimension measurement. The fractal dimension (FD) [40], as described by utilizing Lyapunov exponents, is given by the following equation:(33)dL=j+∑i=1jhihj,where Lyapunov exponents h1,h2,…,hn satisfy ∑i=1jhi≥0 and ∑i=1j+1hi<0.For system (4), the FD takes the form(34)dL=1+h1h2,h1>0>h2,h1+h2<0.The FDs of system (4) are computed numerically with parameter values as in Table 1 and are plotted in Figures 2, 4, and 6. The strange attractors of system (4) (see Figure 5) and its associated FD (see Figure 4(d)) demonstrate that the rise in parameter δ values results in chaotic dynamics for the predator-prey system (4).


### 4.2. Power Spectral Density (PSD) and Recurrence Plot (RP) of the System (4)

To quantify the irregular behavior of system (4), the numerical diagnostics, power spectral density, and recurrence plot [41] are also generated. PSD, which depicts how a signal’s strength varies with frequency. A recurrence plot (RP) in chaos theory is a graph that displays the iterations during which phase space trajectories reach nearly the same region of the phase space for a given instant in time.Example 5.
PSD and RP of system (4)
For the values of the parameter in case (ii), the PSD and RP of the prey population are displayed in Figure9. This figure’s irregular wide peaks are an indication of chaos and randomness. The chaotic character of the system is demonstrated by the random points on the time plane, which guarantee that the almost identical iterative values occur without any rhythm.Figure 9
PSD and Recurrence plot of system (4).: (a) PSD and (b) Recurrence plot.
(a)(b)


## 4.1. Fractal Dimension of System (4)

Strange attractors in a system are characterized by their fractal dimension measurement. The fractal dimension (FD) [40], as described by utilizing Lyapunov exponents, is given by the following equation:(33)dL=j+∑i=1jhihj,where Lyapunov exponents h1,h2,…,hn satisfy ∑i=1jhi≥0 and ∑i=1j+1hi<0.For system (4), the FD takes the form(34)dL=1+h1h2,h1>0>h2,h1+h2<0.The FDs of system (4) are computed numerically with parameter values as in Table 1 and are plotted in Figures 2, 4, and 6. The strange attractors of system (4) (see Figure 5) and its associated FD (see Figure 4(d)) demonstrate that the rise in parameter δ values results in chaotic dynamics for the predator-prey system (4).


## 4.2. Power Spectral Density (PSD) and Recurrence Plot (RP) of the System (4)

To quantify the irregular behavior of system (4), the numerical diagnostics, power spectral density, and recurrence plot [41] are also generated. PSD, which depicts how a signal’s strength varies with frequency. A recurrence plot (RP) in chaos theory is a graph that displays the iterations during which phase space trajectories reach nearly the same region of the phase space for a given instant in time.Example 5.
PSD and RP of system (4)
For the values of the parameter in case (ii), the PSD and RP of the prey population are displayed in Figure9. This figure’s irregular wide peaks are an indication of chaos and randomness. The chaotic character of the system is demonstrated by the random points on the time plane, which guarantee that the almost identical iterative values occur without any rhythm.Figure 9
PSD and Recurrence plot of system (4).: (a) PSD and (b) Recurrence plot.
(a)(b)


## 5. Controlling Chaos

In order to regulate the chaos that exists in system (4), the technique of state feedback control [39] is applied. System (4) may be made into a controlled version by introducing a controlled force un.(35)xn+1=xn+δrxnlnKxn−α1−e−axnyn+unyn+1=yn+δβ1−e−axnyn−dyn,un=−k1xn−x∗−k2yn−y∗,where k1andk2 stand for the feedback gains and x∗,y∗ stands the fixed point of system (4). The controlled system’s Jacobian matrix Jc is provided by the following equation:(36)Jcx∗,y∗=j11−k1j12−k2j21j22,where jpq,p,q=1,2 given in (6), are evaluated at x∗,y∗. The characteristic equation of (36) is(37)λ2−trJcλ+detJc=0,where trJc=j11+j22−k1 and detJc=j22j11−k1−j21j12−k2. Let λ1andλ2 be the solutions of (37).Then,(38)λ1+λ2=j11+j22−k1,and(39)λ1λ2=j22j11−k1−j21j12−k2.Solving the equationsλ1=±1andλ1λ2=1 yields the marginal stability lines. These circumstances demonstrate that λ1,2<1. If we assume that λ1λ2=1, then from (39) we obtain(40)l1:j22k1−j21k2=j11j22−j12j21−1.Assume thatλ1=1, then from (38) and (39) we get(41)l2:1−j22k1+j21k2=j11+j22−1−j11j22+j12a21.Next, forλ1=−1, equations (38) and (39) yield(42)l3:1+j22k1−j21k2=j11+j22+1+j11j22−j12j21.Therefore, the triangular area enclosed by the linesl1,l2,andl3 (see Figure 10(a)) in the k1,k2-plane maintains eigenvalues satisfy λ1,2<1.Figure 10
Chaos control in system (15): (a) stability region in k1,k2-plane and (b) controlled phase diagram.
(a)(b)We fix parameterδ=0.785 and rest as in Example 2 and let x0,y0=0.2462,0.6793. We choose the feedback gains as k1=1.0andk2=1.7 from the stable region (triangular area) in k1,k2-plane plotted in Figure 10(a). The chaotic trajectory is then demonstrated statistically to be stabilized at the fixed point 0.2562,0.6893, see Figure 10(b).


## 6. Discussion

By using a functional response of the Ivlev-type and a prey growth rate of the Gompertz type, we analyze the dynamics of a discrete predator-prey system. With the aid of the center manifold theory, we deduce the existence conditions and directions of flip and Neimark–Sacker bifurcations close to the interior fixed point of system (4) when the bifurcation parameter exceeds a certain threshold value. By presenting numerical simulations that demonstrate unpredictable behaviors of the system via periodic orbits with periods of 2, 4, 8, and 16 through flip bifurcation and periodic orbits with periods of 26, 27, 52, 53, and 80 through NS bifurcations, as well as quasiperiodic and trigger routes to chaos, we corroborate our analytical conclusions. These show that when the dynamics of the prey are chaotic, the system is unstable and the predator either becomes dies out or reaches a stable fixed point. The existence of instability in the system is further supported by numerical estimations of the maximal Lyapunov exponents and the fractal dimension. Additionally, 2D parametric basins of attraction are shown to demonstrate how system dynamics change qualitatively as parameter values increase. Finally, the feedback control approach has been used to govern the system’s chaotic trajectory. The system’s many parameter bifurcations remain a difficult challenge to solve, nevertheless. As we predict, more studies on this subject should produce more analytical findings.



---

*Source: 4555469-2023-02-13.xml*