Bifurcation analysis and chaos control of a discrete fractional-order Leslie-Gower model with fear factor

1.   Introduction

The well-known Leslie-Gower model with a Holling Ⅱ functional response is expressed as follows:

where the biological interpretation of each parameter is provided in Table 1. It is noteworthy that the Leslie-Gower model incorporates a response function indicating the predator's carrying capacity, which is proportional to the prey population, a feature absent in the Lotka-Volterra model. This distinction has been extensively investigated in many references, such as ^[1,2,3,4,5].

Previous studies predominantly focused on direct predation effects. In 2016, Wang et al. ^[6] demonstrated through experiments that predator-induced fear (indirect effects) in prey leads to a reduction in prey birth rates. They introduced the fear factor $F(k, y) = \frac{1}{1+ky}$, where $k$ represents the intensity of fear driving antipredator behavior. To further investigate the influence of fear on population dynamics, we enhance the system (1.1) by incorporating the fear factor $F(k, y)$. The modified model is formulated as follows:

For additional significant findings on the impact of fear in predator-prey models, refer to ^[7,8,9].

There has been increasing recognition that traditional integer-order differential equations may not sufficiently capture the complexity of biological systems ^{[10,11,12,13,14]}. In reality, the behaviors of most organisms in nature are influenced by their historical context. Fractional-order derivatives, which extend the concept of integer-order differentiation, offer a more flexible and accurate framework for modeling memory and hereditary properties in population dynamics. Therefore, the authors of this paper aim to apply the Caputo fractional derivative to (1.2), thereby extending it into a fractional model. The Caputo fractional derivative of a function $u(t)$ of order $\alpha\in(0, 1]$ is given in ^[15] as follows

where $\Gamma$ is the Gamma function. By replacing the integer-order derivative with the Caputo fractional derivative, the following model is obtained:

Specifically, when $\alpha = 1$, the fractional-order model (1.3) reduces to the integer-order model (1.2), demonstrating that the fractional-order model serves as a generalization of the integer-order model. Moreover, in fractional calculus, the rate of change at any given moment, expressed by the fractional-order derivative, depends on the population density over a specified time interval. This feature gives the fractional-order model (1.3) a distinct advantage in capturing memory effects within populations.

When studying populations with nonoverlapping generations or small sizes, mathematical models are often expressed in discrete terms. Although some variables may change continuously in real-world scenarios, the recording of these changes typically happens at specific time intervals during data collection. Thus, employing discrete systems to examine the dynamic behavior of biological populations is highly practical and significant. In references ^[16,17], the authors employed the piecewise constant approximation method to discretize continuous fractional predator-prey models and explored their dynamical properties. Singh and Sharma ^[18] examine a discrete prey-predator model with Holling type Ⅱ functional response and prey refuge, identifying bifurcations and controlling chaos through state feedback, pole placement, and hybrid techniques. Berkal and Almatrafi ^[19] used the exponential piecewise constant argument to discretize continuous fractional activator-inhibitor system. Their analysis includes stability assessments, investigations of Neimark-Sacker and period-doubling bifurcations, and numerical simulations that validate the theoretical findings on the system's dynamics. For a more comprehensive exploration of discrete model studies, readers are advised to consult references ^{[20,21,22,23,24,25]} and the related literature cited therein. However, current literature lacks studies on fractional discrete-time predator-prey Leslie-Gower systems that incorporate the fear effect in prey population. Using the same method as in references ^[16,17] to discretize model (1.3), we can obtain the following discrete model:

Here $h > 0$ represents the time interval of production.

Furthermore, the key contributions and findings of this study are summarized as follows:

● The Caputo fractional derivative of order $(0, 1]$ is utilized to incorporate the memory effect into the dynamical behavior of the proposed model.

● The existence and stability of fixed points are investigated.

● Conditions for the occurrence and direction of period-doubling bifurcation and Neimark–Sacker bifurcation at the positive fixed point are established.

● State feedback and hybrid control strategies are employed to manage bifurcations and chaotic behavior in the model.

