Stability and bifurcation analysis of a discrete Leslie predator-prey system via piecewise constant argument method

1.   Introduction

As a fundamental component of ecological investigation, the interaction between predator and prey populations in biological ecosystems has long been a subject of interest and analysis among biologists and mathematicians. Predator-prey systems involve the dynamic interaction of two populations, in which the survival of one species (the predators) is dependent on the consumption of another species (the prey). A number of researchers have conducted thorough investigations and analyses of diverse mathematical frameworks that explicate the interrelationships between species. The fundamental objective of these inquiries has been to examine the complex and varied dynamics that are intrinsic to these systems. The primary areas of concern include stability analysis, limit cycles, bifurcations, chaotic phenomena, and associated phenomena ^{[1,2,3,4,5,6,7,8]}.

The Lotka-Volterra system, which is the classic predator-prey framework, was initially proposed and established by Volterra and Lotka ^[9,10]. Within the conceptual framework of this system, it is postulated that all components adhere to linear functions, including the growth rate, predator mortality rate, and the transformation of prey biomass into predator reproduction. However, it is crucial to acknowledge that nonlinear components exert a substantial influence on predator-prey dynamics that occur in nature. Scholars have suggested numerous improvements to address the system's inability to incorporate specific real-world events. Leslie and Gower ^[11,12] introduced the classical Leslie-Gower predator-prey system to provide a comprehensive representation of the real-world scenario. This system incorporates the notion that the predator's environmental carrying capacity is directly proportional to the population density of the prey. Britton ^[13] suggested the following Leslie's system:

where $x$ and $y$ represent the population densities of prey and predator, respectively. The parameter $a$ denotes the intrinsic growth rate of prey, while $b$ quantifies the strength of competition among individuals of prey. The rate of change due to interaction is represented by the parameter $c$. Additionally, the intrinsic growth rate of the predator is denoted by $d$ and $\alpha$ signifies the density of prey required to sustain a single predator. All parameters $a, b, c, d$ and $\alpha$ are positive constants.

There are two distinct categories of modeling methodologies, namely, discrete-time systems and continuous-time systems. Much research has been dedicated to examining the nonlinear dynamical characteristics of continuous systems. In contrast, the investigation of discrete systems has been relatively limited. Nevertheless, over time scholars have come to recognize that discrete systems possess distinct dynamical qualities when compared to their continuous counterparts. The discrete-time systems provide the most accurate description of the dynamics exhibited by animals that engage in seasonal reproduction and have nonoverlapping generations. Furthermore, as compared to continuous-time systems, these discrete systems demonstrate much more complicated dynamical patterns ^{[14,15,16,17,18,19,20,21]}. Hence, discrete systems possess more attraction in comparison to continuous systems. It is crucial to use suitable numerical techniques for discretizing the continuous systems. The forward Euler strategy ^{[22,23,24,25,26]} has been widely used due to its simplicity. However, there needs to be more dynamic consistency when comparing the discrete version to its continuous equivalents. The discrete system produced by the use of the Euler approach lacks realism. Certain parameters and initial values may attribute the presence of negative values for the prey and predator population sizes. Nevertheless, using the piecewise constant argument approach ^{[27,28,29,30]} eliminates negative values. Khan et al. ^[31] investigated the discrete form of the system (1.1) as follows:

They obtained the discretization using the forward Euler technique. They investigated stability analysis and demonstrated the occurrence of the Neimark-Sacker (NS) bifurcation in the system. The authors used the hybrid control methodology to regulate the chaotic behavior in the system. In this study, we use the piecewise constant argument as an alternative method for discretizing the continuous-time system (1.1). By utilizing piecewise constant arguments to solve nonlinear differential equations and considering the regular time interval for the average growth rate in both populations, we can rewrite system (1.1) as follows:

where $[t]$ represents the integer part of $t$ and $0 < t < \infty$. Furthermore, integrating system (1.3) on an interval $[n, n+1)$ with $n = 0, 1, 2, \cdots$ yields the following system:

Taking $t \to n+1$, we obtain the following discrete-time system:

The subsequent sections of the paper are structured in the following manner: The objective of Section 2 is to investigate the presence and topological classification of fixed points. Section 3 provides a comprehensive examination of the bifurcation analysis pertaining to the period-doubling (PD) and NS bifurcations occurring at the positive fixed point. Section 4 employs two control methodologies in order to effectively regulate bifurcations and chaos. In order to validate and explicate the theoretical findings, numerical examples are presented in Section 5. The analysis performed in this research is ultimately summarized in Section 6.

