Abstract
This study examines the discrete prey-predator model in the sense of Caputo fractional derivative by incorporating harvesting on the predator population and immigration on the prey population. We establish the topological categories of the model’s fixed points. We show analytically that a fractional order prey-predator model supports both a Neimark–Sacker (NS) bifurcation and a period-doubling (PD) bifurcation under specific parametric circumstances. Using the central manifold and bifurcation theory, we provide evidence for NS and PD bifurcations. It has been discovered that the parameter values and the initial conditions have a significant influence on the dynamical behavior of the fractional order prey-predator model. Furthermore, two chaos management strategies are applied to eliminate the chaos that objectively exists in the model. Finally, numerical simulations are used to demonstrate complicated and chaotic behavior in order to support our theoretical and analytical discussions.
1. Introduction
In ecological science, the investigation of the predator-prey relationship with various ecological phenomena has gained a lot of popularity. An essential model in population dynamics, the dynamics of interacting populations [1] are examined using the prey-predator paradigm. Continuous-time population models, like the Lotka-Volterra model, have been utilized in population dynamics to comprehend the interaction between ecological species [2–8]. On the other hand, discrete-time population models have also gained attention recently [9–11] because these can create more complicated and interesting dynamical behaviors than continuous-time models and are better suited to model populations with non-overlapping generations. For example, a 1-dimensional discrete-time autonomous system can display chaos, but a continuous-time setup requires chaos in at least a 3-dimensional autonomous system [12, 13].
The traditional predator-prey relationship always takes the following shape:withwhere the time-dependent functions and stand in for the prey and predator population densities, respectively. Every constant is assumed to be positive. The carrying capacity is indicated by parameter . The predator’s mortality rate is represented by constant . represents the functional response, while represents the uptake functions.
Each population system in an ecological system employs a different strategy, such as refuging and grouping, to look for food sources and to defend itself. Therefore, a variety of ecological criteria and elements are used in the creation of mathematical models. The functional response in population dynamics, which refers to the amount of prey consumed by a predator dependent on the density of the prey per unit of time, is a crucial element in every prey-predator encounter. The functional response of the Holling type II [14] is commonly employed and well researched to those from the Holling types I, III, and IV for the majority of arthropod predators. To explore the dynamical interplay between prey and predator species, Iviev [15] devised a novel functional response, known as the Iviev functional response:where the positive constants and represent the highest rate of predation and the decrease in the urge to hunt, respectively. Numerous investigations of the predator-prey relationship with Ivlev-type functional responses have been conducted. The results suggested that Iviev-type relationships between the species have a number of models in ecological applications, including dynamics in host-parasite models [16], predator-prey models [17–25], animal coat patterns [26], and phytoplankton-zooplankton model [27]. The following predator-prey model with Iviev functional [28] reaction will be taken into consideration.
The parameter refers to the prey’s growth rate, while represents the rate at which prey is turned into a predator after being digested.
Global dynamics, Neimark–Sacker bifurcation, and hybrid control in a Leslie’s prey-predator model were all examined by Khan et al. in [22]. The discrete-time predator-prey model’s stability, bifurcation, and chaos control are examined by the authors in [21], along with the Allee impact on the predator. Santra et al. systematically examined the bifurcation analysis and chaos management of a discrete prey-predator model using a unique prey-refuge idea in [24]. Evaluating the impact of harvesting is realistic when studying the prey-predator paradigm. Fisheries, forestry, and wildlife management all regularly engage in population harvesting [29]. Everyone is aware that one of the most significant systemic changes is immigration. Immigration is a manifestation of an external factor that affects an organism’s ability to create a specific habitat and affects the rate of population increase. Different climatic change responses could disrupt interactions, particularly those between predators and their prey. Occasional roving and nomadic lifestyles are adopted by living species as a result of frequent and seasonal movement and immigration [30, 31]. Additionally, immigration is a significant development that helps to stabilize the environment. In order to better understand this discrete-time predator-prey model with continuous prey population immigration with harvesting on predators, let us look at it.where the time-dependent functions and stand in for the prey and predator population densities, respectively. The highest rate of predation and the declining desire to hunt are represented by the positive constants and , respectively. The predator’s mortality rate is represented by the constant . The prey’s growth rate is indicated by the parameter , while the rate at which the prey becomes a predator after being digested is indicated by the value . Also, and represent the immigration rate on prey population and harvesting rate on predator population, respectively.
Any degree of classical differentiation and integration is generalized in fractional calculus. It is used in a variety of scientific and engineering domains, including biology, fluid mechanics, and medicine, which sparks a great deal of curiosity among academics. Due to its use in numerous domains [32–37] during the last two decades, fractional order calculus has drawn the attention of researchers. Numerous writers have recently studied biological models [38–45] with fractional order. The primary factor is that fractional order models are inherently connected to memory-based systems, which are present in the majority of biological systems [46, 47]. A fractional-order prey-predator model was presented by Javidi and Nyamoradi [48] and its biological behaviors were discussed. The dynamics in parameter spaces of two logistic population maps that are linearly connected and have the same growth rate have recently been studied by Layek and Pati [49].
We apply the Caputo fractional derivatives to the continuous system (5) in the current study and provide a theoretical explanation of the bifurcation occurrences. There are various definitions for fractional derivatives. Caputo’s definition of fractional derivatives, which is frequently used in real-world settings, is one of the most well-known definitions.
Definition 1. Considerwhere denotes the derivative of in the -order, is the value of rounded up to the nearest integer, and is the operator for the Riemann–Liouville integral of -order.where is the gamma function of Euler. The operator is also known as the -order Caputo differential operator.
The following is the model (5)’s fractional order form
There are a variety of techniques for discretizing the model, like a model (8). The piece-wise constant approximation (PCA) [50–52] is one of them. Using the PCA technique, the model is discretized. The steps are as follows:
Assume that model (8) initial conditions are . The discretized version of model (8) is given as follows:
First, let , so . Thus, we obtain
The answer to (10) is simplified to
Second, let , so . Then,which have the following solution:where . After times repetition, we obtainwhere . For , model (14) becomes
The following are a few contributions that this research makes:(1)There are two interdependent species in the planned model, each of which is a source of sustenance for the other. In this study, we looked at how immigration affected the community of prey and harvesting affected the community of predator in the model.(2)Potential fixed points are looked for when evaluating the stability of the model in question.(3)It has been demonstrated that the proposed model can undergo PD and NS bifurcations.(4)The model has become chaotic as a result of the Neimark–Sacker bifurcation, so the OGY (Ott, Grebogi, and Yorke) and state feedback control approaches have been used to manage it.(5)In order to verify the accuracy of our theoretical findings, some numerical examples for our fractional order discrete-time predator-prey model with immigration and harvesting have been provided.
The following sections of this paper are structured as follows: The fixed point topological classes are examined in Sect. 2. In Sect. 3, we investigate analytically the chance of a PD or NS bifurcation of the model (15) under a given parametric condition. In Sect. 4, we numerically show model dynamics with bifurcation diagrams and phase portraits to support our analytic results. In Sect. 5, we apply the OGY and state feedback management techniques to stabilize the chaos of the chaotic model. A short discussion is provided in Sect. 6.
2. Stability of Fixed Point
The three fixed points of model (15) areand , where , which, for any valid parameter value except few parameters must satisfy , always exist. The fixed points and are nonnegative.
Model (15)’s fixed point ’s variational matrix is provided bywhere
At , the characteristic equation becomeswhere
So, and . The eigenvalues of the characteristic equation (19) are , which are real with . The following lemma is one we make in relation to the stability criterion of .
Lemma 2. For the axial fixed point , the topological classification listed below is appropriate:(i)Source if (ii)Sink if (iii)Non-hyperbolic if (iv)Saddle if otherwise
At , the characteristic equation becomeswhere
So, and . The following lemma is one we make in relation to the stability criterion of .
Lemma 3. For the coexistence fixed point , the topological classification listed below is appropriate:(i)Source if(i.i) and (i.ii) and (ii)Sink if(ii.i) and (ii.ii) and (iii)Nonhyperbolic if(iii.i) and (iii.ii) and .(iv)Saddle if otherwise
Letwith
Also, let
The model (15)’s positive fixed point’s topological categorization for with and is depicted in Figure 1. The eigenvalues of the associated Jacobian matrix at the fixed point are complex in the left side of the curve in the green color-shaded section. On the opposite side, the eigenvalues are real. Figure 2 shows the maximum Lyapunov exponents and the bifurcation diagram.
(a)
(b)
(c)
3. Analysis of Bifurcation
This section introduces a study that will look at the NS and PD bifurcations at the fixed point of the model using parameter as a bifurcation parameter.
3.1. Neimark–Sacker Bifurcation
For the analysis of NS bifurcation for the positive coexistence equilibrium point , we consider the following study.
Assume is the small perturbation of where . Therefore, the model becomes using perturbation
If , then the fixed point becomes the origin, and by using Taylor series at expanding and to the third order, model (25) becomeswhere
The characteristic equation of the model (26) is , where and , with .
The characteristic equation’s roots are .
From , and , we have and .
Furthermore, it is imperative that when , , , which is equivalent to .
To study the normal form, let and . We ascertain , and serving the transformation , the model (26) becomeswhere the functions and , respectively, denote the terms in the model (28) for the variables with the order at least two.
To transit through NSB, the following discriminating amount must be nonzero:where
We came to the following conclusion as a result of the analysis indicated above.
Theorem 4. “If , the model proceeds through NS bifurcation at for the parameter to vary in the vicinity of . If , a smooth closed invariant curve with a positive fixed point can bifurcate and the bifurcation is subcritical (resp. supercritical).”
3.2. Period-Doubling Bifurcation
The positive fixed point ’s one eigenvalue is , and the other one is neither 1 nor , if the the model’s parameters vary around the set .
The PD bifurcation is examined using parameter . Furthermore, (,) is this model perturbation caused by the fluctuations of .
If , then equilibrium becomes the origin, and by using Taylor series about expanding to the third order of and , model (31) becomeswhere
We assume , which is invertible. Now, using the transformation , model (32) becomeswhere the functions and , respectively, denote the terms in model (34) for the variables with the order at least two.
System (34) has a center manifold at (0, 0) in a highly closed neighborhood of , which can be deduced using the center manifold theorem and is essentially expressed as follows.where
The model (34)’s restrained center manifold, , has the following form:where
In order for PD bifurcation to occur, the two differentiating quantities and must both be nonzero,
Finally, the outcome of the analysis above is as follows.
Theorem 5. The model experiences PD bifurcation at for varying amounts of in a limited neighborhood of if and . Furthermore, for , the period-two orbits that bifurcate from is stable (unstable).”
4. Numerical Simulations
This section uses the bifurcation diagram, phase portrait, Lyapunov exponent, and fractal dimension to demonstrate the qualitative dynamical properties of the discrete fractional-order model for various parameter values. Numerical simulations will be performed to support the theoretical findings we reached for the model (15). The following parameter values were chosen: and varies between . We observe a positive fixed point and the bifurcation point for the model (15) is calculated at .
Figure 3 depicts the model trajectory as evolving from a fixed point to a PD bifurcation and then to a chaotic attractor. Figures 3(c)–3(d) display the computed MLEs and FDs related to Figures 3(a) and 3(b). The phase portraits are shown in Figure 4 with reference to the bifurcation Figure 3, which essentially depicts the bifurcation of a smooth, invariant closed curve into a chaotic attractor from a stable fixed point.
(a)
(b)
(c)
(d)
The orbit diagram of the prey and predator populations is given in Figure 5, along with other fixed parameter values are and varies between . We observe a positive fixed point and the bifurcation point for the model 7 is calculated at . Figure 5 illustrates the trajectory’s progression from a fixed point to an NSB and then to a chaotic attractor. Figures 5(a), 5(b), and 6 illustrate the phase portrait, MLEs, and FD of Figures 5(a) and 5(b), respectively. There are three separate periodic windows for each of the bifurcation processes for both prey and predator. We have also studied NS bifurcation, and the associated NS bifurcation diagram, MLEs, and FDs are shown in Figure 7. This was done by altering the fractional order in the range and fixing all other parameters indicated above for Figure 5 with .
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
The prey-predator model in the NS bifurcation diagram may behave more dynamically if other parameter values vary (for example, parameter . A fresh Neimark–Sacker bifurcation diagram is produced when the parameter values are set to and varies between , as illustrated in Figures 8(a) and 8(b). The model experiences a Neimark–Sacker bifurcation at . Also, another NS bifurcation diagram is produced when the parameters values are set as and varies between , as illustrated in Figures 9(a) and 9(b). The co-dimension-2 bifurcation diagrams in -space are shown in Figure 2(a). Figure 2(b) displays the plot of the maximal Lyapunov exponents for two control parameters through a 2D projection onto the plane.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
4.1. Fractal Dimension
The fractal dimensions (FD) measurement, which is used to determine a model’s chaotic attractors defined by [53]where the biggest integer number is such that and and ’s are Lyapunov exponents.
Now, the structure of the fractal dimensions of model (15) is as follows:
Given that the chaotic dynamics of the model (15) (ref. Figure 6) are quantified with the sign of FD (ref. Figure 5(d)), it is guaranteed that the dynamics of the fractional order prey-predator model become unstable as the value of the parameter rises.
4.2. 0-1 Chaos Test Algorithm
Assume that is a measured discrete set of data for the chaos test technique [54–56], where , and the total amount of the data is .
We generate new coordinates by choosing a random number and doing the following:where .
The following definition applies to mean square displacement as of this moment:
We also define as the modified mean square displacement as follows:
The following description follows for the median correlation coefficient value:wherein which , and covariance and variance of vectors of length are defined as follows:
The output can now be demonstrated as follows.(i)The dynamics remain stable (i.e., periodic or quasi periodic) when , whereas suggests that the dynamics are chaotic.(ii)As opposed to Brownian-like (unbounded) trajectories, which show chaotic dynamics, bounded trajectories on the plane show regular dynamics (i.e., periodic or quasi-periodic dynamics).
Example 1. Select the parameter values so that with , the Brownian-like (unbounded) trajectories in new coordinates -plane displaying in Figure 10(a) are compatible with a chaotic system dynamics. The curve verses plotted in Figure 10(b) represent the correlation coefficient value.
(a)
(b)
4.3. Co-Dimension-2 Bifurcation
For the parameter values , the parametric space presented in Figure 11(a), we detected the following bifurcations:(i)The fixed point of model (15) tolerates a 1 : 2 resonance bifurcation for ad (ii)The fixed point of model (15) tolerates a 1 : 3 resonance bifurcation for ad (iii)The fixed point of model (15) tolerates a 1 : 4 resonance bifurcation for ad
(a)
(b)
(c)
4.4. Biological Implications
Bifurcations in discrete prey-predator models can have important ecological effects. These models describe the connections between a set of creatures that are preyed upon and another set of species that engage in predation. In these scenarios, the animals that are preyed upon are consumed by the creatures that engage in hunting.
Period-doubling bifurcations and Neimark–Sacker bifurcations are observed in dynamical systems, including ecological models. These divides have important biological effects and can provide insight into the stability and complexity of ecological systems. Period-doubling bifurcations occur when the oscillation period in population dynamics doubles, leading to important ecological effects. Period-doubling bifurcations often signal the transition from regular, repeating patterns to random behaviour in a system. In the context of ecological models, this change could indicate a reduction in the capacity to generate precise forecasts and the emergence of complex, unpredictable fluctuations in population levels. Period-doubling bifurcations are associated with the creation of periodic orbits, which consist of stable cycles of different durations. From an ecological standpoint, this can be viewed as the variation between different population cycles, such as the periodic changes in prey and predator populations that have variable lengths. By studying these models, we can gain knowledge about the basic mechanisms that impact population cycles and other ecological processes and develop more effective methods to improve ecosystem stability and resilience.
Neimark–Sacker bifurcations are associated with the transition from periodic to quasi-periodic behaviour in dynamical systems. In ecological models, this may suggest a shift from simple, regular population cycles to more intricate, nonrepetitive patterns. The occurrence of Neimark–Sacker bifurcations results in quasi-periodic oscillations in the ecological system. The oscillations do not recur exactly, which adds complexity to the temporal dynamics of interacting species. Generally, the Neimark–Sacker bifurcation in discrete prey-predator models highlights the importance of understanding the dynamics of populations and their interactions in ecological systems. By examining these models, we can acquire an understanding of the basic mechanisms that influence population cycles and other ecological processes and create better methods to enhance ecosystem stability and resilience.
Therefore, we may interpret that period-doubling and Neimark–Sacker bifurcations in ecological models as signifying transitions in system behaviour from simple and predictable to complex and potentially chaotic dynamics. These categorizations offer insight into the stability, strength, and flexibility of ecological systems, emphasizing the importance of considering complex dynamics and bifurcation theory when studying population interactions.
Strong resonance bifurcation is a form of bifurcation that can happen in discrete prey-predator models. It is marked by the occurrence of regular or unpredictable changes in population dynamics, which can have significant ecological consequences.
In ecological terms, the strong resonance bifurcation can cause the development of population cycles with high amplitudes, which can lead to more intense boom-and-bust cycles in the ecosystem. This can have significant impacts on the stability of the ecosystem, as the bigger fluctuations in population size can lead to increased competition for resources, predation pressure, and other ecological interactions. Strong resonance Bifurcation can also result in the development of chaotic patterns in population cycles. Unpredictable changes in the ecosystem can arise from chaotic dynamics, making it challenging to establish effective conservation and management methods. In general, the strong resonance bifurcation in discrete prey-predator models emphasizes the need to comprehend the dynamics of populations and their interactions in ecological systems. By examining these models, we can acquire an understanding of the fundamental mechanisms that influence population cycles and other ecological processes, and create more efficient approaches to enhance ecosystem stability and resilience.
5. Chaos Control
According to a performance criterion, dynamical systems are regarded to be the best since they prevent chaos. Chaotic behavior is studied in many disciplines, including physics, biology, ecology, and telecommunications. Additionally, a variety of industries, including communication systems, physics laboratories, biochemistry, turbulence, and cardiology, can benefit from the application of effective chaos management techniques. Recently, the difficulty of controlling chaos dynamics in discrete-time systems has caught the attention of many scholars.
When addressing the issue of managing chaos, the four approaches for investigating chaos control in discrete-time models that are most frequently addressed are the state feedback method, pole-placement methodology, OGY technique, and hybrid control approach. In the model of prey-predator with fractional order, we introduce OGY [57] and state feedback [58] for controlling chaos. The OGY approach does not allow us to employ the control parameter . Using as a control parameter, the OGY technique is put into practice.
We can change the model (15) as illustrated below in order to apply the OGY method.where is the chaotic control parameter. Assume that the chaotic region is defined as , where and symbolize the nominal parameter. The system (48) in the vicinity of the unstable fixed point at can be represented by the following linear map if the model (15) has a chaotic zone created by the expansion of an NS bifurcation at that contains an unstable fixed point.wherewhere and
The controllability matrix of the system (48) is consequently defined as follows:
Consequently, it is easy to deduce that the rank of is 2. We consider that where , then system (48) becomes
Additionally, the suitable controlled system is provided by (15).
If both of the eigenvalues of the matrix’s eigenvalues are situated inside an open unit disk, the fixed point is also locally asymptotically stable.
Also,where .
Furthermore,where
The equations and can then be solved to yield the lines of marginal stability. These limitations also guarantee that the open unit disc has both eigenvalues. Taking into account the case , , and consecutively, we obtain the following equations, respectively, from (56):
The triangular region in the plane circumscribed by the straight lines and then has stable eigenvalues for a given parametric value.
State feedback control, a method, is used to stabilize chaos at the moment where the system’s (15) unstable paths begin. The system (15) can be made to take on a controlled form by introducing a feedback control law as the control force and using the following formula.where represents the nonnegative fixed point of the system (15). The values and indicate the feedback gains.
Example 2. To talk about the system (15)’s OGY feedback control mechanism, we set . In this situation, the unstable system (15) has a single non-negative fixed point . Then, we present the controlled system below based on these parametric parameters.where . We also obtainThen, it is easy to confirm that the matrix’s rank is 2. As a result, the system (60) can be controlled and provides the managed system’s Jacobian matrix.For marginal stability, it offers the lines and .We conducted numerical simulations (see Figures 12–14) to study how the state feedback control influence functions as a chaos controller in an unstable environment. The parameters will be set to the same values as the OGY method that we choose except . The selected feedback gains are and .
(a)
(b)
(c)
(d)
6. Conclusions
In the current study, the dynamics of a fractional-order prey-predator model are investigated, and three fixed points are found under particular parametric conditions. Our findings offer a thorough analysis of the stability of these fixed points, which is given in the article in great detail. Additionally, we both analytically and quantitatively show that the model system can experience period-doubling and Neimark–Sacker bifurcations under specific circumstances. We have numerically computed strong resonance , and bifurcations respectively. Notably, our results show that the system becomes unstable as the parameters , and increase, leading to a bifurcation from a stable state to chaotic behaviour. We see in the simulations the ensuing chaotic behaviour. Furthermore, we notice the effects of harvesting and immigration of the behaviour of the model. For example, low immigration on prey and higher harvesting on predators cause unstable model dynamics while high immigration on prey and lower harvesting on predators stabilize the model dynamics. We also show, numerically and analytically, that the OGY method can be used to regulate chaotic behaviour.
Our major discovery is that the behaviour of the system is strongly influenced by the amount of memory represented by the parameter . Our findings specifically show that weak memory, which corresponds to nearing one, causes chaotic behaviour while strong memory, which corresponds to approaching zero, stabilizes the system. These results underline how crucial memory is to the model system’s behaviour.
In summary, this study provides a thorough analysis of the dynamics of a model system and demonstrates the occurrence of bifurcations and chaos under specific parametric conditions. We also emphasize the influence of memory on the behaviour of the system and demonstrate the efficiency of the OGY method in controlling chaotic behaviour. Our research advances knowledge of the model system’s dynamics and sheds light on the function of memory in the system’s behaviour.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
The main findings were verified by all authors, who also gave their unanimous approval to the final manuscript.