● To validate the accuracy of our theoretical findings, numerical examples for the fractional-order discrete-time Leslie-Gower model with a fear factor are provided.

The remainder of this paper is structured as follows: In Section 2, we investigate the existence and stability of the equilibrium points of model (1.4). Section 3 analytically demonstrates that model (1.4), under specific parametric conditions, undergoes period-doubling or Neimark-Sacker bifurcation. Section 4 explores the control of chaos toward an unstable equilibrium point using feedback control or hybrid control approaches. Section 5 presents a quantitative analysis of the dynamics of model (1.4) to validate our analytical findings. Finally, Section 6 offers brief conclusions.

2.   Existence and stability of the equilibrium points

This section explores the existence of equilibrium points in model (1.4) and assesses their stability by evaluating the eigenvalues of the Jacobian matrix at these points. The following definition and lemma are introduced to assist in this stability analysis of equilibrium points.

Definition 1. ^[26] Let $\lambda_{1}$ and $\lambda_{2}$ denote the two roots of the characteristic equation $F(\lambda) = \lambda^2+p\lambda+q = 0$ associated with the Jacobian matrix $J(x, y)$. The equilibrium point $(x, y)$ is termed

(1) sink if $|\lambda_{1}| < 1$ and $|\lambda_{2}| < 1$, and the sink is locally asymptotically stable;

(2) source if $|\lambda_{1}| > 1$ and $|\lambda_{2}| > 1$, and the source is locally unstable;

(3) saddle if $|\lambda_{1}| > 1$ and $|\lambda_{2}| < 1$ (or $|\lambda_{1}| < 1$ and $|\lambda_{2}| > 1$);

(4) non-hyperbolic if either $|\lambda_{1}| = 1$ or $|\lambda_{2}| = 1$.

Lemma 1. ^[26] Let $F(1) > 0$ in $F(\lambda) = \lambda^2 - M\lambda + N$, where $\lambda_{1}, \lambda_{2}$ are the two roots of $F(\lambda) = 0$. Then, the following results hold true:

(1) $|\lambda_{1}| < 1$ and $|\lambda_{2}| < 1$ if, and only if, $F(-1) > 0$ and $N < 1$;

(2) $|\lambda_{1}| < 1$ and $|\lambda_{2}| > 1$ (or $|\lambda_{1}| > 1$ and $|\lambda_{2}| < 1$) if, and only if, $F(-1) < 0$;

(3) $|\lambda_{1}| > 1$ and $|\lambda_{2}| > 1$ if, and only if, $F(-1) > 0$ and $N > 1$;

(4) $\lambda_{1} = -1$ and $|\lambda_{2}| \neq 1$ if, and only if, $F(-1) = 0$ and $M \neq 0, -2$;

(5) $\lambda_{1}$ and $\lambda_{2}$ are conjugate complex and $|\lambda_{1, 2}| = 1$ if, and only if, $M^2 < 4N$ and $N = 1$.

2.1.   Existence of equilibrium points

It is evident that the equilibrium points of model (1.4) satisfy the following equations:

This algebraic system is satisfied if $x = K$ and $y = 0$, indicating that the model (1.4) has a boundary equilibrium point $E_0(K, 0)$ for all model parameters. To find the interior equilibrium point $E_{1}$, we will solve the following system simultaneously:

From the second equation, we find $y = x$. Substituting $y = x$ into the first equation yields:

Given that all model parameters are positive, we obtain the following results through direct calculations.

Theorem 1. The model (1.4) always has a boundary equilibrium point $E_0(K, 0)$ and a positive equilibrium point $E_{1} = (x^*, y^*)$, where

2.2.   Stability of equilibrium points

The Jacobian matrix of model (1.4) evaluated at any point $(x, y)$ is as follows:

Theorem 2. The boundary equilibrium point $E_0(K, 0)$ exhibits the following behaviors:

(1) It is a source point if $r > \frac{2\Gamma(\alpha+1)}{h^{\alpha}}$;

(2) It is a saddle point if $0 < r < \frac{2\Gamma(\alpha+1)}{h^{\alpha}}$;

