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

Authors: Zhihua Chen; Qamar Din; Muhammad Rafaqat; Umer Saeed; Muhammad Bilal Ajaz
Journal: Journal of Mathematics
 (2020)

---

## Abstract

Selective harvesting plays an important role on the dynamics of predator-prey interaction. On the other hand, the effect of predator self-limitation contributes remarkably to the stabilization of exploitative interactions. Keeping in view the selective harvesting and predator self-limitation, a discrete-time predator-prey model is discussed. Existence of fixed points and their local dynamics is explored for the proposed discrete-time model. Explicit principles of Neimark–Sacker bifurcation and period-doubling bifurcation are used for discussion related to bifurcation analysis in the discrete-time predator-prey system. The control of chaotic behavior is discussed with the help of methods related to state feedback control and parameter perturbation. At the end, some numerical examples are presented for verification and illustration of theoretical findings.


---

## Body

## 1. 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. Through 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]. The prey-predator 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–16].Taking into account the selective harvesting and predator self-limitation [17], the following predator-prey model was proposed:(1)dxdt=rx1−xk−αxyβ+x−axz,dydt=by1−ycx,dzdt=δzaxm−η−d,where xt, yt, and zt 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. The 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:(2)xn+1=xnexpr1−xnk−αynβ+xn−azn,yn+1=ynexpb1−yncxn,zn+1=znexpδaxnm−η−d.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.


## 2. Existence of Equilibria and Stability

