Discrete-Time Predator-Prey Interaction with Selective Harvesting and Predator Self-Limitation

Institute of Computing Science and Technology, Guangzhou University, Guangzhou 510006, China Department of Mathematics, University of Poonch, Rawalakot, Azad Kashmir, Pakistan Department of Mathematics and Statistics, University of Lahore, Lahore 54000, Pakistan NUST Institute of Civil Engineering, School of Civil and Environmental Engineering, National University of Sciences and Technology, Sector H-12, Islamabad, Pakistan


Introduction
In recent years, the prey-predator system as the fundamental structure in dynamics has got much attention from researchers [1,2]. In these studies, there are numerous components that affect dynamical properties of the prey-predator systems.
rough mathematical modeling, study of a dynamical behavior among prey and predator has a lot of importance to get general perspicacity during the last few years [3]. e preypredator interaction with selective harvesting has been investigated by many authors in differential and difference equations.
In addition, the harvesting has a significant impact on the dynamics of the system [4,5]. In the harvesting system, the main goal is to establish how much we can harvest without altering dangerously the harvested population. Peng et al. [6] examined local stability and Hopf bifurcation in a delayed predator-prey system with a prey refuge and selective harvesting model. Chakraborty et al. [7] investigated global stability, bifurcation, and optimal control for a stage structured prey-predator fishery model with harvesting. Xiao et al. [8] discussed coexistence, bifurcation, and limit cycle in a ratio-dependent prey-predator interaction with predator harvesting model. Similarly, for some other discussions associated with qualitative behavior of selective harvesting models, we refer the interested reader to [9][10][11][12][13][14][15][16].
Taking into account the selective harvesting and predator self-limitation [17], the following predator-prey model was proposed: where x(t), y(t), and z(t) represent population densities for prey population, predator population, and harvesting effort at time t, respectively. Furthermore, it is supposed that prey population grows logistically with intrinsic growth rate r and k is used for carrying capacity of environment. e availability of prey for predation is taken into account by implementing Holling type II functional response (αx/(β + x)) in which α denotes maximum harvesting rate and α is the half-saturation constant. Moreover, a represents the catchability coefficient, b is the intrinsic growth rate of predator, and (1/c) is quantity of prey required to support one consumer at equilibrium when y equals cx. On the other hand, δ is a parameter for stiffness which is used to measure the intensity of reaction between the perceived rent and the effort, m and d are the price for unit catch and the cost of unit effort, respectively, and η represents tax on per unit biomass of the landed fish.
In order to investigate more richer and chaotic behavior of (1), we implement piecewise constant argument method [18] to system (1) as follows: In this paper, qualitative behavior of model (2) is discussed. Remaining discourse of this paper is summarized as follows. Existence of equilibria for system (2) and local stability analysis for these equilibria are presented in Section 2. An explicit criterion related to emergence of Neimark-Sacker bifurcation in system (2) is discussed in Section 3. In Section 4, an explicit criterion related to occurrence of period-doubling bifurcation is implemented. We discuss chaos control techniques for the discrete-time model under the influence of bifurcations in Section 5. Finally, in Section 6, numerical simulations are provided to illustrate all theoretical discussions.