(3) It is a non-hyperbolic point if $r = \frac{2\Gamma(\alpha+1)}{h^{\alpha}}$.

Proof. The Jacobian matrix (2.1) evaluated at the boundary equilibrium point $E_{0}(K, 0)$ is given by:

The eigenvalues of $J(E_{0})$ are $\lambda_{1} = 1 - \frac{h^{\alpha} r}{\Gamma(\alpha+1)}$ and $\lambda_{2} = 1 + \frac{h^{\alpha} s}{\Gamma(\alpha+1)}$. Clearly, $|\lambda_{2}| > 1$ and:

By applying Definition 1, the proof is complete. □

Next, we analyze the local dynamics of model (1.4) at the positive equilibrium point $E_1(x^*, y^*)$. The Jacobian matrix (2.1) evaluated at $E_1(x^*, y^*)$ is expressed as:

where

The characteristic equation of $J(E_{1})$ is given by:

where

Let $F(\lambda) = \lambda^{2} - M \lambda + N$, then:

As a result, we establish the following theorem.

Theorem 3. Let $E_{1}$ be the unique positive equilibrium point of model (1.4). Then:

(1) $E_{1}$ is a sink point if $|M| < 1 + N < 2$;

(2) $E_{1}$ is a source point if $|M| < |1 + N|$ and $|N| > 1$;

(3) $E_{1}$ is a saddle point if $M^2 > 4N$ and $|M| > |1 + N|$;

(4) $E_{1}$ is non-hyperbolic if $|M| = |1 + N|$ or $N = 1$ and $|M| < 2$.

Proof. (1) According to Lemma 1, $E_{1}$ is a sink point if and only if $F(1) > 0$, $F(-1) > 0$, and $N < 1$. This condition is satisfied when $|M| < 1 + N < 2$. Similarly, the proofs for Theorem 3 (2)–(4) follow straightforwardly from the definitions and properties of $M$ and $N$. □

3.   Bifurcation analysis

In this section, we will analyze the period-doubling and Neimark-Sacker bifurcation behaviors of model (1.4) at the positive equilibrium point $E_{1}(x^*, y^*)$. The reason for not analyzing the boundary equilibrium point is that, at $E_{0}(K, 0)$, the predators have become extinct, leaving only the prey.

3.1.   Period-doubling bifurcation

Assume $M^2 > 4N$ and

and we can ascertain that the characteristic equation (2.3) satisfies

Therefore, the eigenvalues of the characteristic equation (2.3) are

Let

The condition $|\lambda_{2}|\neq1$ implies that $s_{1}\neq s_{2}, s_{3}$.

Based on the above analysis, we can conclude that when the parameters vary within the set

model (1.4) will undergo period-doubling bifurcation at $E_{1}(x^*, y^*)$.

Select parameters $(h, \alpha, r, k, K, b, a, s) \in P.D$ and consider $s^*$ as a small perturbation of $s$, i.e., $s^* = s - s_1$, where $|s^*|\ll 1$. After this perturbation, model (1.4) can be represented as follows:

Let $u_n = x_n - x^*$ and $v_n = y_n - y^*$, transforming the equilibrium point $E_1(x^*, y^*)$ into the origin $O(0, 0)$. The model (3.2) can thus be expressed as:

The Taylor expansion of model (3.3) around $(u_n, v_n) = (0, 0)$ yields the following form:

where

and

Define

then

Using the transformation:

model (3.4) is transformed into the following form:

where

and

Next, we apply the center manifold theorem ^[26] to analyze the dynamics around the equilibrium point $(\widetilde{u}_n, \widetilde{v}_n) = (0, 0)$ at $s^* = 0$. According to the theorem, the model (3.5) has a center manifold, which can be represented as:

Assume that

Then, the center manifold must satisfy

By comparing the coefficients, it can be obtained that

Thus, the model (3.5), when restricted to the center manifold $W^{c}(0, 0, 0)$, is given by:

where