The fixed points of system (2) satisfy the following algebraic system:(3)x=xexpr1−xk−αyβ+x−az,y=yexpb1−ycx,z=zexpδaxm−η−d.Simple computation yields the following fixed points of system (2):(4)E1=k,0,0,E2=dam−η,0,ra2a−dkm−η,E3=ξ+ξ2+4βkr22r,cξ+ξ2+4βkr22r,0,and an interior fixed point is given by(5)E∗=dam−η,cdam−η,ar−αcd/aβm−η+d−dr/km−ηa2,where(6)ξ≔kr−βr−αck.Furthermore, existence conditions for these fixed points of system (2) are summarized in Table 1.Table 1
Existence conditions for fixed points of system (2).
Equilibrium pointExistence conditionsE1Exists for all of the biologically feasible parameter valuesE2η<m and d/km−η<aE3Exists for all of the biologically feasible parameter valuesE∗η<m and ar>aαcd/d+aβm−η+dr/km−ηForr=5.2, k=1.01, α=4.3, β=2.9, a=0.1, b=1.49, c=0.45, d=0.11, δ=2.8, m=3.9, and η=1.5, the nullclines for positive equilibrium are depicted in Figure 1. In Figure 1, x-nullclines, y-nullclines, and z-nullclines are represented by red, blue, and black portions, respectively, and intersection of these nullclines is 0.458333,0.20625,25.7618.Figure 1
Existence of unique positive equilibriumE∗.Now, we state the following lemma related to local asymptotic stability of 3-dimensional discrete-time systems.Lemma 1 (see [19]).
Assume thatA2, A1, and A0 are real constants. Then, the necessary and sufficient conditions for all roots of the equation(7)λ3+A2λ2+A1λ+A0=0,to lie inside the unit disk are(8)A2+A0<1+A1,A2−3A0<3−A1,A02+A1−A0A2<1.First, we see the dynamics of fixed pointE1 and one can easily observe that this equilibrium point is unstable. Indeed, the variational matrix of system (2) at E1 is given as follows:(9)VE1=1−r−kαk+β−ak0eb000e−δd+akη−m.Here,λ1=1−r, λ2=eb, and λ3=e−δd+akη−m are eigenvalues of VE1 and λ2 always lies outside the unit disk. Furthermore, E1 is a saddle point if 0<r<2 or d+akη−m>0. On the other hand, E1 is a source if d+akη−m<0 and r>2. Assume that η<m and d/km−η<a; then, E2 exists and variational matrix at this fixed point is given as follows:(10)VE2=1−drakm−akη−dαd+aβm−ηdη−m0eb0rδm−dak−η01.Now, characteristic polynomial forVE2 is given by(11)λ−ebλ2−2−drakm−aηkλ+1+drδ−dδ+1akm−aηk=0.Then, obviouslyλ=eb is an eigenvalue of VE2, so E2 is also unstable. On the other hand, E2 is a saddle point if 1−d/akm−aηkdrδ−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:(12)2−drakm−aηk<2+drδ−dδ+1akm−aηk,1+drδ−dδ+1akm−aηk>1.For fixed pointE3, the variational matrix of system (2) evaluated at E3 is given by(13)VE3=a11a12a13bc1−b000a33,where(14)a11≔k22αcαc+1−αcr+r2+βkr3αc+2r+β2r22αck2−2αck+kr+βrαck−kr+βr2+4βkr22αck2,a12≔kr−βr−αck+αck−kr+βr2+4βkr22ck,a13≔−akr−βr−αck+αck−kr+βr2+4βkr22r,a33≔expδam−ηkr−βr−αck+αck−kr+βr2+4βkr22r−d.Moreover, characteristic polynomial ofVE3 is computed as follows:(15)λ−a33λ2−1−b+a11λ+1−ba11−bca12=0.The following lemma gives dynamics of system (2) with respect to its fixed point E3.Lemma 2.
The following statements hold true for equilibriumE3 of system (2):(i)
Fixed pointE3 of system (2) is a sink if and only if(16)a33<1,1−b+a11<1+1−ba11−bca12<2.(ii)
Fixed pointE3 of system (2) is a saddle point if and only if(17)b1−a11−ca122−b1+a11−bca12<0,a33<1,orb1−a11−ca122−b1+a11−bca12<0,a33>1,ora33>1,and1−b+a11<1+1−ba11−bca12<2,ora33<1,1−b+a11<1+1−ba11−bca12,1−ba11−bca12>1.(iii)
Fixed pointE3 of system (2) is a source if and only if(18)a33>1,1−b+a11<1+1−ba11−bca12,1−ba11−bca12>1.In order to illustrate Lemma2 for particular parametric values, the topological classification of E3 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, β=1.1, a=0.86, c=0.8, δ=1.75, m=0.7, and η=0.1.Assume thatη<m and ar>aαcd/d+aβm−η+dr/km−η; then, variational matrix VE∗ of (2) evaluated at unique positive equilibrium point E∗ is computed as follows:(19)VE∗=cαd2d+aβm−η2−rdakm−akη+1−dαd+aβm−ηdη−mbc1−b0δakr−cα−rβm−ηd+a2krβm−η2−d2rakd+aβm−η01.Moreover, the characteristic equation ofVE∗ is given by(20)Pλ=λ3+B2λ2+B1λ+B0=0,where(21)B2=b−αcd2aβm−η+d2+drakm−aηk−3,B1=−αb−2cd2aβm−η+d2−αcddδ−baβm−η+d−dr−b+dδ+2akm−η−2b+δdr+3,B0=αb−1cd2aβm−η+d2−αcdb−1dδ+baβm−η+d−b−1drdδ+1akm−η+b−1dδr+1.Figure 2
Topological classification forE3.Lemma 3.
Suppose thatη<m and ar>aαcd/d+aβm−η+dr/km−η; then, unique positive equilibrium point E∗ of (2) is a sink if and only if(22)B2+B0<1+B1,B2−3B0<3−B1,B02+B1−B0B2<1,where B2, B1, and B0 are listed in (21).Stability region of positive equilibrium pointE∗ is depicted in Figure 3.Figure 3
Stability region ofE∗ in abc-space.


## 3. Neimark–Sacker Bifurcation

