Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase Volume 3, Number 1, March 2022, pp.50-59 https://doi.org/10.5206/mase/14496 DYNAMICS OF A DISCRETE PREDATOR-PREY SYSTEM WITH FEAR EFFECT AND DENSITY DEPENDENT BIRTH RATE OF THE PREY SPECIES DEBASIS MUKHERJEE Abstract. This paper analyses a discrete predator-prey system with fear effect and density dependent birth rate of the prey species. The fixed points of the system are determined and their stability is examined. The criterion for Neimark-Sacker bifurcation and flip bifurcation is developed. The chaotic orbit at an unstable fixed point can be stabilized by applying the state feedback control method. Numerically, we illustrate our analytical findings and observe the complex behaviour of the system that leads to stable state to chaotic one. 1. Introduction In recent years, it is observed that the predator-prey interaction is not only governed by direct killing of prey by the predator, but also the indirect effect such as fear caused by the predator. The fear factor influences the birth rate of prey [1]. Based on the fact of fear effect on the prey’s growth rate, several research works are explored [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. In predator-prey interaction, it is commonly assumed that the birth rate of the prey species is constant. But in a real ecological system, a birth rate of the prey species is dependent on the density of prey. In [11], the authors considered density dependent birth rate of the prey species and discussed the dynamical behaviour of the predator-prey system. Aforesaid studies are mainly confined on continuous predator-prey models with two variables. Although, discrete time models are more appropriate than the continuous system when the populations have nonoverlapping generations and virtually remain constant over a generation. From a biological point of view, a discrete time model is applied to investigate the taxonomic group of organisms and species with the passage of time. There are some biological situations where a discrete time system is applicable. For example, fish populations reproduce at specific time moments or for insect populations, for which nonoverlapping generations are occurring in real ecosystems. Other examples include monocarpic plants and semelparous animals which have nonoverlapping populations and their births take place in usual breeding seasons. Moreover, dynamics of discrete time predator-prey system can exhibit a richer set of patterns than those found in continuous systems [12, 13, 14, 15]. Also discrete time models can exhibit chaotic dynamics [12, 13]. So chaos control becomes an interesting topic of research in discrete dynamical system. We will show chaos control by the state feedback control strategy. In [15], the authors observed flip bifurcation and Neimark-Sacker bifurcation in a discrete predator prey system with Holling type III functional response. Santra et al. [16] analysed a discrete predator- prey model with Crowley-Martin functional response where prey population takes refuge. They showed the effects of refuge on the stability of the system in the discrete domain. Furthermore, they ob- tained period doubling bifurcation and Neimark-Sacker bifurcation. Elettreby et al. [14] addressed the complex behaviour of a discrete prey-predator model considering mixed functional response of Holling Received by the editors 25 November 2021; accepted 15 January 2022; published online 6 February 2022. 2020 Mathematics Subject Classification. 39A28, 39A30, 92D25. Key words and phrases. Predator-prey system; fear effect; density dependent birth rate; bifurcation; chaos control. 50 https://ojs.lib.uwo.ca/mase https://dx.doi.org/https://doi.org/10.5206/mase/14496 DYNAMICS OF A DISCRETE PREDATOR-PREY SYSTEM 51 type I and III. Din [13] investigated the complex nature and chaos prevention in a discrete model of prey-predator interaction and found period doubling and Neimark-Sacker bifurcation for larger range of bifurcation parameter. Seno [17] remarked that the dynamics of discrete prey-predator system is consistent with continuous counterpart. Agiza et al. [12] discussed the dynamics of a discrete-time prey-predator model with Holling type II response function. They derived bifurcation diagrams, phase portraits and Lyapunov exponents for various system parameters. They also computed the fractional dimension of a strange attractor of the model. In this paper, we propose a discrete predator-prey system with fear effect and density dependent birth rate of the prey species. We study the existence and stability of different fixed points. After then, we identify the system parameters that give Neimark-Sacker bifurcation and flip bifurcation. The paper is formatted as follows. In Section 2, we present a discrete model of predator-prey interaction with fear effect and the density dependent birth rate of the prey species. Dynamical analysis of different fixed points is described in Section 3. Chaos control is shown in Section 4. In Section 5, the behaviour of the system is demonstrated when values of parameters are changed. A short discussion is given in Section 6. 2. Discrete model Now, we present the following discrete time predator-prey model:  xn+1 = xn [ r (b + cxn)(1 + kyn) −α−βxn − pyn 1 + hxn ] , yn+1 = yn [ −d + pqxn 1 + hxn ] (2.1) where xn and yn represent population densities of prey and predator respectively, and r,b,c,k,α,β,p,q, h,d are positive constants. Here r denotes the birth rate of the prey which is affected by the fear factor 1/(1 + kyn) where k is the level of fear factor. The birth rate r is modified by the density of prey in the form of Beverton-Holt function [18] as r/(b + cxn) where, b and c are positive parameters. α is the natural mortality rate of prey. β represents the intraspecific competition among the prey species. p,h,q,d represent consumption rate, handling time, q conversion efficiency and the death rate of the predator respectively. 3. Fixed points and their nature In this section, we determine the fixed points and their dynamics. Evidently, System (2.1) has three fixed points E0 = (0, 0),E1 = (x̄, 0) and E ∗ = (x∗,y∗) where x̄ = −{c(1 + α) + βb} + √ {c(1 + α) + βb}2 − 4βc{b(1 + α) −r} 2βc , and   x∗ = 1 + d pq −h(1 + d) , y∗ = −D + √ D2 − 4pk(b + cx∗)(1 + hx∗){(b + cx∗)(α + βx∗ + 1) −r} 2p(b + cx∗) (3.1) where D = (b + cx∗){p + k(1 + hx∗)(α + βx∗ + 1)}. Now E1 is feasible if b(1 + α) < r; E∗ is feasible if pq > h(1 + d) and (b + cx∗)(α + βx∗ + 1) < r. 52 D. MUKHERJEE To determine the nature of the fixed points, we compute the Jacobian matrix at each fixed point. The Jacobian matrix at an arbitrary fixed point (x,y) is given by J(x,y) =   br (b+cx)2(1+ky) −α− 2βx− py (1+hx)2 −x{ rk (b+cx)(1+ky)2 + p 1+hx } pqy (1+hx)2 −d + pqx 1+hx   Proposition 1. The fixed point E0 = (0, 0) of system (2.1) is stable if b(α − 1) < r < b(α + 1) and 0 < d < 1. Proof. The characteristic equation at E0 is∣∣∣∣ rb −α−λ 00 −d−λ ∣∣∣∣ = 0. Thus the eigenvalues are λ1 = r b −α and λ2 = −d. Then the fixed point E0 is locally asymptotically stable if |λi| < 1, i = 1, 2. Now |λ1| = ∣∣r b −α ∣∣ < 1 then α− 1 < r b < α + 1. Also, |λ2| = |−d| < 1 then −1 < d < 1. As d is the death rate of predator, this implies that 0 < d < 1. This completes the proof. Proposition 2. Assume that b(1 + α) < r holds. Then, the fixed point E1 = (x̄, 0) of system (2.1) is stable if the following inequalities are fulfilled: α + 2βx̄− 1 < br (b + cx̄)2 < α + 2βx̄ + 1 and d− 1 < pqx̄ 1 + hx̄ < d + 1. Proof. The characteristic equation at E1 is∣∣∣∣∣∣∣ br (b+cx̄)2 −α− 2βx̄−λ −x̄( rk b+cx̄ + p 1+hx̄ ) 0 −d + pqx̄ 1+hx̄ −λ ∣∣∣∣∣∣∣ = 0. Hence, the eigenvalues are λ1 = br (b + cx̄)2 −α− 2βx̄, λ2 = −d + pqx̄ 1 + hx̄ . The fixed point E1 is locally stable if |λi| < 1, i = 1, 2. Now |λ1| < 1 is equivalent to α + 2βx̄− 1 < br (b + cx̄)2 < α + 2βx̄ + 1 and |λ2| < 1 is equivalent to d− 1 < pqx̄ 1 + hx̄ < d + 1. This completes the proof. We remark that the above stability conditions imply that the predator goes to extinction while prey is there. Proposition 3. Assume that pq > h(1 + d) and (b + cx∗)(α + βx∗ + 1) < r hold. Then, the fixed point E∗ = (x∗,y∗) of system (2.1) is stable if the following inequalities are fulfilled: pqx∗y∗ (1 + hx∗)2 [ rk (b + cx∗)(1 + ky∗)2 + p 1 + hx∗ ] + br (b + cx∗)2(1 + ky∗) < 1 +α+ 2βx∗ + py∗ (1 + hx∗)2 , (3.2) α+2βx∗+ py∗ (1 + hx∗)2 − br (b + cx∗)2(1 + ky∗) < 1+ pqx∗y∗ 2(1 + hx∗)2 [ rk (b + cx∗)(1 + ky∗)2 + p 1 + hx∗ ] . (3.3) DYNAMICS OF A DISCRETE PREDATOR-PREY SYSTEM 53 Proof. The characteristic equation at E∗ is∣∣∣∣∣∣∣ br (b+cx∗)2(1+ky∗) −α− 2βx∗ − py ∗ (1+hx∗)2 −λ −x∗ [ rk (b+cx∗)(1+ky∗)2 + p 1+hx∗ ] pqy∗ (1+hx∗)2 1 −λ ∣∣∣∣∣∣∣ = 0, which is now written in the form λ2 −ηλ + γ = 0, where η = br (b + cx∗)2(1 + ky∗) −α− 2βx∗ − py∗ (1 + hx∗)2 + 1, and γ = br (b + cx∗)2(1 + ky∗) −α− 2βx∗ − py∗ (1 + hx∗)2 + pqx∗y∗ (1 + hx∗)2 [ rk (b + cx∗)(1 + ky∗)2 + p 1 + hx∗ ] . For stability of E∗, we use Jury Criterion which is given by |η| < 1 + γ < 2. This condition has two parts, namely (i) γ < 1, and (ii) −1 − γ < η < 1 + γ. By the formula for γ given above, Part (i) is precisely (3.2). For Part (ii), the left inequality is nothing by (3.3), while the right inequality reduces to 0 < pqx∗y∗ (1 + hx∗)2 [ rk (b + cx∗)(1 + ky∗)2 + p 1 + hx∗ ] which is always true when the interior fixed point exists. From the above analysis, we infer that E∗ is stable under the conditions of the theorem, and the proof is completed. 3.1. Bifurcation around the interior fixed point. In discrete context, the Neimark-Sacker bifurca- tion is the counterpart of the Hopf bifurcation that takes place in continuous systems. It was explored by Neimark [19] and alone by Sacker [20]. Hopf bifurcation generates limit cycles in the phase plane in the continuous models. Alternately, Neimark-Sacker bifurcation produces dynamically invariant cycles. Subsequently, we may get isolated periodic orbits as well as trajectories that cover the invariant circle densely. Biologically, Neimark-Sacker bifurcation implies that all the populations can oscillate around some mean values. Flip bifurcation is another type of bifurcation which is also recognized as period doubling bifurcation and it occurs when a small changes in bifurcation parameters give rise to a new system that bifurcate twice the period as the original system. This bifurcation indicates the loss of stability of a periodic orbit. System (2.1) has at most one unique fixed point E∗, hence the system does not admit fold bifurcation. So we are interested in examining the Neimark-Sacker bifurcation and flip bifurcation in the sequel. Proposition 4. System (2.1) admits Neimarck-Sacker bifurcation at E∗ if the following conditions are satisfied: br (b + cx∗)2(1 + ky∗) + pqx∗y∗ (1 + hx∗)2 [ rk (b + cx∗)(1 + ky∗)2 + p 1 + hx∗ ] = α+ 2βx∗ + py∗ (1 + hx∗)2 + 1, (3.4) and pqx∗y∗ (1 + hx∗)2 [ rk (b + cx∗)(1 + ky∗)2 + p 1 + hx∗ ] < 4. (3.5) Proof.. If the Jacobian matrix J(E∗) has two complex conjugate eigenvalues with modulus 1, Neimark-Sacker bifurcation appears [21]. This requires that det(J(E∗)) = γ = 1 and −2 < tr(J(E∗)) = η < 2. Replacing η and γ (see the proof of Proposition 3), the first condition (γ = 1) is precisely (3.4); and the second condition on η is (3.5). This completes the proof. 54 D. MUKHERJEE Proposition 5. System (2.1) admits a flip bifurcation at E∗ if the following conditions are satisfied: 2 [ 1 + br (b + cx∗)2(1 + ky∗) ] + pqx∗y∗ (1 + hx∗)2 [ rk (b + cx∗)(1 + ky∗)2 + p 1 + hx∗ ] = 2 [ α + 2βx∗ + py∗ (1 + hx∗)2 ] . (3.6) Proof. System (2.1) admits flip bifurcatiopn when a single eigenvalue is −1. Thus the condition for flip bifurcation can be written in the form 1 +η +γ = 0 where γ and η are as in the proof of Proposition 3. This condition is precisely (3.6), and hence completes the proof. 4. Chaos control Chaos control is a technique of stabilization by means of small perturbation which are used to unstable periodic orbits for a given system. Sometimes bifurcation and chaotic behaviour are really undesirable phenomena in discrete dynamical systems, because there may be an extinction of population due to chaos. So controlling chaos is an important issue. There are different methods for controlling chaos, e.g., feedback control strategy, hybrid control technique and pole-placement method. By applying these methods, one can retard or remove the chaotic behaviour due to appearance of bifurcation in the dynamical systems and rebuilt the stability of the system. In this section, we use mainly the state feedback control technique [13] to stabilize a chaotic orbit at an unstable fixed point of system (2.1). Consider the following controlled system related to (2.1):  xn+1 = xn [ r (b + cxn)(1 + kyn) −α−βxn − pyn 1 + hxn ] −u(xn,yn), yn+1 = yn [ −d + pqxn 1 + hxn ] (4.1) where u(xn,yn) = c1(xn − x∗) + c2(yn − y∗) is a feedback controlling force with c1 and c2 being the feedback gains and (x∗,y∗) being the unique fixed point of system (2.1). The Jacobian matrix of system (4.1) evaluated at (x∗,y∗) is given by J(x∗,y∗) = ( m11 − c1 m12 − c2 m21 m22 ) where m11 = br (b + cx∗)2(1 + ky∗) −α− 2βx∗ − py∗ (1 + hx∗)2 , m12 = −x∗( rk (b + cx∗)(1 + ky∗)2 + p 1 + hx∗ ), m21 = pqy∗ (1 + hx∗)2 , m22 = 1. The characteristic equation of the variational matrix J(x∗,y∗) is λ2 − (m11 + m22 − c1)λ + m22(m11 − c1) −m21(m12 − c2) = 0 (4.2) Suppose λ1 and λ2 are the roots of equation (4.2), then we have λ1 + λ2 = m11 + m22 − c1, λ1λ2 = m22(m11 − c1) −m21(m12 − c2). (4.3) The lines of marginal stability can be derived by the equations λ1 = ±1 and λ1λ2 = 1. These restrictions ensure that the eigenvalues λ1 and λ2 have moduli equal to 1. DYNAMICS OF A DISCRETE PREDATOR-PREY SYSTEM 55 First suppose that λ1λ2 = 1. Then from (4.3), we find l1 : c1m22 − c2m21 = m11m22 −m12m21 − 1. Next assume that λ1 = 1. Then from (4.3), we obtain l2 : c1(1 −m22) + c2m21 = m11 + m22 −m11m22 + m12m21. Lastly, assume that λ1 = −1 and from. (4.3), we get l3 : c1(1 + m22) − c2m21 = m11 + m22 + 1 + m11m22 −m12m21. The stable eigenvalues lie within a triangular region bounded by the lines l1, l2 and l3 in the c1-c2 plane. 5. Numerical Simulation In this section, we present some numerical simulation to illustrate the usefulness of the obtained results as well as for giving direction to find desirable bifurcations and chaos of the discrete time system (2.1). In Fig. 1, we select the parameter values r = 4.5,k = 1,b = 1,c = 1,p = 4,q = 1,h = 1,α = 0.1,β = 0.1. We draw the bifurcation diagram with respect to the parameter d in the interval (1.5, 2.8). As d increases, we observe a transition phase from stability to bifurcation within a limit cycle, to a periodic window and ultimately to chaos. In Fig. 2, we select the parameter values r = 4.5,k = 1,b = 1,c = 1,p = 4,q = 1,h = 1,α = 0.1,β = 0.01. Here β is decreased from 0.1 to 0.01 from previously chosen parameters in Fig. 1. We draw the bifurcation diagram with respect to the parameter d in the interval (1.5, 2.8). As d increases, we observe a transition phase from stability to bifurcation within a limit cycle, to a periodic window and ultimately to chaos. Here we observe flip bifurcation. In Fig. 3, we select the parameter values r = 4.5,b = 1,c = 1,p = 1.8,q = 1,h = 1,α = 0.1,β = 0.01. We draw the bifurcation diagram with respect to the parameter k in the interval (2.1, 2.8). As k increases, we observe a transition phase from stability to bifurcation within a limit cycle, to a periodic window and ultimately to chaos. In Fig. 4, we select the parameter values r = 4.5,k = 1,b = 1,c = 1,p = 4,q = 1,h = 1,d = 1,α = 0.1,β = 0.01. and the initial value is (0.1, 0.1). With the above choice of parameters, we find chaotic behaviour of the system. To avoid the chaotic dynamics, feedback gains c1 = 0.3 and c2 = −1.2. are chosen. The chaotic orbit is stabilized at the fixed point (1, 0.306246). 6. Discussion In this article, we have studied the qualitative behaviour of a discrete predator-prey model with fear effect and density dependent birth rate of the prey species. The predator functional response is taken as Holling type II. Prey’s birth rate is assumed to be as Beverton-Holt type function [18]. We have mainly identified the system parameters that affect the dynamics of the system. We have observed two boundary fixed points and a unique interior fixed point. Stability analysis of these fixed points is examined by Jury technique. The criterion for Neimark-Sacker bifurcation and flip bifurcation are used for examining the bifurcations around the positive fixed point. It is identified that the parameter d, the death rate of predator in the system is more relevant for the appearance of flip bifurcation and Neimark-Sacker bifurcation whenever it is varied in some appropriate interval. We have also found Neimark-Sacker bifurcation by varying the parameter k. In investigating bifurcation, we have noted 56 D. MUKHERJEE Figure 1. Bifurcation diagram for prey and predator populations with d for fixed values r = 4.5,k = 1,b = 1,c = 1,p = 4,q = 1,h = 1,α = 0.1 and β = 0.1. that β, the intraspecific competition among the prey species has an important role. It is checked that smaller values of β may result flip bifurcation while larger values for β may result for Neimark-Sacker bifurcation. In [1], the authors studied system (2.1) with b = 1, and c = 0 and remarked that the cost of fear affect the existence of Hopf bifurcation as well as the direction of Hopf bifurcation in the continuous model. But in our discrete time model (2.1), we observed Neimark-Sacker bifurcation and chaotic behaviour of the system varying the fear factor. Recently, Kundu et al. [22] analysed similar type of system with b = 1,c = 0 and h = 0 without obtaining different type of bifurcations and they also observed that the system with fear effect becomes stable from chaotic dynamics by increasing fear factor which is not so in our system (see Fig. 3). The conditions of Proposition 1, shows that when the intrinsic growth rate of prey lies in a certain interval and the death rate of predator remains below a certain threshold value, both the populations go to extinction. If the restrictions of Proposition 2 are satisfied, the predator population goes to extinction while prey population can sustain there. The stable coexistence of all the populations are possible when all the conditions of Proposition 3 hold. The conditions of Proposition 4 suggests that Neimark-Sacker bifurcation is possible for system (2.1). But it is difficult to interpret biologically these conditions. Numerical simulations indicates that E∗ is stable for d < 2.2 and loses its stability at d = 2.2 and the system undergoes Neimark-Sacker bifurcation when the death rate d exceeds the value 2.2 (see Fig. 1). We have observed that when intraspecific competition among the prey species is low and the death rate of the predator exceeds the value 2.05, system admits flip bifurcation (see Fig. 2) DYNAMICS OF A DISCRETE PREDATOR-PREY SYSTEM 57 Figure 2. Bifurcation diagram for prey and predator populations with d for fixed values r = 4.5,k = 1,b = 1,c = 1,p = 4,q = 1,h = 1,α = 0.1 and β = 0.01. follows from the Proposition 5. We also note that when the fear factor k crosses the critical value 2.15, the system undergoes Neimark-Sacker bifurcation (see Fig. 3). The chaotic nature of the system is nicely controlled by the state feedback control strategy (see Fig. 4). Numerical simulation exhibits that feedback control mechanism can dominate chaos to unstable fixed point strongly and ultimately stability of the system is achieved. Acknowledgment. The author is grateful to the anonymous reviewer and the Associate Editor for their helpful comments for improving this paper. 58 D. MUKHERJEE Figure 3. Bifurcation diagram for prey and predator populations with k for fixed values r = 4.5,d = 1.6,b = 1,c = 1,p = 1.8,q = 1,h = 1,α = 0.1 and β = 0.01. 0 0.5 1 1.5 2 2.5 3 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Prey P re d a to r 0 0.5 1 1.5 2 2.5 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Prey P re d a to r Figure 4. Phase diagram of system (2.1) for r = 4.5,k = 1,b = 1,c = 1,d = 1,p = 4,q = 1,h = 1,α = 0.1 and β = 0.01. with initial values (x0,y0) = (0.1, 0.1) in the left panel and controlled system (4) for c1 = 0.3 and c2 = −1.2 in the right panel. DYNAMICS OF A DISCRETE PREDATOR-PREY SYSTEM 59 References [1] X. Wang, L. Y. Zanette and X. Zou, Modelling the fear effect in predator-prey interactions, Journal of Mathematical Biology 73 (2016), 1179-1204. [2] D. Mukherjee, Study of fear mechanism predator-prey system in the presence of competitor for the prey, Ecological Genetics and Genomics 15 (2020), 100052. [3] D. Mukherjee, Role of fear in predator-prey system with intraspecific competition, Mathematics and Computers in Simulation 177 (2020), 263-275. [4] P. Panday, N. Pal, S. Samanta and J. Chattopadhyay, Stability and bifurcation analysis of a three-species food chain model with fear, International Journal of Bifurcation and Chaos. 28 (2018). [5] N. Pal, S. Samanta and J. Chattopadhyay, Fear effect in prey and hunting cooperation among predators in a Leslie- Gower model, Mathematical Biosciences and Engineering 16 (2019), 5146-5179. [6] S. Sasmal, Population dynamics with multiple Allee effects induced by fear factors- a mathematical study on prey- predator interactions, Applied Mathematical Modelling 64 (2018), 1-14. [7] S. Sasmal and Y. Takeuchi, Dynamics of a predator-prey system with fear and group defence, Journal of Mathematical Analysis and Applications 481 (2020) 123471. [8] J. Wang, Y. Cai, S. Fu and W. Wang, The effect of the fear factor on the dynamics of a predator-prey model incorporating the prey refuge, Chaos 29 (2019) 083109. [9] X. Wang and X. Zou, Modelling the fear effect in predator-prey interactions with adaptive avoidance of predators, Bulletin of Mathematical Biology. 79 (2017), 1-35. [10] H. Zhang, Y. Cai, S. Fu and W. Wang, Impact of the fear effect in a prey-predator model incorporating a prey refuge, Applied Mathematics and Computation 356 (2019), 328-337. [11] F. Chen, X. Guan, X. Huang and H. Deng, Dynamic behaviors of a Lotka-Volterra type predator-prey system with Allee effect on the predator species and density dependent birth rate on the prey species, Open Mathematics 17 (2019), 1186-1202. [12] H. N. Agiza, E. M. Elabbasy, EI-Metwally and A. A. Elasdany, Chaotic dynamics of a discrete predator-prey model with Holling type II, Nonlinear Analysis. Real World Applications 10 (2009), 116-129. [13] Q. Din, Complexity and chaos control in a discrete time prey-predator model, Communications in Nonlinear Science and Numerical Simulation 49 (2017), 113-134. [14] M. E. Elettreby, T. Nabil and A. Khawagi, Stability and bifurcation analysis of a discrete predator-prey model with mixed Holling interaction, Computer Modeling in Engineering and Sciences 122 (2020), 907-921. [15] Z. He and B. Li, Complex dynamic behavior of a discrete-time predator-prey system of Holling-III type. Advances in Difference Equations 2014 (2014), 1-12. [16] P. Santra, G. S. Mahapatra and G. Phaijoo, Bifurcation and chaos of a discrete predator-prey model with Crowley- Martin functional response incorporating proportional prey refuge, Mathematical Problems in Engineering 2020 (2020), 1-18. [17] H. Seno, A discrete prey-predator model preserving the dynamics of a structurally unstable Lotka-Volterra model, Journal of Difference Equations and Application 13 (2007), 1155-1170. [18] L. Berezansky, E. Braverman and L. Idels, Nicholson’s blowflies differential equations revisited: main results and open problems, Applied Mathematical Modelling 34 (2010), 1405-1417. [19] Y. Neimark, On some cases of periodic motions depending on parameters. Dokl. Acad. nauk. SSSR 129 (1959), 736-739. [20] R. J. Sacker, On invariant surfaces and bifurcation of periodic solutions of ordinary differential equations, Ph.D. thesis, IMM-NYU 333, Courant Institute, New York University, 1964. [21] S. Elaydi, Discrete Chaos with Applications in Science and Engineering, Boca Raton, Chapman and Hall, 2008. [22] K. Kundu, S. Pal, S. Samanta, A. Sen and N. Pal, Impact of fear effect in a discrete-time predator-prey system, Bulletin of Calcutta Mathematical Society. 110 (2018), 3-12. Department of Mathematics, Vivekananda College, Thakurpukur, Kolkata-700063, India. Email address: mukherjee1961@gmail.com 1. Introduction 2. Discrete model 3. Fixed points and their nature 3.1. Bifurcation around the interior fixed point 4. Chaos control 5. Numerical Simulation 6. Discussion References