Academic Journal of Science and Technology ISSN: 2771-3032 | Vol. 3, No. 3, 2022 119 Complicate Dynamics of a Discrete Predator‐prey Model with Competition Between Predators Shaosan Xia1, a, Xianyi Li1, b, * 1Department of Big Data Science, School of Science, Zhejiang University of Science and Technology, Hangzhou, 310023, China aE-mail: ssxia96@163.com, bmathxyli@zust.edu.cn *Correspondence author: Xianyi Li Abstract: We consider a discrete predator-prey model with competition between predators in this paper. By simplifying the corresponding continuous predator-prey model, and using the semidiscretization method to obtain a new discrete model, we discuss the existence and local stability of nonnegative fixed points of the new discrete model. What’s more important, by using the bifurcation theory, we derive the sufficient conditions for the occurrences of Neimark-Sacker bifurcation and the stability of closed orbit bifurcated. Finally, the numerical simulation are presented, which not only illustrate the existence of NeimarkSacker bifurcation but also reveal some new dynamic phenomena of this model. Keywords: Discrete predator-prey system with competition, Semidiscretization method, Neimark-Sacker bifurcation. 1. Introduction The dynamic relationship between predator and prey [1-4] is a kind of interaction between different species in the ecosystem. In population biology, Lotka Volterra (L-V) predator-prey model [5,6] is one of the most famous models, which was respectively proposed by A. Lotka in 1925 and V. Volterra in 1926. This model has become the standard basis for many subsequent models in many different fields. Since Lotka and Volterra’s pioneering work, the predator-prey model in the ecosystem has been an important research role from an ecological point of view. In random and deterministic environments, mathematical modeling has an increasing impact on theoretical ecology. In the past few decades, a lot of empirical and theoretical work on predator-prey models has been carried out. Functional response plays an important role in system modeling. Functional response is a function of the number of prey consumed by each predator per unit time. This nutrient function determines the dynamics of the system, such as the stability and bifurcation, etc. The predator-prey models studied in the past half century are given many different functional responses. Enrichment paradox has also become one of the important topics in predator-prey model [7-9]. It is pointed out that enriching the predatorprey system by increasing the carrying capacity of a predator-prey system to prey will lead to the increase of predator equilibrium density, but not lead to the increase of prey equilibrium density, and will eventually destroy the positive balance. Intuitively, this increases the probability of random extinction of predator and prey species. The interspecific and intraspecific interactions in a L-V model can not reflect the real experimental phenomena or the specific characteristics of different groups. Therefore, the original L-V model has been developed and improved by combining more realistic and important biological factors and relationships. The following improved Lotka Volterra predator-prey system is also called Wangersky-Cunningham model [10] (1 ) , . du u u muv d k dv emuv dv d             (1) Most predator-prey models are developed on the basis of trophic function, which is a function of prey density. By introducing intraspecific competition in the predator population, the above model (1) is further modified and studied; for instance, Pielou [11] modified it and took the following form: 2 (1 ) , , du u u muv d k dv emuv dv hv d              (2) where u(τ) and v(τ) stand for the prey and predator density, respectively, at time τ, and the parameters γ,k,m,e,d,h are positive constants, which respectively stand for prey intrinsic growth rate, carrying capacity of the environment to prey, consumption rate, conversion factor of biomass due to change of food level, predator death rate, predator interspecies competition degree. In order to simplify the analysis of system (2), we use the scaling x=u/k, y = mv/γ and t = γτ to nondimensionalize the system (2) to derive the following dimensionless system: (1 ), ( ), dx x x y dt dy y b ax ry dt            (3) where 120 , , . emk d h a b r m     We now use the semidiscretization method, which has been applied in many studies [12-15], to derive the discrete model of system (3). For this, suppose that [t] denotes the greatest integer not exceeding t. Consider the following semidiscretization version of system (3). (4) It is easy to see that the system (4) has piecewise constant arguments, and that a solution (x(t),y(t)) of the system (4) for t ∈ [0,+∞) possesses the following natures: 1. on the interval [0,+∞), x(t) and y(t) are continuous; 2. when t ∈ [0,+∞) except for the points and exist everywhere. The following system can be obtained by integrating (4) over the interval [n,t] for any t∈[n,n+1) and n = 0,1,2,··· 1( ) ( ), ( ) ( ), n n n n x y n b ax ry n x t x e t n y t y e t n            (5) where xn = x(n) and yn = y(n). Letting ( 1)t n   in (5) produces 1 1 1 , , n n n n x y n n b ax ry n n x x e y y e            (6) where a,b,r > 0. We mainly study the properties of system (1.6) in the sequel. Before we analyze the fixed points of the system (6), we recall the following lemma (see [13, pp1628], [16, pp422]). Lemma 1.1. Let F(λ) = λ2+Pλ+Q, where P and Q are two real constants. Suppose λ1 and λ2 are two roots of F(λ) = 0. Then the following statements hold. If F(1) > 0, then |λ1| < 1 and |λ2| < 1 if and only if F(−1) > 0 and Q < 1; λ1 = −1 and λ2 ≠ −1 if and only if F(−1) = 0 and P ≠ 2; |λ1| < 1 and |λ2| > 1 if and only if F(−1) < 0; |λ1| > 1 and |λ2| > 1 if and only if F(−1) > 0 and Q > 1; λ1 and λ2 are a pair of conjugate complex roots with |λ1| = |λ2| = 1 if and only if −2 < P < 2 and Q = 1; λ1 = λ2 = −1 if and only if F(−1) = 0 and P = 2. If F(1) = 0, namely, 1 is one root of F(λ) = 0, then the another root λ satisfies |λ| = (<,>)1 if and only if |Q| = (<,>)1. If F(1) < 0, then F(λ) = 0 has one root lying in (1,∞). Moreover, (iii.1) the other root λ satisfies λ < (=)−1 if and only if F(−1) < (=)0; (iii.2) the other root −1 < λ < 1 if and only if F(−1) > 0. 2. Existence and Stability of Fixed Points In this section, we first consider the existence of fixed points and then analyze the local stability of each fixed point of the system (6). The fixed points of the system (6) satisfy 1 ,x y b ax ryx xe y ye      Considering the biological meanings of the system (6), one only takes into account its nonnegative fixed points. Thereout, one notices that the system (6) has and only has three nonnegative fixed points E0 = (0,0),E1 =(1,0) and 2 ( , ) b r a b E a r a r      for a-b>0. The Jacobian matrix of the system (6) at any fixed point E(x,y) takes the following form . The characteristic polynomial of Jacobian matrix J(E) reads F(λ) = λ2 − pλ + q, where 1( ( )) (1 ) (1 ) ,x y b ax ryp Tr J E x e ry e         1 ( 1) ( 1)( ( )) [1 ( ) ] .b a x r yq Det J E x ry a r xy e           For the stability of fixed points E0, E1 and E2, we can easily get the following Theorems 2.1-2.3 respectively. Theorem 2.1. The fixed point E0 = (0,0) of the system (6) is a saddle. Proof. The Jacobian matrix J(E0) of the system (6) at the fixed point E0 = (0,0) is given by . Obviously, |λ1| = e > 1 and |λ2| = e‾b < 1, so E0 = (0,0) is a saddle. Theorem 2.2. The following statements about the fixed point E1 = (1,0) of the system (6) are true. If b < a, then E1 is a saddle. If b = a, then E1 is non-hyperbolic. If b > a, then E1 is a stable node. Proof. The Jacobian matrix of the system (1.6) at E1 = (1,0) is . Obviously, λ1 = 0 and λ2 = ea‾b. Note |λ1| < 1 is always true. If b < a, then |λ2| > 1, so E1 is a saddle; if b = a, then |λ2| = 1, therefore E1 is non-hyperbolic; if b > a, implying |λ2| < 1, then E1 is a stable node, namely, a sink. The proof is complete. Theorem 2.3. When 0a b  , 2 ( , ) b r a b E a r a r      is a positive fixed point of the system (1.6). Let 0 2 (2 )( ) 2 a b a b r a b       for 2a b  and r1 = b(a − b − 1), then the following table statements are true about the positive fixed point E2. 121 Table 1. Properties of the fixed point E2. Conditions Eigenvalues Properties 0 < a − b ≤ 1 |λ1| < 1,|λ2| < 1 sink 1 < a − b ≤ 2 r < r1 |λ1| > 1,|λ2| > 1 source r = r1 |λ1| = |λ2| = 1 non-hyperbolic r > r1 |λ1| < 1,|λ2| < 1 sink 2 < a – b 0 < r < r1 |λ1| > 1,|λ2| > 1 source r = r1 |λ1| = |λ2| = 1 non-hyperbolic r1 < r < r0 |λ1| < 1,|λ2| < 1 sink r = r0 λ1 = −1,λ2 ≠−1 non-hyperbolic r > r0 |λ1| < 1,|λ2| > 1 saddle 0 < r < r1 |λ1| > 1,|λ2| > 1 source r = r1 λ1 = λ2 = −1 non-hyperbolic r > r1 |λ1| < 1,|λ2| > 1 saddle 2( 2)b a b   0 < r < r0 |λ1| > 1,|λ2| > 1 source r = r0 λ1 = −1,λ2 ≠ −1 non-hyperbolic r > r0 |λ1|< 1, |λ2| > 1 saddle 3. Bifurcation Analysis In this section, we are in a position to use the bifurcation theorem to analyze the local bifurcation problems of the fixed points E2. For related work, refer to [17-21]. 3.1. For fixed point 2 ( , ) b r a b E a r a r      When r = r1 = b(a − b − 1), Theorem 2.4 with Lemma 1.1 (i.5) shows that F(1) > 0, F(−1) > 0, −2 < p < 2 and q = 1, so λ1 and λ2 are a pair of conjugate complex roots with |λ1| = |λ2| = 1. At this time we derive that the system (6) at the fixed point E2 can undergo a Neimark-Sacker bifurcation in the space of parameters 3( , , ) {( , , ) |1 2, 0}Ea b r S a b r R a b r         . In order to show the process clearly, we carry out the following steps. The first step. Take the changes of variables un = xn − x0,vn = yn −y0, which transform fixed point E2 = (x0,y0) to the origin O(0,0), and the system (1.6) into 0 0 0 0 1 1 0 0 ( ) ( ) 1 0 0 ( ) , ( ) . n n n n u x v y n n b a u x r v y n n u u x e x v v y e y                    (7) The second step. Give a small perturbation r* of the parameter r, i.e., r* = r − r1, then the perturbation of the system (3.1) can be regarded as follows 0 0 * 0 1 0 1 1 0 0 ( ) ( )( ) 1 0 0 ( ) , ( ) . n n n n u x v y n n b a u x r r v y n n u u x e x v v y e y                     (8) The corresponding characteristic equation of the linearized equation of the system (7) at the equilibrium point (0,0) can be expressed as F(λ) = λ2 − p(r*)λ + q(r*) = 0, where , and . It is easy to derive p2(r*) − 4q(r*) < 0 when r* = 0, then the two roots of F(λ) = 0 are as follows , moreover which implies The occurrence of Neimark-Sacker bifurcation requires the following conditions to be satisfied * * 1,2 * 0 | ( ) | ( .1) 0;( ) | r d r H dr    1,2( .2) (0) 1, 1,2,3,4.iH i   Since * * 0 1 ( 1) ( ) 1 1 | r b a b p r b       and * * 0 ( ) 1| r q r   ,we have 2 2 1,2 2 ( 2) 4 ( 2) (0) 2(1 ) b a b i ab b a b b           , then it is easy to derive 1,2 (0) 1m  for all m = 1,2,3,4. According to [22, pp517-522], they satisfy all of the conditions for Neimark-Sacker bifurcation to occur. 122 4. Numerical Simulation In this section, we use the bifurcation diagrams and Lyapunov exponents of the system (6) to illustrate our theoretical results and further reveal some new dynamical behaviors to occur as the parameters vary by Matlab software. Fix the parameter values a = 1.5,b = 0.2, let r ∈ (0,0.3) . Figure 1(a) shows the bifurcation diagram of (r,x)-plane, from which the fixed point E2 is unstable when r < r1 = 0.06 whlie stable when b > b0. Hence, the Neimark-Sacker bifurcation occurs at the fixed point E2 = (0.167,0.833) when r = r1, whose multipliers are λ1,2 = 0.892±0.453i with |λ1,2| = 1.The corresponding maximum Lyapunov exponent diagram of the system (1.6) is plotted in Figure 1(b). (a) r ∈ (0,0.3)(b) r ∈ (0,0.3) Figure 1. Bifurcation of the system (1.6) in (r,x)-plane and Maximal Lyapunov exponent. 5. Discussion and Conclusion In this paper, we discuss the dynamical behaviors of a predator-prey model (1.6) with competition between predators. Under the given parametric conditions, we completely show the existence and stability of three nonnegative equilibria E0 = (0,0), E1 = (1,0) and 2 ( , ) b r a b E a r a r      . Then we derive the sufficient conditions for its Neimark-Sacker bifurcation to occur. Meanwhile, it is clear that the positive equilibrium E2 = (x0,y0) is asymptotically stable when r > r1 = b(a − b − 1) and unstable when r < r1 under the condition 1 < a − b ≤ 2. Hence, the system (1.6) undergoes a bifurcation which has been shown to be a Neimark-Sacker bifurcation when the parameter r goes through the critical value r1. Finally, numerical simulations confirm the theoretical analysis results of the system (1.6). Acknowledgment This work is partly supported by the National Natural Science Foundation of China (61473340), the Distinguished Professor Foundation of Qianjiang Scholar in Zhejiang Province (F703108L02), and the Natural Science Foundation of Zhejiang University of Science and Technology (F701108G14). References [1] Md, S. & Rana, S. [2019] “Dynamics and chaos control in a discrete-time ratio-dependent Holling-Tanner model,” J.Egypt.Math.Soc. 2019, doi:10.1186/s42787-019-0055-4. [2] Khan, A. Q. [2016] “Neimark-Sacker bifurcation of a twodimensional discrete-time predator-prey model,” Adv.Differ.Equ. 2016, doi:10.1186/s40064-015-1618-y. [3] Rodrigo, C., Willy, S. & Eduardo, S. [2017] “Bifurcations in a predatorprey model with general logistic growth and exponential fading memory,” Appl.Math.Model. 45, 134–147. [4] Holling, C. [1959] “The functional response of predator to prey density and its role in mimicry and population regulation,” Mem.Entomol.Soc.Can. 91, 385-398. [5] Berryman, A. A., Gutierrez, A. P. & Arditi, R. [1995] “Credible, Parsimonious and useful predator–prey models - A reply to Abrams, Gleeson, and Sarnelle,” Entomological Society of America. 76, 1980–1985. [6] Saunders, M. C. [2006] “Complex Population Dynamics: A Theoretical/Empirical Synthesis,” Entomological Society of America. 35, 1139-1139. [7] Akcakaya, HR., et al. [1995] “Ratio-dependent predation: an abstraction that works,” Ecology. 76, 995–1004. [8] Cosner, C., et al. [1999] “Effects of spatial grouping on the functional response of predators,” Theor.Popul.Biol. 56, 65–75. [9] Gutierrez, AP. [1992] “Physiological basis of ratio-dependent predatorprey theory: the metabolic pool model as a paradigm,” Ecology. 73, 1552–1563. [10] May, R. [1973] “Stability and complexity in model ecosystems with a new introduction by the author,” Princeton University Press. [11] Pielou, EC. [1969] “An introduction to mathematical ecology,” Biometrical Journal. 13, doi:10.1002/bimj.19710130308. [12] Din, Q. [2017] “Complexity and chaos control in a discrete- time prey-predator model,” Commun.Nonliner.Sci.Numer. Simul. 49, 113-134. [13] Li, W. & Li, X. Y. [2018] “Neimark–Sacker bifurcation of a semi-discrete hematopoiesis model,” J.Appl.Anal.Comput. 8, 1679–1693. [14] Hu, Z. Y., Teng, Z. D. & Zhang, L. [2011] “Stability and bifurcation analysis of a discrete predator-prey model with nonmonotonic func-tional response,” Nonlinear. Anal. Real. World. Appl. 12, 2356–2377. [15] Wang, C. & Li, X. Y. [2015] “Further investigations into the stability and bifurcation of a discrete predator-prey model,” J.Math.Anal.Appl. 422, 920–939. [16] Wang, C. & Li, X. Y. [2014] “Stability and Neimark–Sacker bifurcation of a semi-discrete population model,” J.Appl.Anal.Comput. 4, 419–435. 123 [17] Li, M. S., Zhou, X. L. & Xu, J. M. [2020] “Dynamic properties of a discrete population model with diffusion,” Adv.Differ.Equ. 2020, doi:10.1186/s13662-020-03033-w. [18] Gallay, T. [1993] “A center-stable manifold theorem for differential equations in Banach spaces,” Commun.Math.Phys. 152, 249–268. [19] Liu, W., Cai, D. & Shi, J. [2018] “Dynamic behaviors of a discretetime predator–prey bioeconomic system,” Adv.Differ.Equ. 2018, doi:10.1186/s13662-018-1592-0. [20] Jorba, A. & Masdemont, J. [1999] “Dynamics in the center manifold of the collinear points of the restricted three body problem,” Physica D : Nonlinear Phenomena. 132, 189–213. [21] Zhao, M., Xuan, Z. & Li, C. [2016] “Dynamics of a discrete- time predator-prey system,” Adv.Differ.Equ. 2016, doi:10.1186/s13662016-0903-6. [22] Winggins, S. [2003] “Introduction to Applied Nonlinear Dynamical Systems and Chaos,” Springer-Verlag, New York.