Bifurcation and Chaos of a Discrete Predator-Prey Model with Crowley–Martin Functional Response Incorporating Proportional Prey Refuge

+e paper investigates the dynamical behaviors of a two-species discrete predator-prey system with Crowley–Martin functional response incorporating prey refuge proportional to prey density. +e existence of equilibrium points, stability of three fixed points, period-doubling bifurcation, Neimark–Sacker bifurcation, Marottos chaos, and Control Chaos are analyzed for the discrete-time domain. +e time graphs, phase portraits, and bifurcation diagrams are obtained for different parameters of the model. Numerical simulations and graphics show that the discrete model exhibits rich dynamics, which also present that the system is a chaotic and complex one. +is paper attempts to present a feedback control method which can stabilize chaotic orbits at an unstable equilibrium point.


Introduction
e dynamical behavior of the prey-predator system in ecology has been the subject of study for many researchers and hence they have provided a substantial contribution to the growth of the population models [1][2][3][4][5][6][7][8] in the continuous-time domain. In case of populations that have overlapping generations and the birth processes occuring continuously, the ordinary differential equations are used for modeling of predator-prey interaction . Many species, such as monocarpic plants, and semelparous animals have discrete non overlapping generations, and their births occur in regular breeding seasons. eir interactions are characterized by difference equations or formulated as discretetime mappings. Discrete-time predator-prey models are capable of exhibiting much more complicated dynamics in comparison to the corresponding continuous-time models [31][32][33][34][35][36][37][38][39][40][41][42][43][44][45]. e dynamical study of the prey-predator model with different functional responses is an effective approach for time domain. Upadhyay and Naji [46] presented threespecies food chain model with Crowley-Martin type functional response for the dynamical study. Zhou [47] presented qualitative analysis of the predator-prey model for Crowley-Martin functional responses. Tripathi et al. [48] applied Crowley-Martin functional response on a delayed density dependent predator-prey model. Li et al. [49] presented a predator-prey model with Crowley-Martin functional response for positive periodic solutions.
A certain degree of protection can be offered by nature to a set number of prey populations by providing refuges. e effects of prey refuges on the population dynamics are intricate and complex in nature; however, it can be considered to be constituted by two components for modeling purposes.
e first component of refuge effects, which governs positive growth of prey and negatively that of predators, consists of a reduction in prey mortality owing to decreasing predation success.
e second one may arise as the trade-offs and byproducts arising from the prey's hiding behavior, which could cause advantages or even be detrimental for all the interacting populations. González-Olivares and Ramos-Jiliberto [8] presented prey refuges in a simple prey-predator system. Huang et al. [9] presented stability analysis of a preypredator system incorporating a prey refuge. Ma et al. [12] discussed the effect of prey refuges on a prey-predator system. Chen et al. [13] studied the prey refuge on a Leslie-Gower predator-prey model.
Extensive studies have been carried out by many researchers on refuge concept , and most of these papers consider refuge as a constant number of prey or proportional to prey. ese earlier works have been taken into account that the use of refuges by a fraction of prey, or a constant number prey, is influential in existing a stabilizing effect in the dynamics of the interacting populations. is paper will test the above statement assuming that the quantity of prey in refuges is proportional to the prey. We also conduct an analysis of the dynamic properties of such a system via modification of the well-known Lotka-Volterra predatorprey model with prey self-limitation incorporating Crowley-Martin functional response [46,47]. In this paper, we consider the refuge term proportional to prey density. e rest of the paper is organized as follows. In Section 2, a discrete-time prey-predator model under nonoverlapping generation with refuge is formulated. Analysis of the local stability of the proposed model is dealt with in Section 3. e bifurcation analysis of the proposed model is studied through period-doubling bifurcation and Neimark-Sacker bifurcation in Section 4. Sections 5 and 6 give Marottos chaos and control chaos procedures, respectively. e proposed model is supported with the help of numerical simulations in Section 7. Section 8 makes the conclusion to this paper.