2.   Topolgical classification of fixed points

The understanding of fixed point stability holds significant importance within the context of a predator-prey system. The fixed points represent equilibrium states wherein the populations of predators and prey have achieved a condition of balance. By conducting an analysis of their stability, we are able to make predictions about the long-term patterns exhibited by these ecological systems, gaining a deeper understanding of the various factors that contribute to the overall dynamics of the ecosystem.

2.1.   Existence of fixed points

The fixed points for system (1.5) can be obtained by solving

for $x$ and $y$. It is obtained that system (1.5) has two fixed points $E_1 = (\frac{a}{b}, 0), E_2 = (\frac{a \alpha}{cd+b \alpha}, \frac{ad}{cd+b \alpha})$.

2.2.   Stability of fixed points

To classify the fixed points, we employ the following results.

Lemma 2.1. ^[32] Let $\Lambda(\xi) = \xi^2+K_1 \xi+K_0$ be the characteristic polynomial of Jacobian matrix computed at fixed point $(x, y)$ and $\xi_1, \xi_2$ satisfy $\Lambda(\xi) = 0$, then $(x, y)$ is a

(1) sink (locally asymptotically stable (LAS)) when $|\xi_1| < 1$ along with $|\xi_2| < 1$,

(2) source when $|\xi_1| > 1$ along with $|\xi_2| > 1$,

(3) saddle point (SP) when $|\xi_1| < 1 \ {\mathit{\text{and}}} \ |\xi_2| > 1$ (or $|\xi_1| > 1 \ {\mathit{\text{and}}} \ |\xi_2| < 1$),

(4) non-hyperbolic point (NHP) when the modulus}of either of $\xi_1$ and $\xi_2$ is one.

Lemma 2.2. ^[32] Consider the quadratic function $\Lambda(\xi) = \xi^2+K_1\xi+K_0$. Suppose that $\Lambda(1) > 0$. If $\xi_1$ and $\xi_2$ both satisfy the equation $\Lambda(\xi) = 0$, then

(1) $|\xi_1| < 1$ along with $|\xi_2| < 1$ if $\Lambda(-1) > 0 \ {\mathit{\text{and}}} \ K_0 < 1$,

(2) $|\xi_1| < 1 \ {\mathit{\text{and}}} \ |\xi_2| > 1$ (or $|\xi_1| > 1 \ {\mathit{\text{and}}} \ |\xi_2| < 1$) if $\Lambda(-1) < 0$,

(3) $|\xi_{1, 2}| > 1$ if $\Lambda(-1) > 0 \ {\mathit{\text{and}}} \ K_0 > 1$,

(4) $|\xi_2|\neq 1 \ {\mathit{\text{and}}} \ \xi_1 = -1$ if $\Lambda(-1) = 0 \ {\mathit{\text{and}}} \ K_1 \neq 0, 2$,

(5) $\xi_1, \xi_2 \in \mathbb{C}$ along with $|\xi_{1, 2}| = 1$ if $K_1^2-4K_0 < 0 \ {\mathit{\text{and}}} \ K_0 = 1$.

Through simple computations, one can obtain the Jacobian matrix at an arbitrary fixed point $(x, y)$ as follows:

The Jacobian matrix of system (1.5) computed at $E_1$ is shown as follows:

Also, the Jacobian matrix at $E_2$ is given by

Proposition 2.1. The boundary fixed point $E_1$ is

(1) never LAS,

(2) a source if $a > 2$,

(3) SP if $0 < a < 2$,

(4) NHP if $a = 2$.

Next, we classify the positive fixed point $E_2$ of system (1.5) according to the above Jacobian matrix and Lemma 2.2. We obtain

The characteristic polynomial of $J(E_2)$ is

It can be obtained through calculations that

Theorem 2.1. The following holds for the interior fixed point $E_2$ of the system (1.5):

(1) $E_2$ is a sink if one of the following conditions holds:

(a) $d < 2$ and if one of the following conditions holds:

(i) $b+\frac{c(-4+d)d^2}{(-2+d)^2\alpha}\leq 0$ and $a < \frac{d(cd+b \alpha)}{cd^2+b(-1+d)\alpha}$,

(ii) $b+\frac{c(-4+d)d^2}{(-2+d)^2\alpha} > 0$ and $a < \frac{2(-2+d)(cd+b\alpha)}{cd^2+(-2+d)\alpha}$,

(b) $d = 2$ and $a < \frac{4c+2b\alpha}{4c+b\alpha}$,