In order for Eq (3.6) to undergo period-doubling bifurcation, it is necessary that the following two quantities possess nonzero values.

Summarize the above analysis into the following theorem.

Theorem 4. If $\beta_{1}\beta_{2}\neq0$, then model (1.4) undergoes a period-doubling bifurcation at the positive equilibrium point $E_1(x^*, y^*)$ when parameters vary in a small neighborhood of $P.D$. Additionally, when $\beta_{2} > 0$ (respectively, $\beta_{2} < 0$), model (1.4) bifurcates from the equilibrium point $E_1(x^*, y^*)$ to a stable (respectively, unstable) 2-periodic orbit.

3.2.   Nerimark-Sacker bifurcation

Assume $M^2 < 4N$ and

and we can determine that the eigenvalues of the characteristic equation (2.3) satisfy

Namely, Eq (2.3) has two complex conjugate roots with unit modulus. It is clear from the above discussion that

In addition, it is crucial that when $s = s_{4}$, $\lambda^{\theta}_{1, 2}(s_{4})\neq 1 (\theta = 1, 2, 3, 4)$, which is equivalent to $M(s_{4})\neq -2, -1, 0, 2$. Since $M^2 < 4N$, we deduce $M(s_{4})\neq -2, 2$. Additionally, we necessitate that $M(s_{4})\neq 0, -1$, which leads to

Based on the preceding analysis, we conclude that when the parameters vary within a small neighborhood of the set

model (1.4) undergoes a Neimark-Sacker bifurcation at $E_{1}(x^*, y^*)$.

Next, assuming $(h, \alpha, r, k, K, b, a, s)\in N.S.$ and $|s^*|\ll1$ represents a small perturbation of $s_{4}$, model (1.4) can be described as follows:

Let $u_n = x_n-x^*$ and $v_n = y_n-y^*$, then we have

where

and $c_{13}$, $c_{14}$, $c_{15}$, $c_{16}$, $c_{17}$, $c_{18}$, $c_{19}$, $c_{23}$, $c_{24}$, $c_{25}$, $c_{26}$, $c_{27}$, $c_{28}$, $c_{29}$ are given in (2.3) by substituting $s_{1}$ for $s_{4}+s^*$.

In order to obtain the normal form of model (3.9) at $s^* = 0$, we use the following transformation:

with

Using this transformation, model (3.9) will transform as follows:

where

and

Next, a nonzero real number is defined as follows:

where

Through some complicated calculations, we get

and

Based on the aforementioned calculations, we establish the theorem regarding the existence and direction of Neimark-Sacker bifurcation.

Theorem 5. If $(h, \alpha, r, k, K, b, a, s)\in N.S.$ and $\Omega\neq0$, then model (1.4) undergoes a Neimark-Sacker bifurcation at the equilibrium point $E_{1}(x^*, y^*)$ when the parameter $s$ varies in the vicinity of $s_{4}$. Moreover, if $\Omega < 0$ (respectively, $\Omega > 0$), then an attracting (respectively, repelling) closed invariant curve bifurcates from the equilibrium point for $s^* > 0$ (respectively, $s^* < 0$).

4.   Chaos control

Chaos often has detrimental effects on biological systems, disrupting the ecological balance of populations and directly influencing long-term population growth projections. Implementing effective control policies not only safeguards the size of ecological populations but also establishes a strong foundation for the sustainable exploitation of ecological resources ^[27,28]. This section explores two control methods aimed at effectively managing the chaos generated by model (1.4).

4.1.   State feedback control

In this subsection, the state feedback control method ^[22] will be employed to regulate the chaos exhibited by model (1.4). To achieve this, we introduce the following controlled model.

which corresponds to model (1.4). The feedback controlling force is defined as

where $(x^{*}, y^{*})$ represents the positive equilibrium point of the model (1.4), and $p_{1}$, $p_{2}$ stand for the feedback gains. The Jacobian matrix of the controlled model (4.1) evaluated at the positive equilibrium point $E_{1}(x^{*}, y^{*})$ is given by