Description of the Proposed Predator-
Prey Model e proposed model belongs to a generalized predator-prey model which incorporates logistic growth having no stage structure for neither the prey nor the predator, which is represented by the following differential equations: where x and y denote the density of prey and predator populations at any time t, respectively, and r, k, c, d, and f are all positive constants and have their biological meanings accordingly. e intrinsic per capita growth rate of prey population is denoted by r, k is the environmental carrying capacity of prey population, c is the maximal per capita consumption rate of predators, d is the efficiency with which predators convert consumed prey into new predators, and f is the per capita death rate of predators. e term φ(x, y) represents the functional response of the predator population.
According to Maynard Smith [4], there is a quantity x r of prey population which incorporates refuges. ereafter, the functional response is modified to incorporate prey refuges in the above model, and the term φ(x, y) becomes the functional response in the form of φ(x − x r , y). Based on a laboratory study of dragonfly, Crowley and Martin present functional response φ(x, y) � x/(1 + αx)(1 + βy), where α > 0 describe the handling time and β > 0 is the magnitude of interference among predators. In the proposed model, we consider refuge x r � bx, and to derive a discrete-time model let (dx/dt) � ((x t+h − x t )/h) and (dy/dt) � ((y t+h − y t )/h), where x t and y t are the densities of the prey and predator populations in discrete time t. Let h ⟶ 1 and f � 1, then the (n + 1)th generation of the prey and predator populations by replacing t by n is obtained as follows: Let us consider (r/k(r + 1)) � 1 and (r + 1) � a, then from (2) we obtain the following proposed discrete-time predator-prey system as follows: where a, b, c, d, α, and β are all positive constants. By 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.

Fixed Points and Stability Analysis of the Prey-Predator System
Fixed points of system (3) are determined by solving the following nonlinear system of equations: By simple calculation, we get three nonnegative fixed points as follows: (i) P 0 � (0, 0), (ii) P 1 � ((a − 1/a), 0), a > 1, and (iii) P 2 � (x 2 , y 2 ), where x 2 is positive root of the equation

Lemma 1. Equation
Now, the local behavior of system (3) is studied for each equilibrium point of the prey-predator system. e stability of system (3) is carried out by computing the Jacobian matrix corresponding to each equilibrium point.
Proposition 3. Interior equilibrium point P 2 (x 2 , y 2 ) is Sink if one of the following conditions holds: Proposition 4. Interior equilibrium point P 2 (x 2 , y 2 ) is Source if one of the following conditions holds: Proposition 5. Interior equilibrium point P 2 (x 2 , y 2 ) is Saddle if one of the following conditions holds: is Nonhyperbolic if one of the following conditions holds: e prey-predator model (3) undergoes period-doubling bifurcation at P 2 (x 2 , y 2 ), when parameters vary in a small neighborhood of PD 1 or PD 2 . Model (3) undergoes Neimark-Sacker bifurcation at P 2 (x 2 , y 2 ), when parameters vary in a small neighborhood of NS 1 or NS 2 .

Bifurcation Analysis
is section carries out an investigation of the conditions for the existence of Neimark-Sacker bifurcation (NSB) and period-doubling bifurcation (PDB) at the positive fixed point P 2 (x 2 , y 2 ) of proposed system (3).
ere is an emergence of different kinds of bifurcations from the fixed Mathematical Problems in Engineering point in dynamical systems, when a particular parameter passes through its critical value. Many dynamical properties of a system can be discussed owing to emergence of NSB and PDB. We discuss NSB and PDB for the positive fixed point P 2 (x 2 , y 2 ) of the prey-predator system (3) taking b as a bifurcation parameter.