Existence of Equilibria and Stability
e fixed points of system (2) satisfy the following algebraic system: Simple computation yields the following fixed points of system (2): , and an interior fixed point is given by where Furthermore, existence conditions for these fixed points of system (2) are summarized in Table 1.
Lemma 1 (see [19]). Assume that A 2 , A 1 , and A 0 are real constants. en, the necessary and sufficient conditions for all roots of the equation to lie inside the unit disk are First, we see the dynamics of fixed point E 1 and one can easily observe that this equilibrium point is unstable. Indeed, the variational matrix of system (2) at E 1 is given as follows: Journal of Mathematics Here, λ 1 � 1 − r, λ 2 � e b , and λ 3 � e − δ(d+ak(η− m)) are eigenvalues of V(E 1 ) and λ 2 always lies outside the unit disk. Furthermore, E 1 is a saddle point if 0 < r < 2 or d + ak(η − m) > 0. On the other hand, E 1 is a source if d + ak(η − m) < 0 and r > 2. Assume that η < m and (d/k(m − η)) < a; then, E 2 exists and variational matrix at this fixed point is given as follows: Now, characteristic polynomial for V(E 2 ) is given by en, obviously λ � e b is an eigenvalue of V(E 2 ), so E 2 is also unstable. On the other hand, E 2 is a saddle point if (1 − (d/(akm − aηk)))(dr(δ − ((dδ + 2)/(akm − aηk))) + 4) < 0 or |2 − (dr/(akm − aηk))| < 2 + dr(δ − ((dδ + 1)/(akm− aηk))) < 2, and it is a source if and only if the following conditions are satisfied: For fixed point E 3 , the variational matrix of system (2) evaluated at E 3 is given by where a 11 ≔ k 2 2αc(αc + 1) − αcr + r 2 + βkr(3αc + 2r) + β 2 r 2 2αck 2 − (2αck + kr + βr) Moreover, characteristic polynomial of V(E 3 ) is computed as follows: Exists for all of the biologically feasible parameter values e following lemma gives dynamics of system (2) with respect to its fixed point E 3 .

Lemma 2.
e following statements hold true for equilibrium E 3 of system (2): (2) is a sink if and only if

(iii) Fixed point E 3 of system (2) is a source if and only if
In order to illustrate Lemma 2 for particular parametric values, the topological classification of E 3 is depicted in Figure 2 in rb-plane. Moreover, in Figure 2, the blue region represents sink, the red region represents source, and three remaining white regions represent saddle. On the other hand, other parameters are taken as k � 1.13, α � 1.03, Assume that η < m and ar > (aαc d/(d + aβ(m − η))) +(dr/k(m − η)); then, variational matrix V(E * ) of (2) evaluated at unique positive equilibrium point E * is computed as follows: Moreover, the characteristic equation of V(E * ) is given by where Lemma 3. Suppose that η < m and ar > (aαc d/ (d + aβ(m − η))) + (dr/k(m − η)); then, unique positive equilibrium point E * of (2) is a sink if and only if where B 2 , B 1 , and B 0 are listed in (21).
Stability region of positive equilibrium point E * is depicted in Figure 3.

Neimark-Sacker Bifurcation
is section is related to emergence of Hopf bifurcation about positive fixed E * of system (2). For this, we apply a direct criterion for occurrence of Hopf bifurcation without finding the eigenvalues of the system under consideration (cf. [20][21][22][23][24]). Generally, the explicit criterion for occurrence of Hopf bifurcation in n-dimensional discrete system is given as follows.

Journal of Mathematics
Moreover, it is supposed that the following conditions hold true: . ., 2 (or 1), when n is odd (or even, respectively).

(C2)
Transversality criterion: Choosing n � 3, the following result gives conditions under which (2) undergoes Neimark-Sacker bifurcation when b is taken as bifurcation parameter.

Lemma 5.
e positive fixed point E * of system (2) undergoes Neimark-Sacker bifurcation at b � b 0 if the following conditions hold true: where B 2 , B 1 , and B 0 are given in (21), and b 0 can be any possible real root of Proof. Choosing n � 3 and b as bifurcation parameter, from Lemma 4 and characteristic equation (20), one has

Period-Doubling Bifurcation
is section is related to emergence of flip bifurcation about positive fixed E * of system (2). For this, we apply a direct criterion for occurrence of flip bifurcation without finding the eigenvalues of the system under consideration. Generally, the explicit criterion for occurrence of flip bifurcation in n-dimensional discrete system is given as follows.

then, period-doubling bifurcation occurs at critical value μ 0 .
Choosing n � 3 and μ � b as bifurcation parameter, the following result provides us parametric conditions under which system (2) undergoes period-doubling bifurcation whenever b is taken as bifurcation parameter.

Lemma 7.
e unique positive fixed point E * of system (2) undergoes flip bifurcation at b � b 0 if the following conditions hold: where B 2 , B 1 , and B 0 are given in (21)