where the variables $A$, $a_{11}$, and $a_{12}$ are defined in Eq (2.2). The corresponding characteristic equation of the Jacobian matrix $J_{1}(x^{*}, y^{*})$ is

Let $\lambda_{1}$ and $\lambda_{2}$ be the roots of the Eq (4.3), then

The lines of marginal stability $l_{1}$, $l_{2}$, and $l_{3}$ are derived by solving $\lambda_{1}\lambda_{2} = 1$, $\lambda_{1} = 1$ and $\lambda_{2} = \pm1$, respectively. These conditions ensure that $|\lambda_{1, 2}| = 1$. Then, we derive the marginal stability lines as follows:

Therefore, $l_{1}$, $l_{2}$, and $l_{3}$ in the $(p_{1}, p_{2})$-plane form a triangular region which leads to $|\lambda_{1, 2}| < 1$.

Theorem 6. If $p_{1}$ and $p_{2}$ lie within a triangular region bounded by the lines $l_{1}$, $l_{2}$, and $l_{3}$, it can be concluded that the model (4.1) is stable.

4.2.   Hybrid control

Next, we apply the hybrid control approach proposed by ^[23] to control chaos. The controlled model of (1.4) with the hybrid control approach is depicted below:

where $0 < \rho < 1$, and the controlled strategy in (4.9) combines feedback control and parameter perturbation. By appropriately choosing the controlled parameter $\rho$, the chaotic behaviors of the equilibrium point $(x^{*}, y^{*})$ of the controlled model (4.9) can be accelerated (delayed) or even entirely eliminated. The Jacobian matrix of the controlled model (4.9), evaluated at the positive equilibrium point $(x^{*}, y^{*})$, is given by

where the variables $A$, $a_{11}$, and $a_{12}$ are defined in Eq (2.2). Then, the positive equilibrium point $(x^{*}, y^{*})$ of the controlled model (4.9) is locally asymptotically stable if the roots of the characteristic polynomial of (4.10) lie within the open unit disk. According to the Jury condition, the equilibrium point of the model remains stable if, and only if, the following conditions are met:

Theorem 7. If $|2+A\rho a_{11}-A\rho s| < 1+(1+A\rho a_{11})(1-A\rho s)-(A\rho a_{12})(A\rho s) < 2$ can hold, it can be concluded that the model (4.9) is stable.

5.   Numerical experiments

In this section, by exploring specific cases of model (1.4), we confirm the theoretical analysis above and discover new intriguing complex dynamic behaviors. Additionally, we validate the effectiveness of linear feedback techniques and hybrid control strategies for chaos control through numerical simulations.

5.1.   Period-doubling bifurcation at positive equilibrium point

Take

After straightforward calculations, we determine the positive equilibrium point $E_{1} = (1.7498, 1.7498)$ and the critical bifurcation value $s_{1} = 2.8472$. The Jacobian matrix of model (1.4) evaluated at $E_{1}$ with $s = s_{1}$ is given by

whose characteristic equation is

The roots of (5.2) are $\lambda_{1} = -1$ and $\lambda_{2} = 0.4768$. Moreover, we have $\beta_{1} = -0.7959 \neq 0$ and $\beta_{2} = 1.0080 > 0$. According to Theorem 4, model (1.4) undergoes period-doubling bifurcation at $E_{1}$ as $s$ passes through $s_{1}$. This behavior, verified through corresponding bifurcation diagrams shown in Figure 1, utilizes initial conditions $(x_0, y_0) = (1.74, 1.74)$ and varies $s$ in the range $[2.7, 3.7]$.

5.2.   Neimark-Sacker bifurcation at positive equilibrium point

Take

After straightforward calculations, we determine the positive equilibrium point $E_{1} = (2.1746, 2.1746)$ and the critical bifurcation value $s_{4} = 1.1872$. The Jacobian matrix of model (1.4) evaluated at $E_{1}$ with $s = s_{4}$ is given by

whose characteristic equation is