This section is related to emergence of Hopf bifurcation about positive fixedE∗ 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–24]). Generally, the explicit criterion for occurrence of Hopf bifurcation in n-dimensional discrete system is given as follows.Lemma 4 (see [25]).
Consider the followingn-dimensional discrete system:(23)Yk+1=FμYk,where μ∈ℝ represents some bifurcation parameter. Furthermore, take into account the following characteristic polynomial for the Jacobian matrix JY∗=mijn×n evaluated about fixed point Y∗∈ℝn of map Fμ:(24)Pμλ=λn+τ1λn−1+⋯+τn−1λ+τn,where τi=τiμ,v, i=1,2,…,n, in which v denotes control parameter or some other parameter to be determined. Take into account a sequence of determinants of the form Di±μ,vi=0n such that D0±μ,v=1, and(25)Di±μ,v=detT1±T2,where(26)T1=1τ1τ2⋯τi−101τ1⋯τi−2001⋯τi−3⋯⋯⋯⋯⋯000⋯1,T2=τn−i+1τn−i+2⋯τn−1τnτn−i+2τn−i+3⋯τn0⋯⋯⋯⋯⋯τn−1τn⋯00τn0⋯00.Moreover, it is supposed that the following conditions hold true:C1 Eigenvalue criterion: Dn−1−μ0,v=0, Dn−1+μ0,v>0, Pμ01>0, −1nPμ0−1>0, Di±μ0,v>0, for i=n−3,n−5,…, 2 (or 1), when n is odd (or even, respectively).C2 Transversality criterion: d/dμDn−1−μ,vμ=μ0≠0.C3 Nonresonance or resonance criterion: cos2π/l≠φ, or resonance condition cos2π/l=φ, where l=3,4,5,…, and φ=1−0.5Pμ01Dn−3−μ0,v/Dn−2+μ0,v; then, Neimark–Sacker bifurcation occurs at critical value μ0.Choosingn=3, the following result gives conditions under which (2) undergoes Neimark–Sacker bifurcation when b is taken as bifurcation parameter.Lemma 5.
The positive fixed pointE∗ of system (2) undergoes Neimark–Sacker bifurcation at b=b0 if the following conditions hold true:(27)1−B1+B0B2−B0=0,1+B1−B0B2+B0>0,1+B2+B1+B0>0,1−B2+B1−B0>0,ddb1−B1+B0B2−B0b=b0≠0,cos2πl≠1−1+B2+B1+B021+B0,l=3,4,5,…,where B2, B1, and B0 are given in (21), and b0 can be any possible real root of 1−B1b+B0bB2b−B0b=0.Proof.
Choosingn=3 and b as bifurcation parameter, from Lemma 4 and characteristic equation (20), one has(28)D2−b=1−B1+B0B2−B0=0,D2+b=1+B1−B0B2+B0>0,Pb1=1+B2+B1+B0>0,−13Pb−1=1−B2+B1−B0>0,ddbD2−bb=b0=ddb1−B1+B0B2−B0b=b0≠0,1−0.5Pb1D0−bD1+b=1−1+B2+B1+B021+B0.


## 4. Period-Doubling Bifurcation

This section is related to emergence of flip bifurcation about positive fixedE∗ 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.Lemma 6 (see [26]).
Consider the followingn-dimensional discrete system:(29)Yk+1=FμYk,where Yk∈ℝn and μ∈ℝ represents some bifurcation parameter. Moreover, assume that conditions (24)–(26) of Lemma 4 are also satisfied. Moreover, it is supposed that the following conditions hold true:H1 Eigenvalue criterion: Pμ0−1=0, Dn−1±μ0,v>0, Pμ01>0, Di±μ0,v>0, i=n−2,n−4,…, 1 (or 2), when n is even (or odd, respectively).H2 Transversality criterion: ∑i=1n−1n−iτi′/∑i=1n−1n−in−i+1τi−1≠0, where τi′ denotes derivative of τμ at μ=μ0; then, period-doubling bifurcation occurs at critical value μ0.Choosingn=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.
The unique positive fixed pointE∗ of system (2) undergoes flip bifurcation at b=b0 if the following conditions hold:(30)1−B1+B0B2−B0>0,1+B1−B0B2+B0>0,1+B2+B1+B0>0,1−B2+B1−B0=0,1±B0>0,B2′−B1′+B0′3−2B2+B1≠0,where B2, B1, and B0 are given in (21), Bi′ is derivative of Bib at b=b0, and b0 can be any possible real root of algebraic equation 1−B2b+B1b−B0b=0.


## 5. Chaos Control