Neimark-Sacker Bifurcation.
is section presents NSB of the proposed prey-predator model (3) at P 2 (x 2 , y 2 ) for the parameters are located in the following set: Similar arguments can be applied to the other case: In analyzing the NSB, b is used as the bifurcation parameter. Furthermore, b * (|b * | ⋘ 1) is the perturbation of b, and we consider a perturbation of the model as follows: Let u n � x n − x 2 and v n � y n − y 2 , then the equilibrium point P 2 (x 2 , y 2 ) is transformed into the origin, and further expanding f and g as a Taylor series at (u n , v n ) � (0, 0) to the third order, the proposed model (10) becomes e characteristic equation associated with the linearization of model (11) at (u n , v n ) � (0, 0) is given by λ 2 − Tr(J 1 (b * ))λ + Det(J 1 (b * )) � 0. Roots of the characteristic equation are Now to study the normal form, let c � Im(λ 1,2 ) and δ � Re(λ 1,2 ). We define T � 0 1 c δ , and using the transformation u n v n � T x n y n , model (11) becomes x n+1 � δx n − cy n + f 1 x n , y n , y n+1 � cx n + δy n + g 1 x n , y n , (14) where the functions f 1 and g 1 denote the terms in model (14) in variables (x n , y n ) with the order at least two. In order to undergo Neimark-Sacker Bifurcation, we require the following discriminatory quantity Ω to be nonzero: where From the above analysis we have the following result.

Period-Doubling Bifurcation.
It is clear that one of the eigenvalues of the positive fixed point P 2 (x 2 , y 2 ) is λ 1 � − 1 and the other λ 2 is neither 1 nor − 1 if parameters of the model are located in the following set: Here, we discuss PDB of the proposed model (3) at P 2 (x 2 , y 2 ), when parameters vary in a small neighborhood of PD 1 . Similar arguments can be applied to the other cases.
, and it is obvious that T is invertible. Using the transformation u n v n � T x n y n , then model (19) becomes where the functions f 1 and g 1 denote the terms in model (20) in variables (u n , v n , b * ) with the order at least two. From the center manifold theorem, we know that there exists a center manifold W c (0, 0, 0) of model (20) at (0, 0) in a small neighborhood of b * � 0, which can be approximately described as follows: where Model (20) which is restricted to center manifold W c (0, 0, 0) has the following form: where Mathematical Problems in Engineering For flip bifurcation, we require the two discriminatory quantities ξ 1 and ξ 2 to be nonzero: Finally, from the above analysis, we have the following result.

Proof. For map,
e eigenvalues corresponding with the fixed point P 2 � (x 2 , y 2 ) are given by where erefore, the fixed point P 2 (x 2 , y 2 ) has a pair of complex eigenvalues, and the norm of them exceeds unity if us, we can state the following theorem. Proof. According to the definition of a snap-back repeller, one needs to find one point P * (x 0 , y 0 ) ∈ U 2 such that P * ≠ P 2 , F M (P * ) � P 2 , and |DF M (P * )| ≠ 0, for some positive integer M, where Map F is defined by (3).
To proceed, notice that Mathematical Problems in Engineering Now, a map F 2 has been constructed to map the point P * � (x 0 , y 0 ) to the fixed point P 2 (x 2 , y 2 ) after two iterations if there are solutions different from P 2 for equations (29) and (30). e solutions different from P 2 for equation (30) satisfy the following equation: Substituting x 1 and y 1 into equation (29) and solving x 0 , y 0 , we have By simple calculations, we obtain where e solutions of equations (31) and (32) will farther subject to (x 0 , y 0 ), (x 1 , y 1 ) ≠ (x 2 , y 2 ), (x 0 , y 0 ) ∈ U 2 and |DF 2 (P * )| ≠ 0, then P 2 (x 2 , y 2 ) is a snap-back repeller in U 2 .
us, the following theorem is established.

Chaos Control of the Proposed System
Several methods can be used for obtaining chaos control in discrete-time models. A few of these methods are the state feedback method, pole-placement technique, and hybrid control method. We give a feedback control method to stabilize chaotic orbits at an unstable positive fixed point of the prey-predator system (3), with the following controlled form: 8 Mathematical Problems in Engineering with the following feedback control law as the control force: where q 1 and q 2 are the feedback gain and (x 2 , y 2 ) is a positive fixed point of the model. e Jacobian matrix J for system (35) at (x 2 , y 2 ) is where e corresponding characteristic equation of matrix J is Let λ 1 and λ 2 are the eigenvalues of the above characteristic equation: e lines of marginal stability are determined by solving the equations λ 1 � ± 1 and λ 1 λ 2 � 1.
ese conditions guarantee that the eigenvalues λ 1 and λ 2 have modulus less than 1.