The roots of (5.4) are $\lambda_{1, 2} = 0.02930 \pm 0.9996i$. Moreover, we have $d = 0.3912 \neq 0$ and $\Omega = -1392.0464 < 0$. According to Theorem 5, model (1.4) undergoes Neimark-Sacker bifurcation at $E_{1}$ as $s$ passes through $s_{4}$. This behavior, verified through corresponding bifurcation diagrams shown in Figure 2, utilizes initial conditions $(x_0, y_0) = (2.2, 2.2)$ and varies $s$ in the range $[1.1, 1.6]$.

5.3.   Chaos control

We chose the parameter values as follows:

By Figure 3, we can get that the variables $x_{n}$ and $y_{n}$ in the model (1.4) are in a chaotic state.

5.3.1.   State feedback control

We use the same parameter values as in (5.5). In Figure 4, the triangular region defined by Theorem 6 bounds the parameters $p_{1}$ and $p_{2}$. Inside this region, the chaotic behavior produced by model (1.4) is effectively managed, resulting in asymptotic convergence toward the equilibrium point $E_{1} = (2.1746, 2.1746)$.

In this case, with feedback gains set to $p_{1} = 0.57$ and $p_{2} = -0.38$ and the controller activated at the 3000th iteration of the model, Figure 5 demonstrates the control effect. A chaotic trajectory is successfully stabilized at the equilibrium point $E_{1} = (2.1746, 2.1746)$. This indicates that the feedback control approach is effective in mitigating bifurcation and chaos.

5.3.2.   Hybrid control

Finally, using the same parameter values as in (5.5), the Jacobian matrix of the controlled model (4.9), evaluated at $E_{1}$, is given by:

The characteristic polynomial of (5.6) is given by

The roots of (5.7) lie within the open unit disk if, and only if, $0 < \rho < 0.881058$. Additionally, the plots for $x_{n}$ and $y_{n}$ of the controlled model (4.9) are shown in Figure 6 with $\rho = 0.87$. From Figure 6, it is clear that the positive equilibrium point $E_{1}$ is stable. Therefore, it can be concluded that employing the hybrid control approach is effective in mitigating bifurcation and chaos.

6.   Conclusions

In this study, we proposed a fractional Leslie-Gower model with a Holling type Ⅱ functional response and antipredator behavior. By employing the piecewise constant approximation method, we derived the discrete model (1.4) and analyzed its dynamical behavior, including the existence and stability of equilibrium points and the possibility of local bifurcations. The following conclusions can be drawn from our research:

(1) The discrete model (1.4) has two equilibrium points: $E_{0}(K, 0)$ and $E_{1}(x^*, y^*)$. The only positive coexistence equilibrium point, $E_{1}(x^*, y^*)$, reflects the coexistence of predators and prey.

(2) Our theoretical analysis and numerical simulations of the positive equilibrium point $E_{1}(x^*, y^*)$ indicate that the model undergoes period-doubling bifurcation and Neimark-Sacker bifurcation under specific parameter conditions. Figures 1 and 2 illustrate how these bifurcations can lead to chaotic behavior at $E_{1}$.

(3) By applying state feedback control and hybrid control methods, we effectively managed the chaotic behavior generated by the model (1.4), as shown in Figures 3–6. These interventions mitigated the adverse effects of chaos and bifurcations, consequently enhancing ecosystem resilience.

This study advances our understanding of the complex dynamics in ecological models influenced by fear effects and provides practical techniques for controlling chaotic behavior in such models. Future work could explore different discretization methods for the model, and new parameters could be chosen to study the influence of various ecological effects on population dynamics.

Author contributions

Yao Shi: Formal analysis, validation, writing-original draft & editing, visualization, methodology; Zhenyu Wang: Conceptualization, investigation, writing-review & editing, software, supervision. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

This work was supported by the Innovation Foundation of Hebei University of Engineering (Grant Nos. SJ2401002097), the National Natural Science Foundation of China (Grant Nos. 12401519), and the Central Guidance on Local Science and Technology Development Fund of Hebei Province (Grant Nos. 246Z1825G).

Conflict of interest

The authors declare no conflict of interest.