(c) $2 < d < 4, b+\frac{c(-4+d)d^2}{(-2+d)^2\alpha} < 0$, and

(2) $E_2$ is a saddle point if one of the following conditions holds:

(a) $0 < d < 2, b > \frac{cd^2}{(2-d)\alpha}$, and $a > \frac{2(-2+d)(cd+b\alpha)}{cd^2+b(-2+d)\alpha}$,

(b) $d > 2$ and $0 < a < \frac{2(-2+d)(cd+b\alpha)}{cd^2+b(-2+d)\alpha}$,

(3) $E_2$ is a source if one of the following conditions holds:

(a) $d = 2$ and $a > \frac{4c+2b\alpha}{4c+b\alpha}$,

(b) $d > 4$ and $a > \frac{2(-2+d)(cd+b\alpha)}{cd^2+b(-2+d)\alpha}$,

(c) $a > \frac{d(cd+b \alpha)}{cd^2+b(-1+d)\alpha}$ and if one of the following conditions holds:

(i) $0 < d < 2$ and $b \geq \frac{cd^2}{(2-d)\alpha}$,

(ii) $2 < d\leq 4$ and $b+\frac{c(-4+d)d^2}{(-2+d)^2\alpha}\leq 0$,

(d) $d < 2, \frac{cd^2}{(2-d)\alpha} < b < -\frac{c(-4+d)d^2}{(-2+d)^2\alpha}$ and

(e) $2 < d\leq 4, b+\frac{c(-4+d)d^2}{(-2+d)^2\alpha} > 0$, and $a > \frac{2(-2+d)(cd+b\alpha)}{cd^2+b(-2+d)\alpha}$,

(4) $E_2$ is non-hyperbolic point and experiences PD bifurcation if $a = A_1 = \frac{2(-2+d)(cd+b\alpha)}{cd^2+b(-2+d)\alpha}$ and if one of the requirements listed below is satisfied:

(a) $b+\frac{c(-4+d)d^2}{(-2+d)^2\alpha} < 0$ and if one of the following conditions holds:

(i) $2 < d\leq 4$,

(ii) $d < 2$ and $b > \frac{cd^2}{(2-d)\alpha}$,

(b) $d \in (0, 2)\cup(2, 4]$ and $b+\frac{c(-4+d)d^2}{(-2+d)^2\alpha} > 0$,

(c) $d > 4$,

(5) $E_2$ is non-hyperbolic point and experiences NS bifurcation if $a = A_2 = \frac{d(cd+b \alpha)}{cd^2+b(-1+d)\alpha}$ and if one of the requirements listed below is satisfied:

(a) $d = 2$,

(b) $d \in (0, 2)\cup(2, 4)$ and $b+\frac{c(-4+d)d^2}{(-2+d)^2\alpha} < 0$.

3.   Bifurcation analysis

This section is focused on conducting a comprehensive investigation of the bifurcation phenomenon involving PD and NS bifurcation in the system (1.5) at the positive fixed point $E_2$. In order to obtain a thorough treatment of bifurcation analysis, we recommend that readers refer to ^[33,34]. The bifurcation holds significant implications for the dynamics of the system, shedding light on scenarios where even slight modifications to parameters yield substantial alterations in the dynamics of predator-prey relationships. In addition, gaining knowledge about PD and NS bifurcations contributes to a deeper comprehension of ecosystem dynamics. This understanding, in turn, facilitates the formulation of effective conservation and management strategies aimed at sustaining the enduring coexistence of predator and prey populations. This study initiates by investigating the PD bifurcation at $E_2$ based on condition (4-c) as presented in Theorem 2.1. By applying a small perturbation $\delta$ ($|\delta|\lll 1$) to the bifurcation parameter around the critical value $A_1$, the system (1.5) is changed to

We transform the fixed point $E_2$ to the origin by considering the change of variables $u_n = x_n-\frac{(a+\delta) \alpha}{cd+b \alpha}, \ v_n = y_n-\frac{(a+\delta)d}{cd+b \alpha}$. As a result, the system (3.1) is transformed into the following form:

where

Next, the system (3.2) is diagonalized through the consideration of the following transformation:

Upon applying the mapping (3.3), the system (3.2) undergoes the alteration as follows:

where

where the values of coefficients are given in Appendix A. Next, assume that $Q^C$ is the center manifold of (3.4) intended at origin in a close neighborhood of $\delta = 0$. It can be approximated as follows:

where