This 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]. The OGY control method is one of the pioneer methods for controlling chaos in discrete-time 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–40].For the application of the OGY method to system (2), the corresponding control system is written as follows:(31)xn+1=xnexpr1−xnk−αynβ+xn−azn,yn+1=ynexpb˜−k1xn−x∗−k2yn−y∗−k3zn−z∗1−yncxn,zn+1=znexpδaxnm−η−d,where x∗=d/am−η, y∗=cd/am−η, and z∗=ar−αcd/aβm−η+d−dr/km−η/a2. Moreover, b˜ is some nominal value of b which is located in bifurcating or chaotic region and k1, k2, and k3 are control parameters. Variational matrix of (31) evaluated at x∗,y∗,z∗ is given by(32)cαd2d+aβm−η2−rdakm−akη+1−dαd+aβm−ηdη−mbc1−b0δakr−cα−rβm−ηd+a2krβm−η2−d2rakd+aβm−η01.This 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 k1, k2, and k3, 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:(33)xn+1=xnexpr1−xnk−αynβ+xn−azn+s1xn−x∗,yn+1=ynexpb1−yncxn+s2yn−y∗,zn+1=znexpδaxnm−η−d+s3zn−z∗,where x∗=d/am−η, y∗=cd/am−η and z∗=ar−αcd/aβm−η+d−dr/km−η/a2. Moreover, s1, s2, and s3 are control parameters. We compute variational matrix of (33) about its positive fixed point as follows:(34)cαd2d+aβm−η2−rdakm−akη+1+s1−dαd+aβm−ηdη−mbc1−b+s20δakr−cα−rβm−ηd+a2krβm−η2−d2rakd+aβm−η01+s3.Keeping in view this variational matrix of control system (33), its characteristic equation is calculated as follows:(35)λ3+C2λ2+C1λ+C0=0,where(36)C2=b−m11−s2−s3−2,C1=−bcm12−m11b−s2−s3−2−s3+1b−s2−1−m13m31,C0=bcm12s3+1−m13m31b−s2−1+m11s3+1b−s2−1,m11=αcd2aβm−η+d2−drakm−aηk+s1+1,m12=−αdaβm−η+d,m13=dη−m,m31=δa2βkrm−η2+adm−ηkr−αc−βr−d2rakaβm−η+d.Then, 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.
The unique fixed pointE∗ of system (33) is a sink if and only if(37)C2+C0<1+C1,C2−3C0<3−C1,C02+C1−C0C2<1.At the end of this section, we present a simple and effective method for larger classes of discrete-time systems. This 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:(38)xn+1=ϑxnexpr1−xnk−αynβ+xn−azn+1−ϑxn,yn+1=ϑynexpb1−yncxn+1−ϑyn,zn+1=ϑznexpδaxnm−η−d+1−ϑzn,where 0<ϑ<1 is control parameter. Meanwhile, the Jacobian matrix for system (38) at its positive fixed point is given by(39)JHE∗=dcdαd+aβm−η2−rakm−akηϑ+1−dαϑd+aβm−η−dϑm−ηbc−b−ϑ+20δakr−cα−rβm−ηd+a2krβm−η2−d2rϑakd+aβm−η01.So, the characteristic polynomial ofJHE∗ is given by(40)μ3+K2μ2+K1μ+K0=0,where(41)K2=−h11+b+ϑ−3,K1=−h12bc−h11b+ϑ−3−h13h31−b−ϑ+2,K0=h12bc+h11b+ϑ−2−h13h31b+ϑ−2,h11=dcdαd+aβm−η2−rakm−akηϑ+1,h12=−dαϑd+aβm−η,h13=−dϑm−η,h31=δakr−cα−rβm−ηd+a2krβm−η2−d2rϑakd+aβm−η.Then, 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.
The unique fixed pointE∗ of system (38) is a sink if and only if(42)K2+K0<1+K1,K2−3K0<3−K1,K02+K1−K0K2<1.


## 6. Numerical Simulation

