
Citation: Brion Woroch, Alex Konkel, Brian D. Gonsalves. Activation of stimulus-specific processing regions at retrieval tracks the strength of relational memory[J]. AIMS Neuroscience, 2019, 6(4): 250-265. doi: 10.3934/Neuroscience.2019.4.250
[1] | Ceyu Lei, Xiaoling Han, Weiming Wang . Bifurcation analysis and chaos control of a discrete-time prey-predator model with fear factor. Mathematical Biosciences and Engineering, 2022, 19(7): 6659-6679. doi: 10.3934/mbe.2022313 |
[2] | Xiaoling Han, Xiongxiong Du . Dynamics study of nonlinear discrete predator-prey system with Michaelis-Menten type harvesting. Mathematical Biosciences and Engineering, 2023, 20(9): 16939-16961. doi: 10.3934/mbe.2023755 |
[3] | Jinxing Zhao, Yuanfu Shao . Bifurcations of a prey-predator system with fear, refuge and additional food. Mathematical Biosciences and Engineering, 2023, 20(2): 3700-3720. doi: 10.3934/mbe.2023173 |
[4] | Mianjian Ruan, Xianyi Li, Bo Sun . More complex dynamics in a discrete prey-predator model with the Allee effect in prey. Mathematical Biosciences and Engineering, 2023, 20(11): 19584-19616. doi: 10.3934/mbe.2023868 |
[5] | Yuanfu Shao . Bifurcations of a delayed predator-prey system with fear, refuge for prey and additional food for predator. Mathematical Biosciences and Engineering, 2023, 20(4): 7429-7452. doi: 10.3934/mbe.2023322 |
[6] | Xiaoyuan Chang, Junjie Wei . Stability and Hopf bifurcation in a diffusivepredator-prey system incorporating a prey refuge. Mathematical Biosciences and Engineering, 2013, 10(4): 979-996. doi: 10.3934/mbe.2013.10.979 |
[7] | Christian Cortés García, Jasmidt Vera Cuenca . Impact of alternative food on predator diet in a Leslie-Gower model with prey refuge and Holling Ⅱ functional response. Mathematical Biosciences and Engineering, 2023, 20(8): 13681-13703. doi: 10.3934/mbe.2023610 |
[8] | Christian Cortés García . Bifurcations in a discontinuous Leslie-Gower model with harvesting and alternative food for predators and constant prey refuge at low density. Mathematical Biosciences and Engineering, 2022, 19(12): 14029-14055. doi: 10.3934/mbe.2022653 |
[9] | Rajalakshmi Manoharan, Reenu Rani, Ali Moussaoui . Predator-prey dynamics with refuge, alternate food, and harvesting strategies in a patchy habitat. Mathematical Biosciences and Engineering, 2025, 22(4): 810-845. doi: 10.3934/mbe.2025029 |
[10] | Tingting Ma, Xinzhu Meng . Global analysis and Hopf-bifurcation in a cross-diffusion prey-predator system with fear effect and predator cannibalism. Mathematical Biosciences and Engineering, 2022, 19(6): 6040-6071. doi: 10.3934/mbe.2022282 |
The predator-prey system is one of the most important systems for studying the interaction of two species in ecology. Predator-prey systems have important mathematical consequences because they describe ecosystem dynamics and the interactions of diverse species. These systems assist in the estimation of population dynamics by depicting the interaction between predators and prey using mathematical equations. Furthermore, they can inform ecologists by providing insights into the complex interactions between various species within an ecosystem. This knowledge can then be used to develop strategies aimed at the conservation of biodiversity and the effective management of resources [1]. Lotka [2] and Volterra [3] established a fundamental predator-prey system consisting of two species. Over time, several scholars have made modifications to this system to offer a more realistic explanation and improve understanding, as it fails to account for numerous real-world scenarios and complexities. To increase the predator-prey dynamic's authenticity, several ecological principles have been incorporated. These principles include the Allee effect, functional response, refuge-seeking behavior, cannibalism, harvesting impact, and interactions between predators and prey that are mediated by fear [4,5,6,7,8,9,10,11,12].
Numerous researchers have utilized the logistic map to illustrate the prey's growth [13,14,15,16]. Nonetheless, there is a lack of research on the stability analysis of a discrete predator-prey system that takes into account the growth of the prey using a Ricker map [17,18,19]. The logistic map in a one-dimensional population growth model is given by $ x_{n+1} = rx_n(1-\frac{x_n}{k}) $, while the Ricker map is defined as $ x_{n+1} = rx_n e^{1-\frac{x_n}{k}} $. One evident unrealistic feature of the logistic map is that $ 1-\frac{x_n}{k} $ is negative for $ x_n > k $, implying that large populations become negative at the next time step. In contrast, the Ricker map is preferable, as large values of $ x_n $ result in extremely small (but still positive) values of $ x_{n+1} $. Thus, if a population exceeds its carrying capacity, it will fall to extremely low levels, while some of the population survives. Another advantage of the Ricker map is that the exponential component $ e^{1-\frac{x_n}{k}} $ provides a nonlinear response to population density changes, mimicking instances where prey populations might face abrupt declines due to predation pressure.
It is typical to represent dynamical systems in one of two ways when modeling them: i) either as continuous-time systems [20], which are described using differential equations, ii) or as discrete-time systems, which are described by difference equations. Throughout the years, scholars have conducted thorough investigations into the nonlinear dynamic properties exhibited by continuous systems. Recently, numerous researchers have paid significant attention to discrete-time systems [21,22,23,24,25,26,27]. This is because discrete systems are much more effective at facilitating nonoverlapping generations than continuous systems. Discrete-time systems have the advantage of making numerical solutions easy to obtain. The study in [28] describes a precise discrete-time analytical (DTA) signal processing method for estimating frequency and phasor that works well with real-time computing requirements. The study in [29] is primarily concerned with determining ways to compute state and output bounding sets for uncertain discrete-time systems with pointwise-bounded, persistent inputs. The authors in [30] proposed the k-symbol discrete-time fractional Lozi system (FLS). Several critical dynamics of these systems are examined. They also investigate the necessary and sufficient requirements for stable and asymptotically stable k-symbol fractional dynamical systems. Furthermore, substantial research suggests that discrete-time systems may display more complex dynamics than corresponding continuous-time systems [31,32,33,34,35,36,37,38,39,40].
There are two distinct approaches to obtaining a discrete system. One way is to start with a continuous system and then use different techniques, such as the Euler technique [41,42,43,44,45] and the piecewise constant argument method [46,47,48,49,50,51], to turn it into a discrete system. On the other hand, we begin the analysis directly with the discrete system. Hamada et al. [52] studied the following discrete predator-prey system with the Ricker-type growth function:
$ {xn+1=rxne1−xnk−bxnyn,yn+1=dxnyn, $
|
(1.1) |
where $ x_n $ denotes prey density, $ y_n $ is predator density, $ r $ is the intrinsic growth rate of the prey, $ k $ is the environmental carrying capacity of prey, and $ b x_n y_n $ and $ d x_n y_n $ represent the predator-prey confrontation, respectively, which are useful for predators and harmful for prey. The parameters $ r, k, b, $ and $ d $ are positive constants.
To enhance their chances of survival and minimize predation risks, prey species often engage in active refuge-seeking behaviors. This phenomenon has a significant impact on the dynamics of predator-prey interactions because it acts as a crucial mechanism in the preservation and mitigation of the extinction risk that prey species face. Gonzalez-Olivares and Ramos-Jiliberto [53] presented prey refuges in a simple predator-prey system. Ma et al. [54] investigated the dynamic behaviors of a predator-prey system, considering the mutual interference of a predator and a prey refuge. Chen et al. [55] investigated the prey refuge in a Leslie-Gower predator-prey model. Molla et al. [56] investigated the stability and Hopf bifurcation of the predator-prey system with refuge on prey. Numerous researchers have conducted extensive research on the refuge effect and acquired some fascinating results [57,58,59,60,61,62,63]. According to literature studies, the change in prey refuge positively affects prey density and negatively influences predator density. For instance, increasing prey refuge leads to a rise in the prey population, while decreasing prey refuge results in a decrease in the prey population. Negative influence signifies an inverse correlation, with one quantity increasing as the other decreases. Moreover, it has both stabilizing and destabilizing effects. Our work supports previous studies [64,65,66] by demonstrating that both prey and predators benefit from a moderate refuge level.
Thus, motivated by the above discussion, we naturally want to know: When a refuge effect is added to the prey population in system (1.1), what will happen to the dynamical properties? Hence, we extend the system (1.1) by adding the refuge effect to the prey population. Thus, the following modified system is obtained:
$ {xn+1=rxne1−xnk−b(1−m)xnyn,yn+1=d(1−m)xnyn. $
|
(1.2) |
Here $ (1-m)x_n $ represents the quantity of prey available for predation, where $ 0 < m < 1 $ is the protection rate of the prey refuge for prey.
The remainder of the paper is formatted as follows: Section 2 investigates the presence and topological classification of fixed points. Section 3 explores the period-doubling (PD) and Neimark-Sacker (NS) bifurcation analysis at the positive fixed point. Section 4 applies two control methods to regulate bifurcations and chaos. To verify and describe the theoretical results, Section 5 presents some numerical examples. Section 6 discusses the influence of refuge on system (1.2). Lastly, our analysis is summarized in Section 7.
Understanding the stability of fixed points is critical in a predator-prey system. These fixed points depict equilibrium states in which predator and prey populations have reached a balance. Analyzing their stability allows us to forecast the long-term behavior of these ecological systems and provide insight into how different elements influence the overall dynamics of the ecosystem.
Proposition 2.1. For system (1.2), we have three types of fixed points:
1) The trivial fixed point $ E_0 = (0, 0) $ always exists.
2) The predator-free fixed point $ E_1 = (k(ln(r)+1), 0) $ exists if $ r > \frac{1}{e} $.
3) The coexistence fixed point $ E_2 = \bigg(\frac{1}{d(1-m)}, \frac{re^{1-\frac{1}{kd(1-m)}}-1}{b(1-m)}\bigg) $ exists if $ r > e^{-1+\frac{1}{kd(1-m)}} $.
Proof. To determine the fixed points of system (1.2), we need to solve
$ x=rxe1−xk−b(1−m)xy, $
|
(2.1) |
$ y=d(1−m)xy. $
|
(2.2) |
From Eq (2.2), it follows that either $ y = 0 $ or $ x = \frac{1}{d(1-m)} $. Substituting $ y = 0 $ into Eq (2.1), we obtain
$ x=rxe1−xk. $
|
(2.3) |
From Eq (2.3), it follows that either $ x = 0 $ or $ x = k(ln(r)+1) $. Next, substituting $ x = \frac{1}{d(1-m)} $ into Eq (2.1), we obtain
$ y = \frac{re^{1-\frac{1}{kd(1-m)}}-1}{b(1-m)}. $ |
The eigenvalues of the Jacobian matrix help determine the stability of fixed points. If $ \xi_1, \xi_2 $ are eigenvalues of the Jacobian matrix, then $ (x, y) $ is a sink (locally asymptotically stable (LAS)) when $ |\xi_1| < 1 $ along with $ |\xi_2| < 1 $. The fixed point $ (x, y) $ is a source when $ |\xi_1| > 1 $ along with $ |\xi_2| > 1 $. The fixed point $ (x, y) $ is a saddle point (SP) when $ |\xi_1| < 1 \wedge |\xi_2| > 1 $ (or $ |\xi_1| > 1 \wedge |\xi_2| < 1 $). Moreover, the fixed point $ (x, y) $ is a non-hyperbolic point (NHP) when the absolute value of either $ \xi_1 $ and $ \xi_2 $ is one. Classifying the positive fixed point directly using eigenvalues is not easy. Thus, we employ the following result:
Lemma 2.2. [67]
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 \wedge K_0 < 1 $,
2) $ |\xi_1| < 1 \wedge |\xi_2| > 1 $ (or $ |\xi_1| > 1 \wedge |\xi_2| < 1 $) if $ \Lambda(-1) < 0 $,
3) $ |\xi_{1, 2}| > 1 $ if $ \Lambda(-1) > 0 \wedge K_0 > 1 $,
4) $ |\xi_2|\neq 1 \wedge \xi_1 = -1 $ if $ \Lambda(-1) = 0 \wedge 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 \wedge K_0 = 1 $.
Through simple computations, one can obtain that:
$ J(x,y)=[e1−xkr(k−x)k+b(−1+m)yb(−1+m)x−d(−1+m)y−d(−1+m)x]. $
|
Proposition 2.3. The trivial fixed point $ E_0 $ is a
1) LAS if $ 0 < r < \frac{1}{e} $,
2) SP if $ r > \frac{1}{e} $,
3) NHP if $ r = \frac{1}{e} $.
Proof. We obtain
$ J(E0)=[er000]. $
|
(2.4) |
The diagonal entries $ \xi_1 = 0 $ and $ \xi_2 = e r > 0 $ are the eigenvalues $ J(E_0) $. Clearly $ |\xi_1| < 1 $ and
$ er{<1 if 0<r<1e,=1 if r=1e,>1 if r>1e. $
|
Proposition 2.4. The fixed point $ E_1 $ is
1) LAS if $ \frac{1}{e} < r < min\{e, e^{\frac{1}{dk(1-m)}-1}\} $,
2) source if $ r > max\{e, e^{\frac{1}{dk(1-m)}-1}\} $,
3) SP if $ mix\{e, e^{\frac{1}{dk(1-m)}-1}\} < r < max\{e, e^{\frac{1}{dk(1-m)}-1}\} $,
4) NHP if any one of the following satisfies:
(i) $ r = e $,
(ii) $ r = e^{\frac{1}{dk(1-m)}-1} $.
Proof. We obtain
$ J(E1)=[−ln(r)bk(−1+m)(1+ln(r))0dk(1−m)(1+ln(r))]. $
|
(2.5) |
The eigenvalues of $ J(E_0) $ are $ \xi_1 = -ln(r) $ and $ \xi_2 = dk(1-m)(1+ln(r)) > 0 $. One can see that
$ |−ln(r)|{<1 if 1e<r<e,=1 if r=e,>1 if r>e. $
|
Similarly, we obtain
$ dk(1−m)(1+ln(r)){<1 if 1e<r<e1dk(1−m)−1,=1 if r=e1dk(1−m)−1,>1 if r>e1dk(1−m)−1. $
|
Next, we classify the positive fixed point $ E_2 $ of system (1.2) using the Jacobian matrix $ J(x, y) $ and Lemma 2.2.
Theorem 2.5. The positive fixed point
1) $ E_2 $ is LAS if any one of the following satisfies:
(i) $ d < \frac{1}{k(1-m)} $ and
$ e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-1+d k (1-m)}\bigg) < r < -3e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-2+d k (1-m)}\bigg), $
(ii) $ \frac{1}{k(1-m)} < d < \frac{2}{k(1-m)} $ and
$ r < min \bigg\{e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-1+d k (1-m)}\bigg), -3e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-2+d k (1-m)}\bigg)\bigg\} $,
(iii) $ d > \frac{2}{k(1-m)} $ and
$ -3e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-2+d k (1-m)}\bigg) < r < e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-1+d k (1-m)}\bigg), $
2) $ E_2 $ is an SP if one of the following satisfies:
(i) $ d < \frac{2}{k(1-m)} $ and $ r > -3e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-2+d k (1-m)}\bigg) $,
(ii) $ d > \frac{2}{k(1-m)} $ and $ r < -3e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-2+d k (1-m)}\bigg) $,
3) $ E_2 $ is a source if any one of the following satisfies:
(i) $ d > \frac{2}{k(1-m)} $ and
$ r > max \bigg\{e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-1+d k (1-m)}\bigg), -3e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-2+d k (1-m)}\bigg)\bigg\}, $
(ii) $ \frac{1}{k(1-m)} < d < \frac{2}{k(1-m)} $ and
$ e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-1+d k (1-m)}\bigg) < r < -3e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-2+d k (1-m)}\bigg), $
(iii) $ d < \frac{1}{k(1-m)} $ and
$ r < min \bigg\{e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-1+d k (1-m)}\bigg), -3e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-2+d k (1-m)}\bigg)\bigg\} $,
4) $ E_2 $ is NHP and experiences PD bifurcation if
$ r = -3e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-2+d k (1-m)}\bigg) $ and
$ d \neq \frac{2}{k(1-m)}, \ r \neq 2dk(1-m)e^{-1-\frac{1}{d k (-1+m)}}, \ 4dk(1-m)e^{-1-\frac{1}{d k (-1+m)}}. $
5) $ E_2 $ is NHP and experiences NS bifurcation if
$ r = e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-1+d k (1-m)}\bigg), d \neq \frac{1}{k(1-m)} $ and $ 0 < r < 4dk(1-m)e^{-1-\frac{1}{d k (-1+m)}} $.
Proof. We obtain
$ J(E2)=[1+e1+1dk(−1+m)rdk(−1+m)−bdd(−1+e1+1dk(−1+m)r)b1]. $
|
(2.6) |
The corresponding characteristic polynomial is
$ \Lambda(\xi) = \xi^2+K_1 \xi+K_0, $ |
where
$ K1=−2−e1+1dk(−1+m)rdk(−1+m), K0=e1+1dk(−1+m)(1+dk(−1+m))rdk(−1+m). $
|
It can be obtained through calculations that
$ Λ(0)=e1+1dk(−1+m)(1+dk(−1+m))rdk(−1+m),Λ(−1)=3+e1+1dk(−1+m)(1+2dk(−1+m))r,Λ(1)=−1+e1+1dk(−1+m)r. $
|
It is easy to see that the positivity of the $ y $-coordinate of $ E_2 $ implies that $ \Lambda(1) > 0 $. By setting $ \Lambda(-1) = 0 $, one can obtain that:
$ e1+1dk(−1+m)(1+2dk(−1+m))r=−3,(−2+dk(1−m)dk(1−m))r=−3e−1−1dk(−1+m),r=−3e−1−1dk(−1+m)(dk(1−m)−2+dk(1−m)), d≠2k(1−m). $
|
By setting $ \Lambda(0) = 1 $, one can obtain that:
$ e1+1dk(−1+m)(1+dk(−1+m))rdk(−1+m)=1,(−1+dk(1−m))rdk(1−m)=e−1−1dk(−1+m),r=e−1−1dk(−1+m)(dk(1−m)−1+dk(1−m)), d≠1k(1−m). $
|
By setting $ K_1\neq 0, 2 $, we obtain that:
$ r \neq 2dk(1-m)e^{-1-\frac{1}{d k (-1+m)}}, \ 4dk(1-m)e^{-1-\frac{1}{d k (-1+m)}}. $ |
Next, by setting $ K_1^2-4K_0 < 0 $ and $ K_0 = 1 $, we obtain that:
$ 0 < r < 4dk(1-m)e^{-1-\frac{1}{d k (-1+m)}}. $ |
The fixed point categorizations in a discrete-time predator-prey model possess distinct ecological interpretations. A sink represents a state of steady coexistence, a saddle shows a state of intermittent stability, an unstable source implies unexpected shifts in population, and non-hyperbolic points hint at complicated and difficult-to-predict interactions. Understanding these categorizations assists ecologists in comprehending the stability and dynamics of predator-prey interactions, which are vital for efficient ecosystem management and conservation.
This section is dedicated to conducting a thorough investigation of PD and NS bifurcation in system (1.2) at $ E_2 $. To get a comprehensive examination of bifurcation analysis, we suggest the readers to [68,69,70,71,72,73,74,75,76,77,78,79,80,81]. These bifurcations signify important changes in the dynamics of the system, providing insights into situations in which minor changes to parameters result in major changes in the dynamics of predator-prey interactions. Additionally to enhance our understanding of ecosystem dynamics, knowing the roles of PD and NS bifurcations also makes it easier to develop efficient conservation and management methods to maintain the long-term coexistence of predator and prey populations.
In this section, we investigate the PD bifurcation at $ E_2 $ under condition 4) stated in Theorem 2.5. By introducing a minimal perturbation $ \delta $ ($ |\delta|\lll 1 $) to the bifurcation parameter $ r $ in system (1.2), the resulting system is obtained:
$ {xn+1=(r+δ)xne1−xnk−b(1−m)xnyn,yn+1=d(1−m)xnyn. $
|
(3.1) |
Assume that $ u_n = x_n-\frac{1}{d(1-m)}, \ v_n = y_n-\frac{(r+\delta) e^{1-\frac{1}{kd(1-m)}}-1}{b(1-m)} $. After substituting the value of $ r = -3e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-2+d k (1-m)}\bigg) $, the system (3.1) is simplified to
$ [un+1vn+1]=[−1+dk(−1+m)2+dk(−1+m)−bd−2d(1+2dk(−1+m))b(2+dk(−1+m))1][unvn]+[F(un,vn,δ)G(un,vn,δ)], $
|
(3.2) |
where
$ F(un,vn,δ)=a1u2n+a2u3n+a3unvn+a4unδ+a5u2nδ+O((|un|+|vn|+|δ|)4),G(un,vn,δ)=b1unvn+b2unδ, $
|
$ a1=(3+6dk(−1+m))2k(2+dk(−1+m)), a2=−(1+3dk(−1+m))2k2(2+dk(−1+m)), a3=b(−1+m), a4=e1−1dk−dkmdk(−1+m), $
|
$ a5=−e1−1dk−dkm(1+2dk(−1+m))2dk2(−1+m), b1=d(1−m), b2=de1+1dk(−1+m)b. $
|
Next, the system (3.2) is diagonalized through the consideration of the following transformation:
$ [unvn]=[−−2b+bdk−bdkmd(1−2dk+2dkm)−b2d11][enfn], $
|
(3.3) |
Upon applying the mapping (3.3), the system (3.2) undergoes the alteration as follows:
$ [en+1fn+1]=[−1003+3dk(−1+m)2+dk(−1+m)][enfn]+[Γ(en,fn,δ)Υ(en,fn,δ)], $
|
(3.4) |
where
$ Γ(en,fn,δ)=c1e2n+c2enf2n+c3e3n+c4e2nfn+c5f3n+c6enfn+c7f2n+c8e2nδ+c9enδ+c10enfnδ+c11fnδ+c12f2nδ+O((|en|+|fn|+|δ|)4),Υ(en,fn,δ)=d1enf2n+d2e3n+d3e2nfn+d4f3n+d5enfn+d6e2n+d7f2n+d8e2nδ+d9enfnδ+d10f2nδ+d11fnδ+d12enδ+O((|en|+|fn|+|δ|)4), $
|
$ c1=b(2+dk(−1+m))(3+dk(−1+m))dk(5+4dk(−1+m)), c2=−3b2(1+3dk(−1+m))4d2k2(5+4dk(−1+m)), $
|
$ c3=−b2(2+dk(−1+m))2(1+3dk(−1+m))d2k2(1+2dk(−1+m))2(5+4dk(−1+m)), $
|
$ c4=3b2(2+dk(−1+m))(1+3dk(−1+m))2d2k2(1+2dk(−1+m))(5+4dk(−1+m)), $
|
$ c5=b2 (1+2dk(−1+m))(1+3dk(−1+m))8d2k2(2+dk(−1+m))(5+4dk(−1+m)), c6=−3b(2+3dk(−1+m))2dk(5+4dk(−1+m)), $
|
$ c7=−b(1+2dk(−1+m))(−3−2dk(−1+m)+2d2k2(−1+m)2)4dk(2+dk(−1+m))(5+4dk(−1+m)), $
|
$ c8=−b e1−1dk−dkm(2+dk(−1+m))2d2k2(5+4dk(−1+m))(−1+m), c9=e1+1dk(−1+m)(2+dk(−1+m))2dk(5+4dk(−1+m))(−1+m), $
|
$ c10=be1−1dk−dkm(2+dk(−1+m))(1+2dk(−1+m))d2k2(5+4dk(−1+m))(−1+m), $
|
$ c11=−e1+1dk(−1+m)(2+dk(−1+m))(1+2dk(−1+m))2dk(5+4dk(−1+m))(−1+m), $
|
$ c12=−be1−1dk−dkm(1+2dk(−1+m))24d2k2(5+4dk(−1+m))(−1+m), $
|
$ d1=3b2(1+3dk(−1+m))4d2k2(5+4dk(−1+m)), d2=b2(2+dk(−1+m))2(1+3dk(−1+m))d2k2(1+2dk(−1+m))2(5+4dk(−1+m)), $
|
$ d3=−3b2 (2+dk(−1+m))(1+3dk(−1+m))2d2k2(1+2dk(−1+m))(5+4dk(−1+m)), $
|
$ d4=−b2(1+2dk(−1+m))(1+3dk(−1+m))8d2k2(2+dk(−1+m))(5+4dk(−1+m)), $
|
$ d5=3b(1+dk(−1+m)+d2k2(−1+m)2)dk(1+2dk(−1+m))(5+4dk(−1+m)), $
|
$ d6=−3b(2+dk(−1+m))(1+4dk(−1+m)+2d2k2(−1+m)2)dk(1+2dk(−1+m))(5+4dk(−1+m)), $
|
$ d7=3b(−1+4dk(−1+m)+8d2k2(−1+m)2+4d3k3(−1+m)3)4dk(2+dk(−1+m))(5+4dk(−1+m)), $
|
$ d8=b e1−1dk−dkm(2+dk(−1+m))2d2k2(5+4dk(−1+m))(−1+m), $
|
$ d9=−be1−1dk−dkm(2+dk(−1+m))(1+2dk(−1+m))d2k2(5+4dk(−1+m))(−1+m), $
|
$ d10=be1−1dk−dkm(1+2dk(−1+m))24d2k2(5+4dk(−1+m))(−1+m), d11=e1−1dk−dkm(1−d2k2(−1+m)2)dk(5+4dk(−1+m))(−1+m), $
|
$ d12=2e1−1dk−dkm(2+dk(−1+m))(−1+d2k2(−1+m)2)dk(1+2dk(−1+m))(5+4dk(−1+m))(−1+m). $
|
Next, we determine the center manifold denoted by $ Q^C $ for the system (3.4) at the origin, in a close neighborhood to $ \delta = 0 $. Using the center manifold theorem, we can derive the following approximate expression for the center manifold $ Q^C $:
$ QC={(en,fn,δ)∈R3+|fn=p1e2n+p2enδ+p3δ2+O((|en|+|δ|)3)}, $
|
where
$ p1=d61−ξ, p2=−d121+ξ, p3=0, $
|
where $ \xi = \frac{3+3 d k (-1+m)}{2+d k (-1+m)} $. As a result, the system (3.4) is limited to $ Q^C $ in the manner as follows:
$ ˜F:=en+1=−en+c1e2n+c9enδ+(c3−c6d6−1+ξ)e3n−(c11d121+ξ)enδ2+(c8−c11d6−1+ξ−c6d121+ξ)e2nδ+O((|en|+|δ|)4). $
|
(3.5) |
For the function (3.5) to go through PD bifurcation, the following two quantities must possess nonzero values:
$ l1=˜Fδ˜Fenen+2˜Fenδ|(0,0)=2c9, $
|
(3.6) |
$ l2=12(˜Fenen)2+13˜Fenenen|(0,0)=2(c3+c21+c6d61−ξ). $
|
(3.7) |
Based on the aforementioned study, the following result is obtained:
Theorem 3.1. Assume that condition 4) of Theorem 2.5 is satisfied. The system (1.2) experiences PD bifurcation at $ E_2 $ if $ l_1, l_2 $ given in (3.6) and (3.7) are nonzero and $ r $ changes in a close neighborhood of $ r = -3e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-2+d k (1-m)}\bigg) $. Moreover, if $ \ l_2 > 0 $ (respectively $ l_2 < 0 $), then a period-2 orbit of the system (1.2) emerges and it is stable (respectively, unstable).
The above result demonstrates how small changes may produce a significant change in the system's behavior, resulting in a doubling of population oscillation periods. This result discloses an important component of the predator-prey relationship, revealing a transition point in the ecosystem from orderly and predictable cycles to chaotic and unpredictable dynamics.
In this section, we investigate the NS bifurcation at $ E_2 $ under condition (5) stated in Theorem 2.5. By introducing a minimal perturbation $ \delta $ ($ |\delta|\lll 1 $) to the bifurcation parameter $ r $ in system (1.2), the resulting system is obtained:
$ {xn+1=(r+δ)xne1−xnk−b(1−m)xnyn,yn+1=d(1−m)xnyn. $
|
(3.8) |
Assume that $ u_n = x_n-\frac{1}{d(1-m)}, \ v_n = y_n-\frac{(r+\delta) e^{1-\frac{1}{kd(1-m)}}-1}{b(1-m)} $. After substituting the value of $ r = e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-1+d k (1-m)}\bigg) $, the system (3.8) is simplified to
$ [un+1vn+1]=[b11−bdd(−1+e1+1dk(−1+m)(1+dk(−1+m))δ)b+bdk(−1+m)1][unvn]+[F(un,vn)G(un,vn)], $
|
(3.9) |
where
$ b11=e−1dk(1−m)(d2e1dk(1−m)k2(−1+m)2+eδ+dk(−1+m)(2e1dk(1−m)+eδ))dk(1+dk(−1+m))(−1+m), $
|
$ F(un,vn)=b(−1+m)unvn−e−1dk(1−m)(1+2dk(−1+m))(eδ+dk(−1+m)(e1dk(1−m)+eδ))2dk2(1+dk(−1+m))(−1+m)u2n+e−1dk(1−m)(1+3dk(−1+m))(eδ+dk(−1+m)(e1dk(1−m)+eδ))6dk3(1+dk(−1+m))(−1+m)u3n+O((|un|+|vn|)4),G(un,vn)=d(1−m)unvn, $
|
The characteristic equation of the linearized system (3.9) is
$ ξ2−α(δ)ξ+β(δ)=0, $
|
(3.10) |
where
$ α(δ)=−e−1dk−dkm(−2d2e1dk−dkmk2(−1+m)2−eδ−dk(−1+m)(3e1dk−dkm+eδ))dk(1+dk(−1+m))(−1+m),β(δ)=1+e1+1dk(−1+m)(δ+δdk(−1+m)). $
|
The solutions of (3.10) are
$ ξ1,2=α(δ)2±i2√4β(δ)−α2(δ). $
|
(3.11) |
Moreover, we obtain
$ \bigg(\frac{d|\xi_1|}{d \delta}\bigg)_{\delta = 0} = \bigg(\frac{d|\xi_2|}{d \delta}\bigg)_{\delta = 0} = \frac{1}{2} e^{1+\frac{1}{d k (-1+m)}} \bigg(\frac{d k (1-m)-1}{d k (1-m)}\bigg) > 0. $ |
Additionally, it is required that $ \xi_{1, 2}^i \neq 1 $ when $ \delta = 0 $ for $ i = 1, 2, 3, 4, $ which corresponds to $ \alpha(0) \neq -2, 2, 0, 1 $. We obtain
$ α(0)=3+2dk(−1+m)1+dk(−1+m)=2−1−1+dk(1−m)<2. $
|
Moreover, $ \alpha(0) \neq -2, 0, 1 $ is equivalent to
$ d≠54k(1−m),32k(1−m),2k(1−m). $
|
(3.12) |
Next, to change (3.9) into normal form at $ \delta = 0 $, we use the following similarity transformation:
$ [unvn]=[−bd0−12+2dk(−1+m)−√−5+4dk(1−m)2+2dk(−1+m)][enfn]. $
|
(3.13) |
Upon application of the mapping (3.13), the system (3.9) takes the following form:
$ [en+1fn+1]=[3+2dk(−1+m)2+2dk(−1+m)−√−5+4dk(1−m)2+2dk(−1+m)√−5+4dk(1−m)2+2dk(−1+m)3+2dk(−1+m)2+2dk(−1+m)][enfn]+[Γ(en,fn)Υ(en,fn)], $
|
(3.14) |
where
$ Γ(en,fn)=b2dke2n+b2(1+3dk(−1+m))6d2k2(1+dk(−1+m))e3n−b√−5+4dk(1−m)(−1+m)2+2dk(−1+m)enfn+O((|en|+|fn|)4),Υ(en,fn)=b(−1+2dk(−1+m))2dk√−5+4dk(1−m)e2n−b2(1+3dk(−1+m))6d2k2√−5+4dk(1−m)(1+dk(−1+m))e3n+b(3+2dk(−1+m))(−1+m)2+2dk(−1+m)enfn+O((|en|+|fn|)4). $
|
Next, we need the following discriminatory value $ L $ to be not zero to make sure that system (1.2) undergoes NS bifurcation.
$ L=(−Re((1−2ξ1)ξ221−ξ1τ20τ11)−12|τ11|2−|τ02|2+Re(ξ2τ21))δ=0, $
|
(3.15) |
where
$ τ20=18(Γee−Γff+2Υef+i(Υee−Υff−2Γef)), τ11=14(Γee+Γff+i(Υee+Υff)),τ02=18(Γee−Γff−2Υef+i(Υee−Υff+2Γef)),τ21=116(Γeee+Γeff+Υeef+Υfff+i(Υeee+Υeff−Γeef−Γfff)). $
|
Therefore, the result derived from the above analysis is as follows:
Theorem 3.2. Suppose that condition 5) of Theorem 2.5 is satisfied. If the condition (3.12) is satisfied and $ L $ given in (3.15) holds a nonzero value, then system (1.2) experiences NS bifurcation at $ E_2 $ as long as $ r $ varies in a close neighbourhood of $ r = e^{-1-\frac{1}{d k (-1+m)}}\bigg(\frac{d k (1-m)}{-1+d k (1-m)}\bigg) $. Furthermore, in instances where $ L $ is negative (alternatively, positive), the NS bifurcation encountered in system (1.2) 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).
The above result illustrates that, under certain conditions, the predator-prey system experiences an NS bifurcation at point $ E_2 $. This finding indicates a transition in the ecosystem from simple to more complex patterns, resulting in the presence of consistent, non-repeating cycles. Understanding the NS bifurcation enables ecologists to identify the start of enduring, nonlinear fluctuations in the ecosystem, hence facilitating the assessment of long-term population dynamics and ecological stability.
In a predator-prey model, real-world factors serve as control parameters, influencing population dynamics. Environmental changes, such as changes in vegetation or landscape, have an impact on both predators and prey by influencing shelter, food availability, and reproductive success. Introducing a competitor species influences both populations, whether it be new prey for the predator or a competing predator for the prey. Human activities such as hunting rules, conservation initiatives, and harvesting have a direct impact on population size and relationships. Climate elements, such as temperature and precipitation, operate as control variables, influencing birth, mortality, and migration patterns.
Control theory may be employed to control population dynamics in a predator-prey model. It is possible to avoid overpopulation and the extinction of species by maintaining a sustainable and balanced ecosystem by the adjustment of factors such as hunting limits or habitat protection. White-tailed deer populations in the US are managed by hunting limitations to minimize overpopulation and habitat destruction [82]. The Great Barrier Reef Marine Park Authority in Australia prioritizes coral ecosystem maintenance and habitat protection for marine biodiversity [83]. Community-based natural resource management in Namibia promotes sustainable activities like controlled hunting, benefitting wildlife and livelihoods [84]. These examples demonstrate how hunting limitations and habitat conservation affect ecological balance and biodiversity globally.
The objective of control theory is to create management plans that guarantee the populations of prey and predators will coexist in the long run. Bifurcations and unstable oscillations have historically been thought of negatively in mathematical biology since they harm the biological population's ability to reproduce. One can create a controller that may alter the bifurcation characteristics for some non-linear systems to obtain certain desired dynamical properties and manage chaos under the impact of PD and NS bifurcations. There are several strategies for chaos control in a discrete-time system. This section focuses on two different types of control strategies: state feedback control and hybrid control approaches. Both methods are effective in controlling bifurcation and chaos. The hybrid control method is easy to implement. The controlled system in the hybrid control method preserves the fixed points of the original system, while in feedback control, the controlled system may preserve only one fixed point at which we want to control bifurcation and chaos. There is only one control parameter $ \rho \in (0, 1) $ in the hybrid control method, while there are two control parameters ($ \kappa_1, \kappa_2 \in \mathbb{R} $) in the feedback control method.
The feedback control technique [85,86] involves transforming the chaotic system into a piecewise linear system to derive an optimal controller that reduces the upper limit. Subsequently, the optimization issue is performed subject to specified constraints. The aforementioned technique is employed to achieve stabilization of chaotic orbits located at an unstable fixed point inside the system (1.2). The controlled system under consideration for this purpose is as follows:
$ {xn+1=rxne1−xnk−b(1−m)xnyn−Un,yn+1=d(1−m)xnyn, $
|
(4.1) |
where $ U_n = \kappa_1 \bigg(x_n-\frac{1}{d(1-m)}\bigg)+\kappa_2 \bigg(y_n-\frac{re^{1-\frac{1}{kd(1-m)}}-1}{b(1-m)}\bigg) $ is the feedback controlling force, $ \kappa_1 $ and $ \kappa_2 $ are feedback gains. Through simple calculations, it is obtained that for system (4.1), we have
$ J(E2)=[1−κ1+e1−1dk(1−m)rdk(−1+m)−b+dκ2dd(−1+e1+1dk(−1+m)r)b1]. $
|
(4.2) |
The matrix $ J(E_2) $ has the following characteristic equation:
$ ξ2+K1ξ+K0=0, $
|
(4.3) |
where
$ K1=−2+κ1−e1+1dk(−1+m)rdk(−1+m),K0=−κ1+e1+1dk(−1+m)(1+dk(−1+m))rdk(−1+m)+dκ2(−1+e1+1dk(−1+m)r)b. $
|
Let $ \xi_1 $ and $ \xi_2 $ are the roots of (4.3), then we have
$ ξ1+ξ2=2−κ1+e1+1dk(−1+m)rdk(−1+m), $
|
(4.4) |
$ ξ1ξ2=−κ1+e1+1dk(−1+m)(1+dk(−1+m))rdk(−1+m)+dκ2(−1+e1+1dk(−1+m)r)b. $
|
(4.5) |
The marginal stability lines may be found by solving the systems of equations $ \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
$ L1:−κ1+(d(−1+e1+1dk(−1+m)r)b)κ2−1+e1+1dk(−1+m)(1+dk(−1+m))rdk(−1+m)=0. $
|
(4.6) |
Next, we take $ \xi_1 = 1 $ and utilizing Eqs (4.4) and (4.5), we obtain
$ L2:(d−de1+1dk(−1+m)rb)κ2+1−e1+1dk(−1+m)r=0. $
|
(4.7) |
Next, we take $ \xi_1 = -1 $ and utilizing Eqs (4.4) and (4.5), we obtain
$ L3:−2κ1+(d(−1+e1+1dk(−1+m)r)b)κ2+3+e1+1dk(−1+m)(r−2rdk(1−m))=0. $
|
(4.8) |
The stable eigenvalues are enclosed within the triangular region bounded by $ L_1, L_2 $, and $ L_3 $.
The hybrid control technique [87] is a method that combines state feedback and parameter modification to stabilize unstable periodic orbits contained in the system's chaotic attractor. As a result, the regulated system retains its stability over a wide variety of parameters. We take the following controlled system:
$ {xn+1=ρ(rxne1−xnk−b(1−m)xnyn)+(1−ρ)xn,yn+1=ρd(1−m)xnyn+(1−ρ)yn, $
|
(4.9) |
where $ \rho \in (0, 1) $. The parameter $ \rho $, acting like a control parameter, balances the impact of the original system (1.2) with the modified system (4.9). If the value of $ \rho $ becomes negative, it might indicate the reverse impact of the original system (1.2). Conversely, if $ \rho $ exceeds $ 1 $, it could indicate an amplified effect of the original system (1.2) beyond its natural influence, perhaps leading to unrealistic or unworkable consequences in the modified system (4.9). The same fixed points are shared by systems (4.9) and (1.2). We obtain
$ J(E2)=[−1+m+e1+1dk(−1+m)rρdk−1+m−bρdd(−1+e1+1dk(−1+m)r)ρb1], $
|
(4.10) |
with corresponding characteristic polynomial
$ Λ(ξ)=ξ2+K1ξ+K0, $
|
(4.11) |
where
$ K1=2−2m−e1+1dk(−1+m)rρdk−1+m,K0=1+e1+1dk(−1+m)rρdk(−1+m)+(−1+e1+1dk(−1+m)r)ρ2. $
|
Theorem 4.1. The fixed point $ E_2 $ of the system (4.9) is LAS if
$ |K_1| < 1+K_0 < 2. $ |
Remark 4.2. These control strategies aim to mitigate bifurcation and chaos in the system (1.2). The mathematical equations in systems (4.1) and (4.9) define parameters $ \kappa_1, \kappa_2 $ and $ \rho $ in the context of control techniques. It is important to note that these specific control methods may not have direct, established parallels in current ecological models or practices. Our approach introduces theoretical modifications, and we acknowledge the need for further research and practical applications within the field of mathematical ecology to fully validate these methods.
In this section, we will corroborate our theoretical findings for system (1.2) by numerical simulations. These numerical simulations will include bifurcation diagrams, phase portraits, time series plots, and maximum Lyapunov exponent (MLE) graphs. We have used MATHEMATICA for computations and MATLAB for graphs.
We assume that $ k = 2.5, b = 1.3, m = 0.5, d = 0.9, x_0 = 2.25, y_0 = 4.45, r \in [3.38, 3.68] $, then, system (1.2) goes through PD bifurcation when $ r \approx 3.451523 $. The positive fixed point is obtained as $ E_2 = (2.222222, 4.395604) $. The eigenvalues of $ J(E_2) $ are $ \xi_1 = -1 $ and $ \xi_2 = -0.428571 $ with $ |\xi_2| \neq 1 $. For these parametric values, we obtain
$ Γ(en,fn,δ)=23.1742e2n−1.39513e3n−55.4374enfn+5.08261e2nfn+33.1465f2n−6.17217enf2n+2.49844f3n+6.23407×10−15δ+31.3355enδ+2.72516e2nδ−38.053fnδ−6.6187enfnδ+4.01878f2nδ+O((|en|+|fn|+|δ|)4),Υ(en,fn,δ)=19.873e2n−1.21015e3n−47.5458enfn+4.4087e2nfn+28.4316f2n−5.35379enf2n+2.16716f3n+5.13358×10−15δ+26.6219enδ+2.36382e2nδ−32.3289fnδ−5.74111enfnδ+3.48592f2nδ+O((|en|+|fn|+|δ|)4). $
|
Thus, we obtain
$ l1=2c9=2×31.3355=62.671>0,l2=2(c3+c21+c6d61−ξ)=2(−1.39513+(23.1742)2+−55.4374×19.8731+0.428571)=−471.0941<0. $
|
The bifurcation diagrams of system (1.2) are given in Figure 1(a), (b), while the MLE is plotted in Figure 1(c).
Next, consider $ k = 2.5, b = 1.3, m = 0.5, d = 1.5, x_0 = 1.3, y_0 = 1.7 $ and varying $ r \in[1.1, 2.1] $. The system (1.2) goes through NS bifurcation at $ r \approx 1.343762 $ and has the positive fixed point $ E_2 = (1.333333, 1.758242) $. The eigenvalues of $ J(E_2) $ are $ \xi_{1, 2} = 0.428571\pm 0.903508 i $ with $ |\xi_{1, 2}| = 1 $. Moreover, some careful calculations give
$ τ20=−0.026310−0.016640i,τ11=0.086667+0.260361i,τ02=0.112976+0.277i,τ21=0.039702+0.02511. $
|
Thus, it is obtained that $ L = -0.078752 < 0 $, which proves the correctness of Theorem 3.2. Bifurcation diagrams are depicted in Figure 2(a), (b), while the MLE is plotted in Figure 2(c). The presence of negative MLEs indicates the presence of chaotic areas.
Next, Figure 3(a)–(h) shows phase portraits of system (1.2) for various various values of $ r $. One can observe that $ E_2 $ is LAS for $ r < 1.343762 $ but loses stability at $ r \approx 1.343762 $ when the system (1.2) goes through NS bifurcation. For $ r \geq 1.343762 $, an invariant curve emerges from $ E_2 $, the radius of which grows as $ r $ grows. Some $ 5-, 10- $ periodic orbits are also plotted in Figure 3(d), (f), Finally, we obtain a strange chaotic attractor given in Figure 3(h).
We assume that $ r = 3.5, k = 2.5, b = 1.3, d = 0.9, x_0 = 1.95, y_0 = 4.50, m \in [0.42, 0.52] $, then, system (1.2) experiences both NS bifurcation and PD bifurcation as $ m $ varies in small neighborhoods of $ m_1 \approx 0.424620 $ and $ m_2 \approx 0.495050 $, respectively. The bifurcation diagrams of system (1.2) are given in Figure 4(a), (b), while the MLE is plotted in Figure 4(c). The presence of negative MLEs indicates the existence of stable fixed points or stable periodic windows, whereas positive MLEs indicate the presence of chaotic areas. Furthermore, the phase portraits of system (1.2) are given in Figure 5(a)–(h) for various values of the parameter $ m $. One can observe that system (1.2) experiences NS bifurcation for small values of refuge. At the NS bifurcation point, an invariant closed curve emerges, representing a repeating pattern in predator-prey populations. Moreover, the system experiences PD bifurcation for large values of refuge. At the PD bifurcation point, the system (1.2) transitions from stable behavior to periodic oscillations, and subsequently, the period of these oscillations doubles.
The existence of two critical values, $ m_1 $ and $ m_2 $, suggests a threshold behavior in the system (1.2). When $ m $ is less than $ m_1 $, the positive fixed point $ E_2 $ is unstable, implying that predator-prey interaction is too skewed in favor of the predators, and the prey population cannot sustain itself. Similarly, when $ m $ is greater than $ m_2 $, then $ E_2 $ is also unstable, indicating that too much refuge availability disrupts the predator-prey balance. This suggests that a moderate level of refuge is beneficial for both predator and prey populations.
The efficacy of the hybrid control approach will next be evaluated. We assume $ \rho = 0.96, r = 3.5, k = 2.5, b = 1.3, d = 0.9, x_0 = 1.95, y_0 = 4.50 $ and vary $ m $ for the controlled system (4.9). If $ 0.405045 < m < 0.500983 $, the positive fixed point $ E_2 $ is LAS. The controlled system's bifurcation diagrams, and Figure 6(a), (b) show that the bifurcation has been postponed in the controlled system (4.9).
Next, we aim to evaluate the efficacy of the feedback control technique. Considering $ r = 3.5, k = 2.5, b = 1.3, d = 0.9 $, and $ m = 0.515 $, as well as the initial conditions $ x_0 = 1.95 $ and $ y_0 = 4.50 $ for the controlled system (4.1), the marginal stability lines are as follows:
$ L_1: \kappa_2 = 0.351067+0.514907 \kappa_1, $ |
$ L_2: \kappa_2 = -1.444444, $ |
and
$ L_1: \kappa_2 = 0.086950+1.02981 \kappa_1. $ |
Figure 7(a) depicts the stability region bounded by lines $ L_1, L_2 $, and $ L_3 $ for system (4.1). The fixed point $ E_2 $ of system (1.2) is shown to be unstable for the given parametric values. The controlled system (4.1) is examined with feedback gains $ \kappa_1 = -2.95 $ and $ \kappa_2 = -1.20 $. Figure 7 illustrates the graph of $ x_n $ as shown in Figure 7(c), $ y_n $ as shown in Figure 7(d), and the phase portrait as presented in Figure 7(b) for the system (4.1). Therefore, it may be deduced that the use of the feedback control methodology seems to be effective in controlling bifurcation and chaos.
Trivial fixed point $ E_0 = (0, 0) $ and boundary fixed point $ E_1 = (k(ln(r)+1), 0) $ are independent of refuge effect $ m $. The positive fixed point $ E_2 = \bigg(\frac{1}{d(1-m)}, \frac{re^{1-\frac{1}{kd(1-m)}}-1}{b(1-m)}\bigg) $ is dependent on refuge effect $ m $. We assume that $ r = 3.5, k = 2.5, b = 1.3, d = 0.9 $. Then, in Table 1, we observe that the density of prey population increases as refuges used by prey increases, while that of predators first increases and then decreases with prey refuges. The same is observed in Figure 8. When the number of prey using refuges is sufficiently large, our findings indicate that the prey population will exceed its maximum environmental carrying capacity, leading to the extinction of predators. Furthermore, the examination of stability and bifurcation demonstrates that an appropriate amount of refuge is advantageous for the populations of both predators and prey.
value of $ m $ | prey population in $ E_2 $ | predator population in $ E_2 $ |
$ 0.1 $ | $ 1.23457 $ | $ 4.10791 $ |
$ 0.3 $ | $ 1.5873 $ | $ 4.44198 $ |
$ 0.5 $ | $ 2.22222 $ | $ 4.47895 $ |
$ 0.7 $ | $ 3.7037 $ | $ 2.98086 $ |
$ 0.9 $ | $ 11.1111 $ | $ -6.83286 $ |
The refuge effect plays a crucial role in determining the stability of predator-prey interactions in an ecosystem. In a predator-prey system, such a refuge could be a physical space, a habitat, or any resource that protects the prey. In this paper, we present and study the complex dynamics of a discrete-time predator-prey system with the refuge effect. The presence and stability of fixed points are investigated. Moreover, a thorough analysis of local bifurcations at the positive fixed point is conducted. The study illustrates that the system (1.2) goes through both PD and NS bifurcation. Moreover, the presence of a positive MLE guarantees the existence of chaos in the system (1.2). Feedback control and hybrid control approaches are used to control bifurcation and chaos. Consequently, effective control is achieved for both types of bifurcation across an extensive range of parameters. Furthermore, numerical simulations are executed to demonstrate the theoretical results that were previously presented. These simulations use several visual representations, including bifurcation diagrams, MLE graphs, phase portraits, and time series plots.
It is observed that a positive fixed point is stable if the refuge parameter $ m $ lies in an optimal range $ m_1 < m < m_2 $. When $ m $ is less than $ m_1 $, the positive fixed point $ E_2 $ is unstable, implying that predator-prey interaction is too skewed in favor of the predators, and the prey population cannot sustain itself. Similarly, when $ m $ is greater than $ m_2 $, then $ E_2 $ is also unstable, indicating that too much refuge availability disrupts the predator-prey balance. This suggests that a moderate level of refuge is beneficial for both predator and prey populations.
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
The authors declare there is no conflict of interest.
This work was supported by Research Funding from Youjiang Medical University for Nationalities, Baise, China under Grant numbers yy2020bsky050 & yy2023rcky002, the National Natural Science Foundation of China under Grant number 62162063, and the Scientific Research and Technology Development Program of Guangxi, China under Grant number 2021AC19308. The funding bodies did not play any role in the design of the study and in writing the manuscript.
[1] | Cohen NJ, Eichenbaum H (1993) Memory, amnesia, and the hippocampal system. Cambridge, MA, US: The MIT Press. |
[2] | Eichenbaum H, Cohen NJ (2001) From conditioning to conscious recollection: Memory systems of the brain. New York: Oxford University Press. |
[3] |
Rissman J, Wagner AD (2012) Distributed representations in memory: insights from functional brain imaging. Annu Rev Psychol 63: 101–128. doi: 10.1146/annurev-psych-120710-100344
![]() |
[4] |
Danker JF, Anderson JR (2010) The ghosts of brain states past: remembering reactivates the brain regions engaged during encoding. Psychol Bull 136: 87–102. doi: 10.1037/a0017937
![]() |
[5] | Rugg MD, Vilberg KL (2014) Brain Networks Underlying Episodic Memory Retrieval. Curr Opin Neurobilogy 23: 255–260. |
[6] | Morcom AM (2014) Re-engaging with the past: recapitulation of encoding operations during episodic retrieval. Front Hum Neurosci 8: 351. |
[7] |
Wheeler ME, Petersen SE, Buckner RL (2000) Memory's echo: Vivid remembering reactivates sensory-specific cortex. Proc Natl Acad Sci 97: 11125–11129. doi: 10.1073/pnas.97.20.11125
![]() |
[8] |
Johnson JD, McDuff SGR, Rugg MD, et al. (2009) Recollection, familiarity, and cortical reinstatement: a multivoxel pattern analysis. Neuron 63: 697–708. doi: 10.1016/j.neuron.2009.08.011
![]() |
[9] | Thakral PP, Wang TH, Rugg MD (2017)Decoding the content of recollection within the core recollection network and beyond. Cortex 91: 101–113. |
[10] |
Johnson JD, Rugg MD (2007) Recollection and the reinstatement of encoding-related cortical activity. Cereb Cortex 17: 2507–2515. doi: 10.1093/cercor/bhl156
![]() |
[11] |
Alvarez P, Squire LR (1994) Memory consolidation and the medial temporal lobe: A simple network model. Proc Natl Acad Sci U S A 91: 7041–7045. doi: 10.1073/pnas.91.15.7041
![]() |
[12] |
McClelland JL, McNaughton BL, O'Reilly RC (1995) Why there are complementary learning systems in the hippocampus and neocortex: insights from the successes and failures of connectionist models of learning and memory. Psychol Rev 102: 419–457. doi: 10.1037/0033-295X.102.3.419
![]() |
[13] |
Moscovitch M, Rosenbaum RS, Gilboa A, et al. (2005) Functional neuroanatomy of remote episodic, semantic and spatial memory: A unified account based on multiple trace theory. J Anat 207: 35–66. doi: 10.1111/j.1469-7580.2005.00421.x
![]() |
[14] |
Norman KA, O'Reilly RC (2003) Modeling hippocampal and neocortical contributions to recognition memory: A complementary-learning-systems approach. Psychol Rev 110: 611–646. doi: 10.1037/0033-295X.110.4.611
![]() |
[15] |
Rosler F, Heil M, Hennighausen E (1995) Distinct Cortical Activation Patterns during Long-Term Memory Retrieval of Verbal , Spatial , and Color Information. J Cogn Neurosci 7: 51–65. doi: 10.1162/jocn.1995.7.1.51
![]() |
[16] |
Khader P, Burke M, Bien S, et al. (2005) Content-specific activation during associative long-term memory retrieval. Neuroimage 27: 805–816. doi: 10.1016/j.neuroimage.2005.05.006
![]() |
[17] |
Nyberg L, Habib R, Herlitz A (2000) Brain activation during episodic memory retrieval: sex differences. Acta Psychol (Amst) 105: 181–194. doi: 10.1016/S0001-6918(00)00060-3
![]() |
[18] |
Vaidya CJ, Zhao M, Desmond JE, et al. (2002) Evidence for cortical encoding specificity in episodic memory: memory-induced re-activation of picture processing areas. Neuropsychologia 40: 2136–2143. doi: 10.1016/S0028-3932(02)00053-2
![]() |
[19] |
Wheeler ME, Buckner RL (2003) Functional dissociation among components of remembering: control, perceived oldness, and content. J Neurosci 23: 3869–3880. doi: 10.1523/JNEUROSCI.23-09-03869.2003
![]() |
[20] |
Wheeler ME, Buckner RL (2004) Functional-anatomic correlates of remembering and knowing. Neuroimage 21: 1337–1349. doi: 10.1016/j.neuroimage.2003.11.001
![]() |
[21] |
Woodruff CC, Johnson JD, Uncapher MR, et al. (2005) Content-specificity of the neural correlates of recollection. Neuropsychologia 43: 1022–1032. doi: 10.1016/j.neuropsychologia.2004.10.013
![]() |
[22] |
Liang JC, Preston AR (2017) Medial temporal lobe reinstatement of content-specific details predicts source memory. Cortex 91: 67–78. doi: 10.1016/j.cortex.2016.09.011
![]() |
[23] | Leiker EK, Johnson JD (2015) Pattern reactivation co-varies with activity in the core recollection network during source memory. Neuropsychologia 75. |
[24] |
Maratos EJ, Dolan RJ, Morris JS, et al. (2001) Neural activity associated with episodic memory for emotional context. Neuropsychologia 39: 910–920. doi: 10.1016/S0028-3932(01)00025-2
![]() |
[25] |
Smith AP, Henson RN, Dolan RJ, et al. (2004) fMRI correlates of the episodic retrieval of emotional contexts. Neuroimage 22: 868–878. doi: 10.1016/j.neuroimage.2004.01.049
![]() |
[26] |
Kahn I, Davachi L, Wagner AD (2004) Functional-neuroanatomic correlates of recollection: implications for models of recognition memory. J Neurosci 24: 4172–4180. doi: 10.1523/JNEUROSCI.0624-04.2004
![]() |
[27] |
Voss JL, Galvan A, Gonsalves BD (2011) Neuropsychologia Cortical regions recruited for complex active-learning strategies and action planning exhibit rapid reactivation during memory retrieval. Neuropsychologia 49: 3956–3966. doi: 10.1016/j.neuropsychologia.2011.10.012
![]() |
[28] |
Dodson CS, Holland PW, Shimamura AP (1998) On the recollection of specific-and partial-source information. J Exp Psychol Learn Mem Cogn 24: 1121–1136. doi: 10.1037/0278-7393.24.5.1121
![]() |
[29] |
Hicks JL, Marsh RL, Ritschel L (2002) The role of recollection and partial information in source monitoring. J Exp Psychol Learn Mem Cogn 28: 503–508. doi: 10.1037/0278-7393.28.3.503
![]() |
[30] |
Simons JS, Dodson CS, Bell D, et al. (2004) Specific- and partial-source memory: effects of aging. Psychol Aging 19: 689–694. doi: 10.1037/0882-7974.19.4.689
![]() |
[31] |
Slotnick SD, Dodson CS (2005) Support for a continuous (single-process) model of recognition memory and source memory. Mem Cognit 33: 151–170. doi: 10.3758/BF03195305
![]() |
[32] |
Mickes L, Wais P, Wixted J (2009) Recollection is a continuous process implications for dual-process theories of recognition memory. Psychol Sci 20: 509–515. doi: 10.1111/j.1467-9280.2009.02324.x
![]() |
[33] |
Slotnick SD (2010) Remember' source memory ROCs indicate recollection is a continuous process. Memory 18: 27–39. doi: 10.1080/09658210903390061
![]() |
[34] |
Vilberg KL, Rugg MD (2008) Memory retrieval and the parietal cortex: A review of evidence from a dual-process perspective. Neuropsychologia 46: 1787–1799. doi: 10.1016/j.neuropsychologia.2008.01.004
![]() |
[35] |
Vilberg KL, Rugg MD (2009) Left parietal cortex is modulated by amount of recollected verbal information. Neuroreport 20: 1295–1299. doi: 10.1097/WNR.0b013e3283306798
![]() |
[36] |
Vilberg KL, Rugg MD (2007) Dissociation of the neural correlates of recognition memory according to familiarity, recollection, and amount of recollected information. Neuropsychologia 45: 2216–2225. doi: 10.1016/j.neuropsychologia.2007.02.027
![]() |
[37] |
Shimamura AP (2011) Episodic retrieval and the cortical binding of relational activity. Cogn Affect Behav Neurosci 11: 277–291. doi: 10.3758/s13415-011-0031-4
![]() |
[38] | Thakral PP, Wang TH, Rugg MD (2015) Cortical reinstatement and the confidence and accuracy of source memory. Neuroimage 109. |
[39] |
Mack ML, Preston AR (2016) Decisions about the past are guided by reinstatement of specific memories in the hippocampus and perirhinal cortex. Neuroimage 127: 144–157. doi: 10.1016/j.neuroimage.2015.12.015
![]() |
[40] |
Ishai A (2008) Let's face it: it's a cortical network. Neuroimage 40: 415–419. doi: 10.1016/j.neuroimage.2007.10.040
![]() |
[41] |
Epstein R, Kanwisher N (1998) A cortical representation of the local visual environment. Nature 392: 598–601. doi: 10.1038/33402
![]() |
[42] |
O'Craven KM, Kanwisher N (2000) Mental imagery of faces and places activates corresponding stiimulus-specific brain regions. J Cogn Neurosci 12: 1013–1023. doi: 10.1162/08989290051137549
![]() |
[43] |
Althoff RR, Cohen NJ (1999) Eye-movement-based memory effect: A reprocessing effect in face perception. J Exp Psychol Learn Mem Cogn 25: 997–1010. doi: 10.1037/0278-7393.25.4.997
![]() |
[44] |
Walther DB, Caddigan E, Fei-Fei L, et al. (2009) Natural scene categories revealed in distributed patterns of activity in the human brain. J Neurosci 29: 10573–10581. doi: 10.1523/JNEUROSCI.0559-09.2009
![]() |
[45] | Brett M, Anton J, Valabregue R, et al. (2002) Region of interest analysis using an SPM toolbox [abstract] Presented at the 8th International Conference on Functional Mapping of the Human Brain, June 2-6, 2002, Sendai, Japan. Available on CD-ROM in NeuroImage, Vol 16. |
[46] |
Warriner AB, Kuperman V, Brysbaert M (2013) Norms of valence, arousal, and dominance for 13,915 English lemmas. Behav Res Methods 45: 1191–1207. doi: 10.3758/s13428-012-0314-x
![]() |
[47] |
Diana RA, Yonelinas AP, Ranganath C (2007) Imaging recollection and familiarity in the medial temporal lobe: A three-component model. Trends Cogn Sci 11: 379–386. doi: 10.1016/j.tics.2007.08.001
![]() |
[48] |
Gonsalves BD, Kahn I, Curran T, et al. (2005) Memory strength and repetition suppression: multimodal imaging of medial temporal cortical contributions to recognition. Neuron 47: 751–761. doi: 10.1016/j.neuron.2005.07.013
![]() |
[49] |
Wais PE (2008) FMRI signals associated with memory strength in the medial temporal lobes: A meta-analysis. Neuropsychologia 46: 3185–3196. doi: 10.1016/j.neuropsychologia.2008.08.025
![]() |
1. | Mohammed Alsubhi, Rizwan Ahmed, Ibrahim Alraddadi, Faisal Alsharif, Muhammad Imran, Stability and bifurcation analysis of a discrete-time plant-herbivore model with harvesting effect, 2024, 9, 2473-6988, 20014, 10.3934/math.2024976 | |
2. | Anum Zehra, Parvaiz Ahmad Naik, Ali Hasan, Muhammad Farman, Kottakkaran Sooppy Nisar, Faryal Chaudhry, Zhengxin Huang, Physiological and chaos effect on dynamics of neurological disorder with memory effect of fractional operator: A mathematical study, 2024, 250, 01692607, 108190, 10.1016/j.cmpb.2024.108190 | |
3. | Aqeel Ahmad, Muhammad Farman, Parvaiz Ahmad Naik, Khurram Faiz, Abdul Ghaffar, Evren Hincal, Muhammad Umer Saleem, Analytical analysis and bifurcation of pine wilt dynamical transmission with host vector and nonlinear incidence using sustainable fractional approach, 2024, 11, 26668181, 100830, 10.1016/j.padiff.2024.100830 | |
4. | Parvaiz Ahmad Naik, Rizwan Ahmed, Aniqa Faizan, Theoretical and Numerical Bifurcation Analysis of a Discrete Predator–Prey System of Ricker Type with Weak Allee Effect, 2024, 23, 1575-5460, 10.1007/s12346-024-01124-7 | |
5. | Kottakkaran Sooppy Nisar, A constructive numerical approach to solve the Fractional Modified Camassa–Holm equation, 2024, 106, 11100168, 19, 10.1016/j.aej.2024.06.076 | |
6. | Chih-Wen Chang, Zohaib Ali Qureshi, Sania Qureshi, Asif Ali Shaikh, Muhammad Yaqoob Shahani, Real-Data-Based Study on Divorce Dynamics and Elimination Strategies Using Nonlinear Differential Equations, 2024, 12, 2227-7390, 2552, 10.3390/math12162552 | |
7. | Lana Abdelhaq, Sondos M. Syam, Muhammad I. Syam, An efficient numerical method for two-dimensional fractional integro-differential equations with modified Atangana–Baleanu fractional derivative using operational matrix approach, 2024, 11, 26668181, 100824, 10.1016/j.padiff.2024.100824 | |
8. | Parvaiz Ahmad Naik, Yashra Javaid, Rizwan Ahmed, Zohreh Eskandari, Abdul Hamid Ganie, Stability and bifurcation analysis of a population dynamic model with Allee effect via piecewise constant argument method, 2024, 70, 1598-5865, 4189, 10.1007/s12190-024-02119-y | |
9. | Mominul Islam, M. Ali Akbar, A study on the fractional-order COVID-19 SEIQR model and parameter analysis using homotopy perturbation method, 2024, 12, 26668181, 100960, 10.1016/j.padiff.2024.100960 | |
10. | Cahit Köme, Yasin Yazlik, Stability, bifurcation analysis and chaos control in a discrete predator–prey system incorporating prey immigration, 2024, 70, 1598-5865, 5213, 10.1007/s12190-024-02230-0 | |
11. | R. Prem Kumar, G.S. Mahapatra, P.K. Santra, Dynamical analysis of SARS-CoV-2-Dengue co-infection mathematical model with optimum control and sensitivity analyses, 2024, 80, 14681218, 104175, 10.1016/j.nonrwa.2024.104175 | |
12. | Parvaiz Ahmad Naik, Zohreh Eskandari, Mehmet Yavuz, Zhengxin Huang, Bifurcation results and chaos in a two-dimensional predator-prey model incorporating Holling-type response function on the predator, 2024, 0, 1937-1632, 0, 10.3934/dcdss.2024045 | |
13. | Parvaiz Ahmad Naik, Muhammad Farman, Anum Zehra, Kottakkaran Sooppy Nisar, Evren Hincal, Analysis and modeling with fractal-fractional operator for an epidemic model with reference to COVID-19 modeling, 2024, 10, 26668181, 100663, 10.1016/j.padiff.2024.100663 | |
14. | Hafiz Muhammad Shahbaz, Iftikhar Ahmad, Muhammad Asif Zahoor Raja, Hira Ilyas, Kottakkaran Sooppy Nisar, Muhammad Shoaib, A novel design of recurrent neural network to investigate the heat transmission of radiative Casson nanofluid flow consisting of carbon nanotubes (CNTs) across a curved stretchable surface, 2024, 104, 0044-2267, 10.1002/zamm.202400104 | |
15. | Ning Tian, Xiaoqi Liu, Rui Kang, Cheng Peng, Jiaxi Li, Shang Gao, Noise-to-State Stability of Random Coupled Kuramoto Oscillators via Feedback Control, 2024, 12, 2227-7390, 3715, 10.3390/math12233715 | |
16. | Yu Mu, Wing-Cheong Lo, Yuanshun Tan, Zijian Liu, Hybrid control for the prey in a spatial prey-predator model with cooperative hunting and fear effect time lag, 2025, 491, 00963003, 129217, 10.1016/j.amc.2024.129217 | |
17. | Allah Ditta, Parvaiz Ahmad Naik, Rizwan Ahmed, Zhengxin Huang, Exploring periodic behavior and dynamical analysis in a harvested discrete-time commensalism system, 2025, 13, 2195-268X, 10.1007/s40435-024-01551-z | |
18. | Subarna Roy, Subhas Khajanchi, Pankaj Kumar Tiwari, Fear and its carry-over effects in a generalist predator–prey system featuring cooperative hunting, 2025, 1598-5865, 10.1007/s12190-024-02346-3 | |
19. | Wei Li, Chunrui Zhang, Codimension-1 and Codimension-2 Bifurcations Analysis of Discrete Predator–Prey Model with Herd Behavior, 2025, 35, 0218-1274, 10.1142/S0218127425500105 | |
20. | Ibraheem M. Alsulami, Rizwan Ahmed, Faraha Ashraf, Exploring complex dynamics in a Ricker type predator–prey model with prey refuge, 2025, 35, 1054-1500, 10.1063/5.0232030 | |
21. | Parvaiz Ahmad Naik, Muhammad Farman, Muhammad Umer Saleem, Zhengxin Huang, Hijaz Ahmad, Muhammad Sultan, 2025, Chapter 11, 978-981-97-6793-9, 163, 10.1007/978-981-97-6794-6_11 | |
22. | Pinar Baydemir, Huseyin Merdan, Bifurcation analysis, chaos control, and FAST approach for the complex dynamics of a discrete-time predator–prey system with a weak Allee effect, 2025, 196, 09600779, 116317, 10.1016/j.chaos.2025.116317 |
value of $ m $ | prey population in $ E_2 $ | predator population in $ E_2 $ |
$ 0.1 $ | $ 1.23457 $ | $ 4.10791 $ |
$ 0.3 $ | $ 1.5873 $ | $ 4.44198 $ |
$ 0.5 $ | $ 2.22222 $ | $ 4.47895 $ |
$ 0.7 $ | $ 3.7037 $ | $ 2.98086 $ |
$ 0.9 $ | $ 11.1111 $ | $ -6.83286 $ |