Chaos Control
is section is dedicated to implementation of chaos control methods to system (2). First we apply a chaos control method based on state feedback control methodology and this method is known as OGY method (that is, Ott-Grebogi-Yorke method) [27]. e OGY control method is one of the pioneer methods for controlling chaos in discretetime systems. On the other hand, it has some deficiencies due to some limitations [28]. In particular, the OGY method may not be applicable to a discrete-time model which is obtained through implementation of Euler approximation [29,30]. For some other methods related to chaos control in discrete-time models, we refer to [29][30][31][32][33][34][35][36][37][38][39][40].
For the application of the OGY method to system (2), the corresponding control system is written as follows: Moreover, b is some nominal value of b which is located in bifurcating or chaotic region and k 1 , k 2 , and k 3 are control parameters. Variational matrix of (31) evaluated at (x * , y * , z * ) is given by is variational matrix is similar to Jacobian matrix of original system (2) at its positive fixed point E * . Since this variational matrix is independent of control parameters k 1 , k 2 , and k 3 , system (31) is not controllable via the OGY method.
Next, we introduce another chaos control method based on random state feedback control strategy as follows: Journal of Mathematics where x * � (d/a(m − η)), y * � (c d/a(m − η)) and z * � ((a(r − (αc d/(aβ(m− η) + d))) − (dr/k(m − η)))/a 2 ). Moreover, s 1 , s 2 , and s 3 are control parameters. We compute variational matrix of (33) about its positive fixed point as follows: Keeping in view this variational matrix of control system (33), its characteristic equation is calculated as follows: where en, it is easy to see that system (33) is controllable as long as its positive fixed point is a sink. Consequently, one has the following result.

Lemma 8. e unique fixed point E * of system (33) is a sink if and only if
At the end of this section, we present a simple and effective method for larger classes of discrete-time systems.
is chaos control method is based on state feedback control and parameter perturbation. We call this method as the hybrid chaos control method [41]. For application of the hybrid control method to system (2), the corresponding control system can be written as follows: where 0 < ϑ < 1 is control parameter. Meanwhile, the Jacobian matrix for system (38) at its positive fixed point is given by So, the characteristic polynomial of J H (E * ) is given by where , en, it is easy to see that system (38) is controllable as long as its positive fixed point is a sink. Consequently, one has the following result.

Lemma 9.
e unique fixed point E * of system (38) is a sink if and only if

Numerical Simulation
First of all, we choose r � 4.  (2) is given by e roots of this characteristic equation are λ 1 � − 1, λ 2 � − 0.634182 and λ 3 � 0.212512. erefore, existence of period-doubling bifurcation with the help of criterion of eigenvalues is justified. Next, all axioms of Lemma 7 at b � b 0 are satisfied as follows: e bifurcation diagrams and maximum Lyapunov exponents (MLE) are depicted in Figure 4.
Furthermore, in order to observe the controllability of system (33), we choose r � 4.8, k � 1.5, α � 2.5, β � 1.7, For these parametric values, the controllable region (green) for system (33) is depicted in Figure 5. Secondly, for similar parametric values, one can see the effectiveness of hybrid control strategy (38). en, some simple calculation yields that system (38)

Concluding Remarks
ere are abundant studies dedicated to the effects of harvesting on population growth. In particular, various classes of prey-predator interaction are studied under the influence of harvesting. Most of these investigations are dedicated to predator-prey interactions with overlapping generations. In order to see effect of selective harvesting on the dynamics of prey-predator interaction with nonoverlapping generation, a discrete-time predatorprey model is studied. Our investigation reals that proposed model undergoes period-doubling and Neimark-Sacker bifurcations around its interior fixed point whenever intrinsic growth rate of predator is taken as bifurcation parameter. Consequently, the discrete-time model has rich dynamical behavior compared to its continuous counterpart. Furthermore, chaotic and bifurcating behavior of the model is controlled with implementation of various chaos control strategies [42][43][44][45].

Data Availability
e data used to support the findings of the study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest. : Controllable region (green) for system (33).