As a result, the system (3.4) is limited to $Q^C$ in the manner as follows:

For the function (3.5) to go through PD bifurcation, it is necessary that the following two quantities possesses nonzero values:

From the previous discussion, we get the following theorem:

Theorem 3.1. Assume that condition (4-c) of Theorem 2.1 holds true. The system (1.5) experiences PD bifurcation at $E_2$ if $l_1, l_2$ given in (3.6) and (3.7) are nonzero and $a$ fluctuates in a close neighborhood of $A_1 = \frac{2(-2+d)(cd+b\alpha)}{cd^2+b(-2+d)\alpha}$. Moreover, if $\ l_2 > 0$ (respectively, $l_2 < 0$), then a period-2 orbit of the system (1.5) emerges and it is stable (respectively, unstable).

Next, we proceed to investigate the NS bifurcation at $E_2$ under condition (5-a) stated in Theorem 2.1. By applying a small perturbation $\delta$ ($|\delta|\lll 1$) to the bifurcation parameter around the critical value $A_2$, the system (1.5) is changed to

We transform the fixed point $E_2$ to the origin by considering the change of variables $u_n = x_n-\frac{(a+\delta) \alpha}{cd+b \alpha}, \ v_n = y_n-\frac{(a+\delta)d}{cd+b \alpha}$. As a result, the system (3.8) is transformed into the following form:

where

The characteristic equation of the Jacobian matrix of system (3.9) estimated at origin is

where

The complex solutions for (3.10) are calculated as:

Moreover, we obtain

Additionally, it is required that $\xi_{1, 2}^k\neq 1$ when $\delta$ = $0$ for $k = 1, 2, 3, 4$, which corresponds to $\alpha(0)$ $\neq$ $-2, 2, 0, 1$. We obtain

Moreover, $\alpha(0) = -2$ is equivalent to $c = 0$, which is not possible. Next, to change (3.9) into normal form at $\delta = 0$, we use the following transformation:

Upon application of the mapping (3.12), the system (3.9) takes the following form:

where

and

The map (3.13) can undergo NS bifurcation if the following quantity is nonzero:

where

Therefore, the result derived from the above analysis is as follows:

Theorem 3.2. Suppose that condition (5-a) of Theorem 2.1 holds true. If $L$ given in (3.14) holds a nonzero value, then system (1.5) experiences NS bifurcation at $E_2$ as long as $a$ fluctuates in a close neighborhood of $A_2 = \frac{d(cd+b \alpha)}{cd^2+b(-1+d)\alpha}$. Furthermore, in instances where $L$ is negative (alternatively, positive), the NS bifurcation encountered in system (1.5) at $E_2$ is categorized as supercritical (subcritical), giving rise to the presence of a unique closed invariant curve originating from $E_2$, that is, attracting (repelling).

4.   Chaos control

The aim of optimizing dynamical systems in order to meet particular performance criteria and minimize chaotic behavior is a highly desirable objective. Chaos control techniques are extensively employed in several fields of applied research and engineering. Historically, bifurcations and unstable oscillations have been regarded in a negative light within the field of mathematical biology due to their detrimental impact on the reproductive capacity of biological populations. It is possible to design a controller that may modify the bifurcation features of a nonlinear system in order to achieve specific desired dynamical attributes and effectively control chaos under the effects of PD and NS bifurcations. Multiple strategies exist for the purpose of managing chaos in a discrete-time system. This section is dedicated to examining two distinct control methods, namely, state feedback control and hybrid control approaches.

Initially, the state feedback control strategy, as described in references ^[35,36], is employed to regulate the chaotic behavior of the system (1.5). The suggested methodology entails the conversion of the chaotic system into a piecewise linear system in order to obtain an optimal controller that effectively reduces the upper limit. Following this, the problem of optimization is carried out subject to certain constraints. The technique described above is utilized in order to attain stabilization of chaotic orbits situated at an unstable fixed point within the system (1.5). The controlled system under consideration for this purpose is as follows:

where $U_n = \kappa_1 \bigg(x_n-\frac{a \alpha}{cd+b \alpha}\bigg)+\kappa_2 \bigg(y_n-\frac{ad}{cd+b \alpha}\bigg)$ is the feedback controlling force with feedback gains $\kappa_1$ and $\kappa_2$. Through simple calculations, it is obtained that for system (4.1) we have

The matrix $J(E_2)$ has the following characteristic equation:

where

Let $\xi_1$ and $\xi_2$ be the roots of (4.3), then we have