First of all, we chooser=4.8, k=1.5, α=2.5, β=1.7, a=0.1, c=0.4, d=0.1, δ=3.2, m=3.5, η=1.6, and bifurcation parameter b∈1.5,3.4. With the variation in bifurcation parameter in interval 1.5,3.4, system (2) undergoes period-doubling bifurcation at b≡b0=2.79335. Moreover, at r=4.8, k=1.5, α=2.5, β=1.7, a=0.1, c=0.4, d=0.1, δ=3.2, m=3.5, η=1.6, and b=2.79335, the unique positive fixed point of (2) is E∗=0.526316,0.210526,28.7938, and characteristic polynomial about this steady state for variational matrix of system (2) is given by(43)λ3+1.42167λ2+0.286899λ−0.134771=0.The roots of this characteristic equation areλ1=−1, λ2=−0.634182 and λ3=0.212512. Therefore, existence of period-doubling bifurcation with the help of criterion of eigenvalues is justified. Next, all axioms of Lemma 7 at b=b0 are satisfied as follows:(44)1−B1+B0B2−B0=0.503337>0,1+B1−B0B2+B0=1.46034>0,1+B2+B1+B0=2.5738>0,1−B2+B1−B0=0,1+B0=0.865229>0,1−B0=1.13477>0,B2′−B1′+B0′3−2B2+B1=2.68723≠0.The bifurcation diagrams and maximum Lyapunov exponents (MLE) are depicted in Figure4.Figure 4
Bifurcation diagrams and MLE for system (2) with r=4.8, k=1.5, α=2.5, β=1.7, a=0.1, c=0.4, d=0.1, δ=3.2, m=3.5, η=1.6, b∈2.5,3.4, and initial conditions x0,y0,z0=0.5,0.2,30. (a) Bifurcation diagram for xn. (b) Bifurcation diagram for yn. (c) Bifurcation diagram for zn. (d) Maximum Lyapunov exponents.
(a)(b)(c)(d)Furthermore, in order to observe the controllability of system (33), we choose r=4.8, k=1.5, α=2.5, β=1.7, a=0.1, c=0.4, d=0.1, δ=3.2, m=3.5, η=1.6, and b=3.2. 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). Then, some simple calculation yields that system (38) is controllable for 0<ϑ<0.755791.Figure 5
Controllable region (green) for system (2).Finally, for verification of Neimark–Sacker bifurcation, we chooser=5.2, k=1.01, α=4.3, β=2.9, a=0.1, c=0.45, d=0.11, δ=2.8, m=3.9, η=1.5, and b∈1,2.5. Then, system (2) encounters Neimark–Sacker bifurcation and b≡b0=1.49673. Consequently, at r=5.2, k=1.01, α=4.3, β=2.9, a=0.1, c=0.45, d=0.11, δ=2.8, m=3.9, η=1.5, and b=1.49673, system (2) has nonhyperbolic equilibrium 0.458333,0.20625,25.7618 since for these parametric values, characteristic polynomial of variation matrix of system (2) is given by(45)λ3+0.820428λ2+0.0258192λ−0.658643=0,with roots λ1,2=−0.739536±0.673117i and λ3=0.658643 such that λ1,2=1. Moreover, at b=b0=1.49673, the conditions of Lemma 5 are fulfilled as follows:(46)D2−b0=1−B1+B0B2−B0=0,D2+b0=1+B1−B0B2+B0=1.13238>0,Pb01=1+B2+B1+B0=1.1876>0,−13Pb0−1=1−B2+B1−B0=0.864034>0,ddbD2−bb=b0=ddb1−B1+B0B2−B0b=b0=−2.94443≠0,1−0.5Pb1D0−bD1+b=1−1+B2+B1+B021+B0=−0.739536.We consider the equationcos2π/l=−0.739536, and one has l=±2.61453. Thus, nonresonance condition is also satisfied. Bifurcation diagrams and MLE are shown in Figure 6. Moreover, in order to check the controllability of system (33), we choose r=5.2, k=1.01, α=4.3, β=2.9, a=0.1, c=0.45, d=0.11, δ=2.8, m=3.9, η=1.5, and b=2.3. For these parametric values, the controllable region (green) for system (33) is depicted in Figure 7. Secondly, for similar parametric values, one can see the effectiveness of hybrid control strategy (38). Then, some simple calculation yields that system (38) is controllable for 0<ϑ<0.812776.Figure 6
Bifurcation diagrams and MLE for system (2) with r=5.2, k=1.01, α=4.3, β=2.9, a=0.1, c=0.45, d=0.11, δ=2.8, m=3.9, η=1.5, b∈1,2.5, and initial conditions x0,y0,z0=0.458,0.206,25.75. (a) Bifurcation diagram for xn. (b) Bifurcation diagram for yn. (c) Bifurcation diagram for zn. (d) Maximum Lyapunov exponents.
(a)(b)(c)(d)Figure 7
Controllable region (green) for system (33).


## 7. Concluding Remarks

There 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 predator-prey 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–45].



---

*Source: 6737098-2020-07-08.xml*