Multiple bifurcations of a discrete modified Leslie-Gower predator-prey model

: In this paper, we work on the discrete modified Leslie type predator-prey model with Holling type II functional response. The existence and local stability of the fixed points of this system are studied. According to bifurcation theory and normal forms, we investigate the codimension 1 and 2 bifurcations of positive fixed points, including the fold, 1:1 strong resonance, fold-flip and 1:2 strong resonance bifurcations. In particular, the discussion of discrete codimension 2 bifurcation is rare and di ffi cult. Our work can be seen as an attempt to complement existing research on this topic. In addition, numerical analysis is used to demonstrate the correctness of the theoretical results. Our analysis of this discrete system revealed quite di ff erent dynamical behaviors than the continuous one.


Introduction
Mathematical biology has gradually become one of the hot spots and frontier topics in recent decades.The advent of qualitative and quantitative analysis of biological models has made it possible to move from understanding the underlying mechanisms of principles to analyzing them scientifically.It provides a guarantee for making predictions about some biological phenomena.Analyzing species' interactions is one of the research directions in mathematical biology.Generally, a network of interacting species, called a trophic web, forms a complex structure.Each population of interacting species is affected by the others [1].Considering the interaction of two species, when the growth rate of one population increases and the other's decreases, these two species are in a predatory relationship.Within a region, the two species with a predatory relationship may coexist or experience the case of species extinction.These situations are of interest to biologists and mathematicians.Thus they construct predator-prey models to research the predatory relationship.
The general Leslie type predator-prey model is one of the most classical predator-prey models; it has the following form where x is the population of the prey and ȳ is the population of the predator.The parameters r and s are intrinsic growth rates for the prey and predator populations, respectively.k 1 reflects the carrying capacity for prey and p( x, ȳ) is the functional response function which represents how the predator deals with changes in the prey.h x is the carrying capacity for the predator which is proportional to the available prey.Different functional response functions p( x, ȳ) have various effects on the dynamical properties of system (1.1).The authors of [2][3][4] have considered that p( x, ȳ) is replaced by Holling type III functional response.The case that p( x, ȳ) is Holling type IV functional response has been discussed by the authors of [5][6][7].The authors of [8,9] have characterized the case that the functional response function p( x, ȳ) is related to the predator and prey populations.In 2003, Aziz-Alaoui and Okiye [10] modified system (1.1) by introducing Holling type II functional response and predator's other food sources.This model is given as follows: where mȳ x+n is Holling type II functional response.h x + k 2 is the carrying capacity for the predator population and the term k 2 stands for other food sources for the predator population.All parameters of this system are positive.
Ever since system (1.2) was proposed, it has attracted many interested researchers.Giné and Valls [11] discussed the nonlinear oscillations in system (1.2).Lin and Jiang [12] applied n = k 2 h to system (1.2), combined with stochastic perturbation.And, Xie et al. [13] investigated the case that n = k 2 h with the linear harvesting of two species.In particular, Xiang et al. [14] applied the following scaling to system (1.2): x = k 1 x, ȳ = hk 1 y, τ = t r and they obtained the following: where k = mh r , b = s r , a 1 = n k 1 and a 2 = k 2 hk 1 .They analyzed the codimension 2 and 3 bifurcations of system (1.3).Moreover, in their work, a lot of valuable results were found by adding the changing environment.
Besides the well-known continuous models, such as those described in [15][16][17], the discrete ones have a profound influence and are equally noteworthy.The discrete systems are more applicable to populations with non-overlapping generations, and they have many unique phenomena in addition to the dynamics corresponding to the continuity.The bifurcation of discrete systems plays a key role.When the bifurcation parameter is slightly perturbed near the critical value, the topology of the discrete system changes.Then, it will exhibit a series of dynamic changes that deserve attention.Especially, the study of the discrete codimension 2 bifurcations is more difficult and should receive attention.In the analysis of codimension 2 bifurcations, two independent coefficients of the difference equation are selected as bifurcation parameters.Crossing the two-dimensional bifurcation curves can cause the occurrence of some corresponding codimension 1 bifurcations.There are many attractive results about the bifurcations in discrete systems, such as those described in [18][19][20][21][22].Among them, the authors of [21] discretized system (1.2) and analyzed its dynamics, which involved the codimension 1 bifurcations and Marotto's chaos.To the best of our knowledge, there are no works about codimension 2 bifurcations of the discrete form of system (1.2).Thus, this issue is the major research topic of our work.
There are many methods of obtaining the discrete form of continuous systems, such as the Runge-Kutta, Taylor series and linear multistep methods.Although these are higher-accuracy methods, they use more calculations, more past values or derivatives [23].However, the Euler method is a traditional and simple way.In particular, the stability of the Euler integrator is associated with the value of the step size.When the step size of the Euler method is large, it may be possible to obtain dynamics that are very different from those of the original continuous system.Chaos is also related to this artificially induced instability [20].
Therefore, we apply the same scaling as in [14] to system (1.2) and use the Euler method.Then the following model is obtained: where d is the step size and all parameters are positive.We provide the stability and bifurcation analysis of system (1.4), and this paper is organized as follows.In Section 2, the existence and local stability of fixed points are investigated mainly through the use of the stability theory and center manifold theorem.
In Section 3, we analyze the occurrence of codimension 1 and 2 bifurcations of the interior fixed points, including fold bifurcation, 1:1 and 1:2 strong resonance bifurcations and fold-flip bifurcation.
Moreover, the results are demonstrated through numerical analysis in Section 4. A brief conclusion is shown in Section 5.