The lines of marginal stability are derived by solving $\xi_1 = \pm 1$ and $\xi_1 \xi_2 = 1$. These conditions ensure that $|\xi_{1, 2}| < 1$. Assume that $\xi_1 \xi_2 = 1$, then Eq (4.5) implies that

Next, we take $\xi_1 = 1$ and, utilizing Eqs (4.4) and (4.5), we obtain

Next, we take $\xi_1 = -1$ and, utilizing Eqs (4.4) and (4.5), we obtain

The stable eigenvalues are enclosed within the triangular region bounded $L_1, L_2$ and $L_3$.

Next, we apply the hybrid control approach ^[37] for controlling chaos through both types of bifurcation effects. The hybrid control technique refers to a methodology that integrates the utilization of state feedback and parameter adjustment in order to achieve stabilization of unstable periodic orbits that are present within the chaotic attractor of a given system. The controlled system of (1.5), when the hybrid control approach is used, becomes

where $\rho \in (0, 1)$. The parameter $\rho$ serves as a control parameter, which balances the influence of the original system (1.5) in comparison to the modified system (4.9). If the value of $\rho$ is negative, it might suggest the reverse impact of the original system (1.5). If the value of $\rho$ is more than $1$, it suggests that the original system (1.5) has an intensifying impact on the modified system (4.9), perhaps resulting in impractical or unfeasible outcomes. Systems (4.9) and (1.5) share the same fixed points. We obtain

then, its characteristic equation is as follows:

where

Theorem 4.1. The fixed point $E_2$ of the controlled system (4.9) is LAS if

5.   Numerical examples

In this section, we provide empirical evidence to support our theoretical conclusions for system (1.5) through the implementation of numerical simulations. The numerical simulations will include the depiction of bifurcation diagrams, phase portraits, time series plots, and graphs illustrating the maximum Lyapunov exponent (MLE).

Example 5.1. We assume that $b = 0.5, c = 0.8, d = 2, \alpha = 2.5$, then system (1.5) experiences NS bifurcation as $a$ varies in small neighborhoods of $a_2 \approx 1.280899$. The positive fixed point is obtained as $E_2 = (1.123596, 0.898876)$ for $a = 1.280899$. The eigenvalues of $J(E_2)$ are $\xi_{1, 2}$ = $-0.280899$$\pm 0.959737i$ with $|\xi_{1, 2}| = 1$. Moreover, some careful calculations give

Thus, it is obtained that $L = -0.257431 < 0$, which proves the correctness of Theorem 3.2. The bifurcation diagrams of system (1.5) are given in Figure 1(a) and (b), while the MLE is plotted in Figure 1(c) by using initial conditions $x_0 = 1.10$ and $y_0 = 0.90$ and varying $a \in [1.05, 2.25]$.

Next, Figure 2(a)–(f) show phase portraits of system (1.5) for various values of $a$. The positive fixed point $E_2$ remains stable when $a < 1.280899$. However, owing to an NS bifurcation, it loses stability at $a = 1.280899$. This results in the formation of an invariant closed curve, the radius of which expands as $a$ grows. Furthermore, the MLE graph demonstrates that negative values of the MLE indicate a range of periodic orbits. For instance, orbit-7 in Figure 2(d), orbit-11 in Figure 2(f), orbit-15 in Figure 2(h), orbit-4 in Figure 2(l) and orbit-8 in Figure 2(m) exhibit this behavior. Conversely, positive MLE values suggest that the system is chaotic. Figure 2(o) and (p) exhibit the presence of strange attractors, indicating the chaotic character that results from positive MLE values.

Example 5.2. We assume that $b = 0.5, c = 0.8, d = 3.15, \alpha = 2.5$, then system (1.5) experiences both PD bifurcation and NS bifurcation as $a$ varies in small neighborhoods of $A_1 \approx 0.924857$ and $A_2 \approx 1.11764$, respectively. The positive fixed point is obtained as $E_2 = (0.613301, 0.772759)$ for $a = A_1$. The eigenvalues of $J(E_2)$ are $\xi_1 = -1$ and $\xi_2 = -0.45665$ with $|\xi_2| \neq 1$. The positive fixed point is obtained as $E_2 = (0.741142, 0.933838)$ for $a = A_2$. The eigenvalues of $J(E_2)$ are $\xi_{1, 2}$ = $-0.760285$$\pm 0.649589i$ with $|\xi_{1, 2}| = 1$. The bifurcation diagrams of system (1.5) are given in Figure 3(a) and (b), while the MLE is plotted in Figure 3(c) by using initial conditions $x_0 = 0.60$ and $y_0 = 0.75$ and varying $a$.