Numerical Simulations
is section presents the qualitative dynamical behaviors of the system by presenting Lyapunov exponent, bifurcation diagram, and phase plane for specific parameter values. We present the numerical simulation to observe rich dynamics of our model, so we consider the hypothetical parameters values of the proposed discrete prey-predator system. e value of the parameters are consider as a � 4.2, b � 0.12, c � 3.0, d � 3.5, α � 0.1, and β � 0.1 with initial population        e orbit diagram of species for a, which is a scaling parameter of the intrinsic per capita growth rate of prey population. Figure 4 shows that the system trajectory evolves    from a fixed point to NSB and finally into a chaotic attractor. ere is one visible periodic windows in the bifurcation diagram for both prey and predator. We also give local magnification of the corresponding bifurcation figure. Figure 10 is the corresponding Lyapunov exponents of Figure 4. Figure 5 shows the orbit diagram of prey and predator population for the refuge parameter (b) with other fixed     parameter values. e system trajectory evolves into a periodic orbit from the chaotic attractor and finally into a fixed point.
ere is only one obvious periodic window in all the bifurcation process for both prey and predator. We also give local magnification of the corresponding bifurcation figure. Figure 11 is the corresponding Lyapunov exponent of Figure 5. Figure 6 shows the orbit diagram of the population for the parameter (c) of per capita consumption rate of predators. e system trajectory evolves from a fixed point to NSB and finally into a chaotic attractor. We also give local magnification of the corresponding bifurcation figure. Figure 12 is the corresponding Lyapunov exponent of Figure 6. Figure 7 shows the orbit diagram of species for the parameter (d) of efficiency with which predators convert the consumed prey into a new predator. e system trajectory evolves from a fixed point to NSB and finally into a chaotic attractor.
ere is one big visible periodic window in the bifurcation diagram for both prey and predator. We also give local magnification of the corresponding bifurcation figure. Figure 13 is the corresponding Lyapunov exponent of Figure 7. Figure 8 shows the orbit diagram of prey and predator population for the parameter α with other fixed parameter values. e system trajectory evolves into a fixed point from the periodic orbit. We also give local magnification of the corresponding bifurcation figure. Figure 14 is the corresponding Lyapunov exponent of Figure 8. Figure 9 shows the orbit diagram of prey and predator population for the parameter β with other fixed parameter values. e system trajectory evolves into a periodic orbit from the chaotic attractor and finally into a fixed point. ere is only one obvious periodic window in all the bifurcation process for both prey and predator. We also give local magnification of the corresponding bifurcation figure. Figure 15 is the corresponding Lyapunov exponent of Figure 9.

Chaos Control.
We observe chaotic behavior of preypredator, as shown in Figure 16(a). In this case fixed point (0.3183, 0.6919) is unstable. In the feedback control method for feedback gain q 1 � 0.5 and q 2 � − 0.5, we observe the fixed point (0.3183, 0.6919) is stable, as shown in Figure 16(b).

Fractal Dimension.
e fractal dimension is defined by using Lyapunov exponents as follows.
Let f be a map on R m . Consider an orbit with Lyapunov exponents h 1 ≥ h 2 ≥ . . . . ≥ h m , and let p denote the largest integer such that  Table 1.

Conclusion
In this paper, we have studied the effects of refuges on a predator-prey interaction, by using the analytical and graphical approach. e refuges are considered as prey refuge proportional to prey density. We evaluate the effects with regard to the local stability of the interior equilibrium point and the long-term dynamics of the interacting populations.
e results show that the effects of refuges can stabilize the interior equilibrium point of the proposed preypredator discrete-time domain model. We have studied bifurcations in a discrete predator-prey model with refuge. Furthermore, it is also shown that the model exhibits various bifurcations of codimension 1, including period-doubling bifurcation, and Neimark-Sacker bifurcation as the values of parameters vary. Numerical simulations and graphical presentation show the rich dynamics in bifurcation and chaotic nature of the system model.

Data Availability
No data were used to support this study.