The existence and stability of fixed points
From the following equations: it is easy to know that system (1.4) has the trivial fixed point P 0 (0, 0) and the semitrivial fixed points P 1 (1, 0), P 2 (0, a 2 ).For the positive fixed points, we have the following assumptions.
Next, we consider the stability of these fixed points and have the following propositions.
Proposition 1.The fixed point P 0 is unstable and P 1 is a saddle.
Proof.J(P 0 ) and J(P 1 ) are the Jacobian matrices of system (1.4) at P 0 and P 1 , respectively, where and Apparently, the eigenvalues associated with P 0 satisfy that k , then P 2 is semi-stable from the left.Proof.The Jacobian matrix of system (1.4) at P 2 is (2.4) a 1 and λ P 2 ,2 = 1 − bd are the associated eigenvalues, where λ P 2 ,2 < 1.It follows from k , then the fixed point P 2 is a saddle (or stable).
When a 2 = a 1 k , let u n = x n and v n = y n − a 1 k .Thus, system (1.4) become as follows: (2.5) Applying the invertible transformation to system (2.5), we get where By the center manifold theorem, we assume that Naturally, we conclude that F ′ (0) = 1 and a 1 > 0. Hence, by the theory in [24], P 2 is an unstable fixed point and it is semi-stable from the left.□ Proposition 3. The fixed point P 3 is non-hyperbolic.
Proof.The Jacobian matrix of system (1.4) at P 5 is According to Lemma 4, we come to the result as follows.If H(−1) = 1 + T r(J(P 5 )) + Det(J(P 5 )) and hold, then P 5 is stable.□ Remark.The example with specific parameters given in Section 4 can intuitively illustrate this proposition.The conditions for P 4 stability are similar, so we omit them.

Bifurcations at P 3
In this subsection, we discuss the codimension 1 and 2 bifurcations at P 3 .The coefficients that are not listed explicitly in Subsection 3.1 will be given in Appendices A, B and C.
Proof. a 2 is chosen as the bifurcation parameter and a new variable.Let a 2 = (1−k−a 1 ) 2 +4a 1 4k + a * f , where a * f is a sufficiently small perturbation.Then we transform P 3 into (0, 0) by taking U n = x n − x 3 and V n = y n − y 3 .Thus, system (1.4) can be expanded at the origin and we obtain the following form where . Applying the following invertible transformation to system (3.1): where According to the center manifold theorem, 3 ) is assumed and the equation 1−k+a 1 0 are calculated.
Theorem 7. If the conditions (k−1) 2 a 1 (3k−1), a 1 (a 1 +4k) and detD γ V(0) 0 hold, then system (1.4) undergoes a 1:1 strong resonance bifurcation at P . When |γ| is sufficiently small, the following local dynamics exist: 1) There is a fold bifurcation that occurs on in the γ-space.2) For one of the fixed points born at the fold bifurcation of 1), there is a Neimark-Sacker bifurcation that occurs on in the γ-space.Moreover, the invariant closed curve of the Neimark-Sacker bifurcation is attracting (repelling 3) There is a homoclinic bifurcation at which the stable and unstable manifolds of the saddle point born at the fold bifurcation of 1) occurs on two curves H 1,2 and has the asymptotic form: in the γ-space.The distance between two homoclinic bifurcation curves H 1,2 is exponentially small with regard to |γ|.The above-described curves and phase portraits are shown schematically in Figure 1.
and vn = y n − y 3 ; then, we can expand system (1.4) at (0, 0) as follows 1−k+a 1 and J R1 (γ) is the Jacobian matrix of system (3.4) at the origin.When b * R1 = a * R1 = 0, the two eigenvalues of J R1 (0) are λ1 = λ2 = 1.We can select the following linearly independent eigenvectors (generalized eigenvectors): which satisfy the following equations where the symbol ⟨ * , * ⟩ stands for the standard scalar product in R 2 .Then, we can construct the invertible transformation ūn vn = ln p0 to simplify the linear part of system (3.4).And, the following equations are deduced: Under these new coordinates, system (3.4) becomes where and ā00 (0 According to Lemma 9.6 [27], system (3.8) can be written in the following form if |γ| is sufficiently small: where φ 1 γ ( ln , mn ) is the flow of the planar system where By calculation, the nondegeneracy condition d20 (0 And, the following equivalence relationship is deduced: Suppose that conditions (3.11) and (3.12) hold.According to Lemma 3.2 [28], we get the new system by applying the analytic changes of coordinates and a scaling of time to system (3.10),where Further, a series of complex calculations provide that If conditions (3.11) and (3.12) and the transversality condition detD γ V(0) 0 hold, then system (3.13) is the versal unfolding of the Bogdanov-Takens singularity of codimension 2. Referring to Proposition 3.1 [28], the existence of 1:1 strong resonance bifurcation of system (1.4) is obtained.
k(a 1 +k−1) and d . The fold-flip bifurcation of system (1.4) occurs at P 3 when {b, a 2 } varies in a sufficiently small neighborhood of . Furthermore, there are the following local dynamical behaviors taking place: 2) There exists a nondegenerate flip bifurcation on the curve ( X, Ŷ, α 1 ) = (0, 0, 0).

Bifurcations at P
, k + a 1 < 1 , then system (1.4) has two positive fixed points P 4,5 (x 4,5 , y 4,5 ).When a 2 decreases to a 1 k (or is less than it), P 4 is not positive and system (1.4) has the unique interior fixed point P 5 .Because of the different levels of details, in this work, we consider the bifurcations at The coefficients that are not listed explicitly in this subsection will be given in Appendix D.
The planar system is smoothly equivalent to the following: where If Ã1 (0), B1 (0) 0, the bifurcations of system (3.30) can be analyzed by using the following system: where s 3 = sign( Ã1 (0)).By applying the theory in [27] and [31], we summarize our analysis into Theorem 9.The schematic diagrams of bifurcation curves and phase portraits can be found on pp.444-446 in [27].□

Numerical analysis
In this section, we use some cases with specific values to explain our theoretical analysis.Let a 1 = 0.1, 2 = 0.4, k = 0.3, b = 1.5 and d = 0.7.We find that the conditions H(1) = 0.330403 > 0, H(−1) = 1.95787 > 0 and Det(J(E 5 )) − 1 = −0.855866< 0 hold.Therefore, P 5 (0.564575, 0.964575) is the stable interior fixed point of system (1.4). Figure 2 is the time-series diagram of the prey and predator populations.As the number of iterations increases, the populations of both species are constant, which implies that P 5 is stable.

Fold bifurcation
Let a 1 = 0.2, k = 0.4, d = 0.8 and b = 0.5.a 2 = 0.6 is the critical value such that the fold bifurcation of system (1.4) occurs.Figure 3 is the fold bifurcation diagram.By analyzing this diagram, we find that there is no positive fixed point of system (1.4) when a 2 > 0.6.If a 2 = 0.6, there is a unique interior fixed point P 3 (0.2, 0.8).And if a 2 < 0.6, there are two interior fixed points P 4 and P 5 bifurcate from P 3 .According to the eigenvalues of P 4 and P 5 , we know that P 5 is stable and P 4 is unstable.

1:1 strong resonance bifurcation
Fixing d = 0.8, k = 0.7 and a 1 = 0.1.The critical values of the bifurcation parameters are a R1 = 0.15714286 and b R1 = 0.35 when λ P 3 ,1 = λ P 3 ,2 = 1.Then we calculate that the conditions (k−1) 2 = 0.03 a 1 (3k−1) = 0.11, (k−1) 2 = 0.03 a 1 (a 1 +4k) = 0.29 and detD γ V(0) = −3770.470 are satisfied.Hence, system (1.4) undergoes a 1:1 strong resonance bifurcation at P 3 (0.1, 0.25714286) as {b, a 2 } varies in a small neighborhood of {0.35, 0.15714286}.The software package MatContM has been applied to confirm our analysis.For more details of this package, [32,33] are available for reference.The Neimark-Sacker (NS) and fold (LP) curves are shown in Figure 4(a) and the symbol R1 denotes the 1:1 strong resonance bifurcation point.In Figure 4(b), we plot the phase portrait with a 2 = 0.13 and b = 0.35.Based on observation of this subfigure, there exists an invariant closed curve.Figure 4(c) is the time series diagram of the prey x and predator y populations with the same specific parameters as Figure 4(b).The blue and magenta points represent the populations of prey and predator, respectively.Both species exist in periodic oscillations due to the invariant closed curve.k(a 1 +k−1) = −16.66667and d a 1 (k+4+4a 1 ) 2k(a 2 1 −(k−1) 2 ) = −5.75.By Theorem 8, the fold-flip bifurcation of system (1.4) occurs at P 3 (0.2, 1.33333) when {b, a 2 } varies in a sufficiently small neighborhood of {2.63, 1.13333}.Figure 5(a) illustrates our analysis and LPPD is the symbol of the fold-flip bifurcation point.Then, let a 2 = 1.15 and b = 3.5.The phase portrait is shown in Figure 5(b).At this point, system (1.4) is in a chaotic state in which the predator and prey populations cannot coexist stably.In Figure 5(c), the positive values of the corresponding maximum Lyapunov exponents also explain the existence of chaos.

1:2 strong resonance bifurcation
Suppose that d = 2.3, a 1 = 0.2 and a 2 = 0.4.k R2 = 0.336956 and b R2 = 1.47749 are the critical bifurcation parameters such that λ1 27343 0, Ã1 (0) = −178.9360 and B1 (0) = 26.76840 hold.According to Theorem 9, the 1:2 strong resonance bifurcation of system (1.4) takes place at P 5 (0.576225, 0.976225), when {k, b} varies in a small neighborhood of {0.336956, 1.47749}.The NS and PD curves are plotted in Figure 6(a) which illustrates the existence of this bifurcation.The symbol R2 stands for the 1:2 strong resonance bifurcation point.In this example, we take the step size of the Euler method as 2.3.The choice of a large step size changes the dynamical behaviors of the original continuous system.Moreover, we can also find fold-flip and generalized period-doubling bifurcation points in Figure 6(a), where the symbol GPD represents the period-doubling bifurcation point.Figure 6(b) is the phase portrait of x and y, when k = 0.33 and b = 1.47.From this subfigure, we can observe that due to the existence of this bifurcation, the dynamical properties of system (1.4) change in the neighborhood of the critical parameters and produce an invariant closed curve.
Then, select d and b as the bifurcation parameters.Let k = 0.6, a 1 = 0.2 and a 2 = 0.3.By numerical simulation, we know that a 1:2 strong resonance bifurcation of system (1.4) exists at P 5 (0.27320508, 0.57320508) when {b, d} varies in a sufficiently small neighborhood of {1.0728398, 4.3176514}.The PD and NS curves in Figure 6(c) intuitively show the existence of 1:2 strong resonance bifurcation.The critical values of bifurcation parameter d are large which illustrates the effect of step size on the dynamics.

Conclusions
In this paper, we discuss the existence and stability of the fixed points of system (1.4).The bifurcations of system (1.4) have also been investigated and numerical analysis has been used to support our work.We provide an analysis of the codimension 1 and 2 bifurcations at P 3 , including fold, 1:1 strong resonance and fold-flip bifurcations.It demonstrates the abundant dynamics of system (1.4).Further, we have analyzed the existence of 1:2 strong resonance bifurcation at P 5 .And, through examples with specific values, we found that selecting another two independent coefficients as bifurcation parameters does not change the existence of this bifurcation.In addition, the continuation curves are used to show the new dynamical behaviors, including the generalized period-doubling and bifurcations at P 5 .For the interior fixed point P 4 , the analysis of bifurcations is similar, so we have omitted the detailed descriptions in this work.The occurrences of these codimension 2 bifurcations imply that system (1.4) undergoes local codimension 1 bifurcations, such as fold, flip, homoclinic and Neimark-Sacker bifurcations.These complex phenomena may lead to the species not being able to coexist in a stable state.
To the best of our knowledge, the content of our analysis for system (1.4) has not been studied.Our work demonstrates that, compared to the continuous system, the dynamical properties of the discrete system are variable.In particular, the dynamics of discrete systems are affected by the choice of a large step size of the Euler method.The difference between continuous and discrete systems is attractive.Additionally, 1:1 strong resonance bifurcation corresponds to Bogdanov-Takens bifurcation in continuous systems.However, generalized period-doubling bifurcation, fold-flip bifurcation and 1:2 strong resonance bifurcation are unique to the discrete system and have many interesting properties.And, the discrete systems can reflect the interactions of species more realistically when the populations of species are small or the processes of birth and death occur at discrete times.Moreover, other codimension 2 bifurcations and the harvesting of species are meaningful topics.We will consider these in the future.

Use of AI tools declaration
The authors declare that they have not used artificial intelligence tools in the creation of this article. .
B. The coefficients in the proof of Theorem 7.

Figure 4 .
Figure 4. (a) The NS and LP curves with d = 0.8, k = 0.7 and a 1 = 0.1; (b) phase portrait with d = 0.8, k = 0.7, a 1 = 0.1, a 2 = 0.13 and b = 0.35; (c) time-series diagram of the prey x and predator y populations with the same specific parameters as Figure 4(b).