Next, Figure 4(a)–(l) show phase portraits of system (1.5) for various values of $a$. One can observe that $E_2$ is LAS for $0.924857 < a < 1.11764$, but loses stability at $a = 0.924857$ when the system (1.5) experiences PD bifurcation. For $a \geq 1.11764$, the system (1.5) experiences NS bifurcation and an invariant curve emerges from $E_2$, the radius of which grows as $a$ grows. It leads to a strange chaotic attractor presented in Figure 4(l).

Therefore, it can be concluded that the positive fixed point $E_2$ of the system (1.5) exhibits stability when the growth rate of the prey population falls within the range $(A_1, A_2)$. When the growth rate of the prey falls below the threshold value of $A_1$, it experiences difficulties in maintaining its population against predation, leading to a reduction that subsequently affects the survival of predators that rely on it. When the growth rate exceeds the value of $A_2$, rapid growth sustains predators initially but carries the potential for resource exhaustion, therefore affecting both prey and predators over a long period of time.

Subsequently, our objective is to assess the effectiveness of the feedback control strategy. Given the values $a = 1.15, b = 0.5, c = 0.8, d = 3.15,$ and $\alpha = 2.5$, together with the initial conditions $x_0 = 0.60$ and $y_0 = 0.75$ for the controlled system (4.1), the marginal stability lines may be determined as follows:

and

The stability region of system (4.1) is shown in Figure 5(a), which is bounded by lines $L_1, L_2$ and $L_3$. The instability of the fixed point $E_2$ of system (1.5) has been shown for the specified parametric values. The system described by Eq (4.1) is analyzed using feedback gains $\kappa_1 = -3.01$ and $\kappa_2 = 1.60$. Figure 5 illustrates the graph of $x_n$ as shown in Figure 5(c), $y_n$ as shown in Figure 5(d), and the phase portrait as presented in Figure 5(b) for the system (4.1). Hence, it may be inferred that the use of the feedback control approach seems to be efficacious in managing bifurcation and chaos.

The effectiveness of the hybrid control strategy will now be assessed. We assume $\rho = 0.95$, $b$ = $0.5$, $c = 0.8, d = 3.15, \alpha = 2.5$ and vary $a$ for the controlled system (4.9). If $0.897012 < a < 1.183794$, the positive fixed point $E_2$ is LAS. One can observe that the stability region has been expanded. The bifurcation diagrams are presented in Figure 6(a) and (b) by using initial conditions $x_0 = 0.60$, $y_0$ = $0.75$ and varying $a \in [0.20, 1.20]$.

6.   Conclusions

Our study focused on a discrete-time system, which provided a more accurate depiction of animal dynamics in seasonal reproduction with nonoverlapping generations. These discrete systems, in comparison to continuous-time systems, displayed more detailed and rich dynamical patterns, making them more attractive. When discretizing continuous systems, it was critical to use proper numerical approaches. In this paper, we used the piecewise constant argument methodology to discretize the system (1.1), leading to a discrete system (1.5) that consistently preserved positive solutions. This approach presented a benefit in comparison to the Euler method, which was susceptible to producing negative solutions, thereby posing challenges for population systems. This study conducted a comprehensive analysis of the existence and stability of fixed points, focusing specifically on local bifurcations that occur at the positive fixed point. It was observed that our system (1.5) undergoes bifurcations in both PD and NS at a positive fixed point, whereas prior investigations ^[31] only reported NS bifurcations. Hybrid control and feedback control techniques were employed to effectively manage bifurcation and disorder within the system denoted as (1.5). In order to validate our theoretical findings, we employed numerical simulations that incorporated a range of graphical representations, such as time series plots, phase portraits, bifurcation diagrams, and MLE graphs.

It can be inferred that the stability of the positive fixed point $E_2$ is maintained when the prey population growth rate lies within the interval $(A_1, A_2)$. In situations where the prey's growth rate fell below a certain threshold value denoted as $A_1,$ it encountered challenges in sustaining its population in the face of predation, resulting in a decline that ultimately impacts the viability of predators dependent on it. When the growth rate exceeded the value of $A_2$, rapid growth sustained predators initially but carried the potential for resource exhaustion, therefore affecting both prey and predators over a long period of time.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This study is supported via funding from Prince Sattam bin Abdulaziz University project number (PSAU/2023/R/1445).

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix A