Dynamics and Control of a Discrete Predator – Prey Model with Prey Refuge: Holling Type I Functional Response

. In this study, we examine the dynamics of a discrete-time predator – prey system with prey refuge. We discuss the stability prerequisite for effective ﬁ xed points. The existence criteria for period-doubling (PD) bifurcation and Neimark – Sacker (N – S) bifurcation are derived from the center manifold theorem and bifurcation theory. Examples of numerical simulations that demonstrate the validity of theoretical analysis, as well as complex dynamical behaviors and biological processes, include bifurcation diagrams, maximal Lyapunov exponents, fractal dimensions (FDs), and phase portraits, respectively. From a biological perspective, this suggests that the system can be stabilized into a locally stable coexistence by the tiny integral step size. However, the system might become unstable because of the large integral step size, resulting in richer and more complex dynamics. It has been discovered that the parameter values have a substantial impact on the dynamic behavior of the discrete prey – predator model. Finally, to control the chaotic trajectories that arise in the system, we employ a feedback control technique.


Introduction
The predator-prey model has both theoretical and real-world uses.A predator-prey model allows for the analysis of potential future events in a dynamic manner.There are several ways that interactions between various species can occur, including competition and predation.The predator-prey relationship is one of the most crucial relationships.Because of its wellknown prevalence and significance, one of the key themes in mathematical ecology is the dynamic interaction of predator and prey.An important population dynamics model that looks at the dynamics of interacting groups [7] is the prey-predator paradigm.The Lotka-Volterra model [1,2], has been employed by population dynamics to comprehend the interaction between ecological species [8][9][10][11][12][13][14].Contrarily, discrete-time models have drawn great interest recently [15-17, 22, 23] because they are better suited to modeling populations with nonoverlapping generations and can produce complex dynamical behaviors than continuous part.The classic predator-prey relationship is given as follows: where the prey and predator population densities are represented by the time-dependent functions xðtÞ and yðtÞ, respectively.All constant is assumed to be positive.The parameter k represents the carrying capacity.The constant δ represents the predator mortality rate.the functional response denoted by ΘðxÞ, whereas ΩðxÞ represents the uptake functions.Each population in an ecological system employs a unique strategy, such as refuging, clustering, etc., to locate food sources and defend itself.Numerous ecological features and elements are employed to build more accurate mathematical models.In population dynamics, the functional response or the ratio of a predator's prey consumption to the density of prey per unit of time must be considered in every prey-predator contact [29,30].For the majority of arthropod predators, the functional response of the Holling type I, II [18], III, and IV is widely used.Later, the Lotka-Volterra model was investigated by Rosenzweig and MacArthur [24] using a logistic growth rate for the prey and a Holling type II functional response to account for the saturation of the predator.A limit cycle emerges when the stable fixed point experiences the Hopf bifurcation, which makes the Rosenzweig and MacArthur model one of the fundamental models since prey-predator cohabitation is not restricted to a stable fixed point.Rosenzweig and MacArthur model's discretetime variant was examined by Hadeler and Gerstmann [25].A straightforward discrete-time prey-predator model with Holling type I incidence was further examined by Danca et al. [10], who showed how chaotic processes may be seen in a straightforward discrete model.In their study of the Rosenzweig and MacArthur prey-predator model with Holling type I, Liu and Xiao [26] provided more evidence that the discrete system displays far richer dynamics than the continuous one.However, a prey refuge offers a more accurate prey-predator model because many prey populations have some sort of available refuge.Maynard [27] demonstrated that while a constant number of refugees of any size changed the neutrally stable behavior into a stable fixed point, the dynamics of the neutrally stable Lotka-Volterra model remained unchanged.In addition, Hassel [28] demonstrated that a big refuge in a model, which displays divergent oscillations in the absence of a refuge, substitutes a stable fixed for the oscillatory behavior.As a result, we note that numerous studies have demonstrated refugia's stabilizing influence on predator-prey relations.Prey refuge has been the subject of certain empirical and theoretical studies, and some of these studies have suggested that prey refuges have a stabilizing effect on predator-prey interactions and can effectively avoid the extinction of prey species [31][32][33][34][35][36][37][38][39][40]47].In a study by Seralan et al. [46], it was explored how the additive type Allee effect in the prey population affected the dynamic difficulties of the Ricker type predator-prey model.The analysis and exploration of the dynamics of a discrete Leslie-Gower predator-prey system with the Allee effect in the predator's population and with fear and Allee effect are observed in a study by Vinoth et al. [48,49], respectively.A detailed exploration of discrete prey-predator model with multistability, torus doubling route to chaos is investigated in a study by Neverova et al. [50], Rajni and Ghosh [51], and Ghosh et al. [52].
According to Pusawidjayanti et al. [19], the following continuous predator-prey model's behavior has been examined as follows: In a study by Khan [20], the author has considered the following discrete-time predator-prey model as follows: where x t and y t , respectively, stand for the prey and predator populations.Parameters β and α are the predator and prey's natural growth rates.By introducing the Allee effect for the prey population in terms of dynamic behavior and Neimark-Sacker (N-S) bifurcation, Kangalgil [21] analyzed the discrete predator-prey model as follows: Allee effect is referred to as x mþx , where m is a positive constant.Utilizing discrete models is another approach to comprehend the challenging issue of prey and predator interaction.In the current work, modification of the model ( 5) is considered by introducing the prey refuge as follows: We focus on the dynamics of model (6) and few of the contributions made by this research include the followings: (1) The suggested discrete-time prey-predator model displays complex dynamics than its continuous equivalent.We looked at the effect of prey refuge on the community of populations in the model.(2) We look for potential fixed points in the stability of the system under study.(3) The analytical result of period-doubling (PD) and N-S bifurcations has been proven.(4) The N-S bifurcation has made the model chaotic, hence the state feedback control procedure has been applied to control it.(5) Some numerical examples for our discrete-time predator-prey model with prey refuge have been supplied in order to confirm the validity of our theoretical results.
The remaining text is organized as follows: The fixed point, topological classes are discussed in Section 2. In Section 3, we analyze the likelihood that the model (6) will exhibit a PD or N-S bifurcation when a particular parametric condition is met.To support the conclusions of our analysis, in Section 4, we quantitatively demonstrate model dynamics 2 Mathematical Problems in Engineering with bifurcation diagrams and phase portraits.In Section 5, we employ state feedback management strategies to control the chaotic model's disorder.In Section 6, a succinct discussion is offered.

!
; where The characteristic equation can be expressed as the following at ζðx * ; y * Þ.
where TrðWζ Þ and DetðWζ Þ are given as follows: Tr The eigenvalues of ( 9) can be derived as One can get the stability requirement of fixed points ζ0 , ζ1 , and ζ2 easily from the Jury's criterion (see [45]) F a ð− 1Þ>0, F a ð0Þ − 1<0, F a ð1Þ>0 and proofs are omitted.Table 2 shows the stability conditions of fixed points, and Figures 1 and 2 show more information.The eigenvalues are real in the right part of the region in ðm a ; rÞ space separated by the dashed line and on the opposite part, the eigenvalues are complex.
Let, m a ¼ M3 ¼ m aPD .Then, the eigenvalues of Wζ are provided as follows: For jλ a2 ð M3 Þ ≠ 1j to be implied as follows: Next, we set Aðm a Þ ¼ Wζ and apply the transformations b x ¼ x − x * ; b y ¼ y − y * .We relocate system (6)'s fixed point to the starting point.Consequently, the system (6) can be expressed as follows: where Y ¼ ðb x; b yÞ T and It is possible to express the system (6) as follows: As symmetric multilinear vector functions on x; y; u 2 So, using straightforward calculation, we arrive at: For ensuring h q1 ; q2 i ¼ 1, where h q1 ; q2 i ¼ q 11 q 21 þ q 12 q 22 , we have to utilize the normalized vector q2 We must examine the sign of s a ðm aPD Þ, the coefficient of the critical standard form [41], to establish the PD bifurcation's direction.
e q 2 ; C e e q 1 ; e q 1 ; e The direction and stability of PD bifurcation can be shown using the following theorem in light of the justification presented above.
3.2.Neimark-Sacker Bifurcation.Next, we take the system (6) at the fixed point ζ2 ðx * ; y * Þ where the parameters ðr; a; The system (6)'s eigenvalues are then complex, λ 1; 2 2 C. Also, Consider the case where q1 ; q2 2 R 2 be two eigenvectors of Aðm aNS Þ and A T ðm aNS Þ for eigenvalue λ a ðm aNS Þ and λa ðm aNS Þ such that: Therefore, by performing simple calculations, we find as follows: To obtain h q1 ; q2 i ¼ 1, where h q1 ; q2 i ¼ q 11 q 21 þ q 12 q 22 , we set the normalized vector q2 By taking into account how m a can fluctuate close to m aNS and for z a 2 C, we can decompose Y 2 R 2 as Y ¼ z a q 1 þ za q1 .z a ¼ hq 2 ; Yi is the exact formulation of z a .Thus, for jm a j close to m aNS , the system (6) switched to the following system as follows: When Taylor expansion is used on the function b h, we get: Symmetric multilinear vector functions can be used to define the Taylor coefficients.b h 20 m aNS ð Þ¼ q 2 ; B e q 1 ; q The first Lyapunov coefficient s 2 ðm aNS Þ sign determines the N-S bifurcation's direction, which is given by the expression as follows: In light of the preceding explanation, the following theorem can be utilized to show the direction and stability of N-S bifurcation.
Theorem 2. Assume that ( 21) is true and that s 2 ðm aNS Þ≠0 is true.If the value of m a fluctuates in a specific area around M4 , system (6) experiences a N-S bifurcation at ζ2 ðx * ; y * Þ.Additionally, if s 2 ðm aNS Þ is negative (resp.positive) and the N-S bifurcation is supercritical (resp.subcritical), a unique invariant closed curve that is attracting (resp.repelling) bifurcates from ζ2 ðx * ; y * Þ as well.

Quantitative Study
In order to support our theoretical findings and demonstrate some novel, intriguing complex dynamical behaviors present in system (6), numerical simulation work has been done to exhibit bifurcation diagrams, phase portraits, Lyapunov exponents, and fractal dimension (FD) of system (6).Before presenting all the scenarios, we discuss the FD first.

Fractal Dimension.
The idea of FD is frequently used in the context of dynamical systems to describe the complexity and self-similarity of the structures inside the system.FD gives an indication of how a set fills space, capturing complex patterns and imperfections that may not be well captured by traditional Euclidean geometry.The FDs measurement, which is defined by Cartwright [42], is used to determine a model's chaotic attractors.
where k is the largest integer number such that ∑ k j¼1 ttt j ≥ 0 and ∑ kþ1 j¼1 ttt j <0 and tt j 's are Lyapunov exponents.Now, the model ( 6)'s fractal dimensions structure is given as follows: Given that the chaotic dynamics of the model (6) (Figure 3) are quantified with the sign of FD (Figure 4(d)), it is inevitable that the dynamics of the model stabilize as the parameter m a increases.
In the following situations, we take into account the bifurcation parameters: Scenario (i) ranging m a between the ranges of 0 and 1 and fixing other parameters as r ¼ 3:55; α a ¼ 3:5; β a ¼ 4:5; δ a ¼ 0:25: The calculated maximum Lyapunov exponent corresponding to Figures 5(a) and 5(b) is shown in Figure 5(c).The chaotic zone has stable fixed points or stable periodic windows as a result of certain Lyapunov exponents being positive and some being negative, as shown in Figure 5(c).The diagrams, as shown in Figures 5(a) and 5(b), demonstrate that the fixed point ζ2 of system ( 6) is unstable up to a scale factor of m a ¼ 0:303922 but becomes stable as the scale factor increases. Figure 6 is arranged with the phase portraits related to Figures 5(a For scenario (ii), Figures 4(a) and 4(b) show the bifurcation diagrams of system (6) in the ðm a − xÞ and ðm a − yÞ planes.We observe that an N-S bifurcation emerges at m a ¼ m aNS ¼ 0:318182 near the fixed point ð0:407407; 0:512169Þ of system (6).At m a ¼ m aNS ; we find λ a1; a2 ¼ The Taylor coefficients are given by b The N-S bifurcation is supercritical as a result, which confirms Theorem 2.
The estimated and displayed maximum Lyapunov exponent for Figures 4(a) and 4(b) is shown in Figure 4(c).The chaotic zone has stable fixed points or stable periodic windows as a result of certain Lyapunov exponents being positive and some being negative, as shown in Figure 4(c).The diagrams, as shown in Figures 4(a) and 4(b), show that the fixed point ζ2 of the system (6) is unstable up to a scale factor of m a ¼ 0:318182 but becomes stable as the scale factor increases.The phase portraits associated with Figures 4(a) and 4(b) for various values of m a are arranged in Figure 3, making it easy to see how a smooth invariant circle separates from the stable fixed point.
For scenario (iii), the bifurcation diagrams of system (6) in the ðβ a − xÞ and ðβ a − yÞ planes are shown in Figures 7(a

Chaos Control
The state feedback method, pole placement methodology, OGY technique, and hybrid control approach are the most frequently used chaos control techniques for discrete-time  Mathematical Problems in Engineering models.The chaos generated by N-S bifurcation, we apply the state feedback approach to control the chaos [43,44] to the system (6) and using m a as a control parameter.We write system (6) as follows: Then, the controlled form of system (32) can be written as follows: where the control force is specified as the feedback gains k 1 and k 2 and ðx * ; y * Þ is the interior fixed point for the system (6).The feedback gains k 1 and k 2 are critical in stabilizing and modifying the behavior of discrete dynamical systems under chaos control (see [44]).The goal of chaos control is to use control techniques to steer a chaotic system toward a desirable state or trajectory.The feedback gains k 1 and k 2 might be compared to regulatory mechanisms that restore the system to a stable state while making sure that critical processes stay within reasonable bounds.Similar to how biological systems are resilient to outside threats, k 1 and k 2 in chaos control can be seen as elements that contribute to the system's capacity to withstand disruptions and return to a desired state.

Conclusions
We examine the dynamics of a discrete predator-prey system utilizing a Holling type I functional response and a prey refuge.We determine the existence conditions and directions of the PD and N-S bifurcations near the interior fixed point of system (6) using the center manifold theory when the bifurcation parameter rises over a predetermined threshold.Notably, our results show that the model exhibits chaotic behavior and that the system becomes unstable when the parameter β a increases, leading to a transition from a stable state to chaotic behavior.Also, the system becomes unstable when the parameter m a increases.These demonstrate that the predator either falls extinct or reaches a stable fixed point when the dynamics of the prey are chaotic.Numerical estimates of the maximal Lyapunov exponents and the FD provide more evidence for the system's instability.Additionally, by varying the two control parameters, the system (6) displays extremely complex nonlinear dynamical behaviors and the chaotic phenomenon may be directly observed in the two-dimensional parameter spaces.Finally, the chaotic trajectory of the system has been controlled using the feedback control strategy.Despite this, resolving the system's numerous parameter, bifurcations remain a challenging task.More research on this topic should lead to more analytical conclusions, as we forecast.
) and 5(b) for different values of m a ; it unmistakably illustrates the procedure by which a smooth invariant circle deviates from the steady fixed point.

22 FIGURE 3 :
FIGURE 3: The phase diagram for altering the input of m a .
) and 7(b) and corresponding MLEs and FDs are shown in Figures 7(c) and 7(d).

Figure 8 (
a) shows the codimension-2 bifurcation diagrams in ðβ a ; m a ; xÞ space.The plot of the maximal Lyapunov exponents for two control parameters is shown in Figure 8(b) through a 2D projection onto the ðβ a ; m a Þ plane.

FIGURE 4 :
FIGURE 4: Visualization of NS bifurcation, MLEs, and FDs of species for changing parameter m a for r = 3.75.

FIGURE 5 : 23 FIGURE 6 :
FIGURE 5: Visualization of NS bifurcation, MLEs, and FDs of species for changing parameter m a for r = 3.55.

FIGURE 7 :
FIGURE 7: Visualization of NS bifurcation, MLEs, and FDs of species for changing parameter β a .

FIGURE 9 :
FIGURE 9: Chaos regulation in system.(a) The stability zone in the ðk 1 ; k 2 Þ plane.(b) The regulated phase diagram.

Table 1
(6)tem(6).It is important to note that, regardless of the magnitude of the predicted eigenvalues at the fixed point ζðx; yÞ, estimated eigenvalues have an effect on the fixed point's local stability.The variational matrix of system (6) is shown as follows:W e ζ x; y