
A predator–prey dynamic reaction model is investigated in a two-layered water body where only the prey is subjected to harvesting. The surface layer (Layer-1) provides food for both species, while the prey migrates to deeper layer (Layer-2) as a refuge from predation. Although the prey is the preferred food for the predator, the predator can also consume alternative food resources that are abundantly available. The availability of alternative food resources plays a crucial role in species' coexistence by mitigating the risk of extinction. The main objective of the work was to explore the effect of different harvesting strategies (nonlinear and linear harvesting) on a predator–prey model with effort dynamics in a heterogeneous habitat. The analysis incorporates a dual timescale approach: the prey species migrate between the layers on a fast timescale, whereas the growth of resource biomass, prey–predator interactions, and harvesting dynamics evolve on a slow timescale. The complete model involving both slow and fast timescales has been investigated by using aggregated model. The reduced aggregated model is analyzed analytically as well as numerically. Moreover, it is demonstrated that the reduced system exhibits the bifurcations (transcritical and Hopf point) by setting the additional food parameter as a bifurcation parameter. A comparative study using different harvesting strategies found that there is chaos in the system when using linear harvesting in the predator–prey model. However, nonlinear harvesting gives only stable or periodic solutions. This concludes that nonlinear harvesting can control the chaos in the system. Additionally, a one-dimensional parametric bifurcation, phase portraits, and time series plots are also explored in the numerical simulation.
Citation: Rajalakshmi Manoharan, Reenu Rani, Ali Moussaoui. Predator-prey dynamics with refuge, alternate food, and harvesting strategies in a patchy habitat[J]. Mathematical Biosciences and Engineering, 2025, 22(4): 810-845. doi: 10.3934/mbe.2025029
[1] | 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 |
[2] | 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 |
[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] | Peter A. Braza . A dominant predator, a predator, and a prey. Mathematical Biosciences and Engineering, 2008, 5(1): 61-73. doi: 10.3934/mbe.2008.5.61 |
[5] | 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 |
[6] | Shuqi Zhai, Qinglong Wang, Ting Yu . Fuzzy optimal harvesting of a prey-predator model in the presence of toxicity with prey refuge under imprecise parameters. Mathematical Biosciences and Engineering, 2022, 19(12): 11983-12012. doi: 10.3934/mbe.2022558 |
[7] | Claudio Arancibia–Ibarra, José Flores . Modelling and analysis of a modified May-Holling-Tanner predator-prey model with Allee effect in the prey and an alternative food source for the predator. Mathematical Biosciences and Engineering, 2020, 17(6): 8052-8073. doi: 10.3934/mbe.2020408 |
[8] | 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 |
[9] | Ilse Domínguez-Alemán, Itzel Domínguez-Alemán, Juan Carlos Hernández-Gómez, Francisco J. Ariza-Hernández . A predator-prey fractional model with disease in the prey species. Mathematical Biosciences and Engineering, 2024, 21(3): 3713-3741. doi: 10.3934/mbe.2024164 |
[10] | Muhammad Shoaib Arif, Kamaleldin Abodayeh, Asad Ejaz . On the stability of the diffusive and non-diffusive predator-prey system with consuming resources and disease in prey species. Mathematical Biosciences and Engineering, 2023, 20(3): 5066-5093. doi: 10.3934/mbe.2023235 |
A predator–prey dynamic reaction model is investigated in a two-layered water body where only the prey is subjected to harvesting. The surface layer (Layer-1) provides food for both species, while the prey migrates to deeper layer (Layer-2) as a refuge from predation. Although the prey is the preferred food for the predator, the predator can also consume alternative food resources that are abundantly available. The availability of alternative food resources plays a crucial role in species' coexistence by mitigating the risk of extinction. The main objective of the work was to explore the effect of different harvesting strategies (nonlinear and linear harvesting) on a predator–prey model with effort dynamics in a heterogeneous habitat. The analysis incorporates a dual timescale approach: the prey species migrate between the layers on a fast timescale, whereas the growth of resource biomass, prey–predator interactions, and harvesting dynamics evolve on a slow timescale. The complete model involving both slow and fast timescales has been investigated by using aggregated model. The reduced aggregated model is analyzed analytically as well as numerically. Moreover, it is demonstrated that the reduced system exhibits the bifurcations (transcritical and Hopf point) by setting the additional food parameter as a bifurcation parameter. A comparative study using different harvesting strategies found that there is chaos in the system when using linear harvesting in the predator–prey model. However, nonlinear harvesting gives only stable or periodic solutions. This concludes that nonlinear harvesting can control the chaos in the system. Additionally, a one-dimensional parametric bifurcation, phase portraits, and time series plots are also explored in the numerical simulation.
The study of ecology is a captivating field of research that focuses on the biological processes occurring among various living organisms. It explores the interactions between different species and their physical environments. In an ecological system, species rely on one another for essential functions, as they do not always exist in isolation. A particularly intriguing aspect of ecology is the predator–prey model, which highlights the dynamic interactions between two populations.
Interactions between species can affect population size in various ways, including competition for limited resources, predation, food scarcity, unregulated harvesting, and environmental factors such as climate change and pollution. Consequently, changes in the population of one species can impact others within the same ecosystem, potentially leading to species extinction. To prevent species extinction, numerous strategies have been proposed and implemented, such as improving the conditions of natural habitats, limiting harvesting, establishing natural reserves, and creating protected areas. These strategies have been incorporated into various mathematical formulations of ecological systems.
Researchers [1,2,3] have shown that refugia can stabilize predator–prey models and may reduce the risk of extinction for prey. The incorporation of prey refuges into mathematical models enhances the accuracy of analyzing predator–prey dynamics. A "prey refuge" refers to habitats that are inaccessible or difficult for predators to reach. The prey population uses these refuges to reduce the predation risk and increase their chances of survival. The concept of prey refuges is well established in ecological studies, as many aquatic prey species, such as fish, seek shelter in deeper water to escape from predators like seabirds. This behavior has significant implications for fishery management, as it influences predator–prey dynamics and the sustainability of fish populations. One study [4] investigated how prey refuges can modulate the chaotic behavior caused by time delays and contribute to a population's persistence within the system. In [5], the researcher demonstrated the effectiveness of refuges in reducing predation as the prey population evolved.
In many ecological systems, predators typically do not rely on a single prey species for sustenance. Instead, they depend on alternative food sources, including other prey species. The availability of additional food sources can prevent predator extinction. The role of supplementary food in sustaining predator populations has been widely explored in the literatures [6,7,8,9,10]. The author of [6] examined the effects of providing additional food to predators on the dynamics of a predator–prey model with prey refuges, using the framework proposed in [7]. Their analysis revealed that while high levels of prey refuges could lead to predator extinction, survival could be ensured with the support of additional food. Building on the model from [6], Samaddar et al. [8] introduced the concept of fear effects and investigated their influence on system's dynamics in the presence of prey refuges and additional food, considering both its quality and quantity. The researcher of [11] examined how the additional food affects the stability of the three-species food chain model combining prey refuges and harvesting.
On the contrary, it is essential to acknowledge that harvesting is an unavoidable aspect of human activity aimed at obtaining economic resources. Many researchers have investigated systems with harvesting functions to enhance the realism of their models. Within the framework of dynamical systems, various types of harvesting functions have been employed to study the effects of harvesting efforts [12]. For instance, in [13,14], the authors analyzed a system with linear harvesting. Researchers [15,16] investigated a system with linear or constant harvesting and demonstrated that it exhibits significantly more intricate and complex dynamics compared with the system without harvesting, while Chen et al. [17] explored a model proposed in [18] by incorporating nonlinear harvesting and delay, demonstrating that nonlinear harvesting induces more complex dynamical behavior. Research by [19,20] further emphasized the significant influence of nonlinear harvesting terms on dynamical systems. Additionally, in [21], the author analyzed a system with nonlinear harvesting due to its greater realism. This present paper aimed to integrate both linear and nonlinear harvesting functions into the proposed predator–prey system, highlighting their realistic features.
Several investigations [22,23] have considered the effects of refuges within a heterogeneous environment consisting of interconnected patches. This spatial heterogeneity necessitates an examination of two distinct dynamics: Local interactions among species on the one hand, and their migration between different patches on the other. Researchers have extensively studied the density-dependent nature of dispersal in the literature [24,25,26]. One author [24] examined how prey's dispersal behavior, influenced by predator density, affects the stability of the predator–prey system in the presence of refuges. Furthermore, the study in [25] explored the impact of dispersal behavior by both predators and prey on the dynamical system, demonstrating that the dispersal of both species can significantly influence the stability of predator–prey interactions.
In an open-access fishery, fishing effort adjusts in response to the perceived rent—whether positive or negative—reflecting the net economic revenue for fishermen. A model that captures this dynamic interaction between perceived rent and fishing effort is known as a dynamic reaction model.
The current study investigated a dynamic reaction model within a predator–prey fishery system situated in a heterogeneous patchy habitat, where only the prey species is subject to harvesting. The model considers a water body with two layers: The surface layer (Layer-1) and the deeper layer (Layer-2). In this system, food is available only in Layer-1, supporting both predator and prey species. The prey species migrate between these two layers. Due to limited food resources in Layer-2, the prey move to Layer-1, where food is abundant, and it is assumed that their population grows exponentially in this region. The predator species depends on prey for their survival and is assumed to always remain in Layer-1 only. However, to avoid predation, the prey species hide themselves in the safest place in the deeper layer, which is known as a prey refuge where predators cannot reach them. For instance, consider seabirds as predators that feed on fish (prey) for survival, while the fish migrate from Layer-1 to Layer-2 to hide from the predators. This migration of prey from Layer-1 to Layer-2 is influenced by predator density: Higher predator numbers in Layer-1 lead to increased prey movement into the deeper layer. Although the prey is the predator's preferred food, the predator can also consume an alternative food source that is abundantly available. To counteract this, an additional food source is introduced for predators, which is an ecologically relevant aspect observed in marine ecosystems where predators can switch to other available food even when their primary prey is abundant. Harvesting plays a significant role in population dynamics, balancing both ecological and economic stability. This study incorporates both linear and nonlinear harvesting functions, acknowledging that prey are harvested at different rates in each layer due to practical fishing constraints. Some prey that seeks refuge in Layer-2 can still be harvested, though at a lower rate than in Layer-1, since fishing nets cover the entire water body but do not reach deep into the ocean. This study primarily examined the impact of linear and nonlinear harvesting on the prey population across both layers. The model further includes economic considerations such as harvesting effort, which evolves dynamically in response to costs and revenues, making it relevant for fishery management. The system exhibits two distinct timescales: a fast timescale and a slow timescale. The migration of the prey population between layers occurs on the fast timescale, while species growth, species interactions, and changes in fleet size occur on the slow timescale. Incorporating these timescales facilitates the simplification of the complex system using aggregation methods [22,27,28,29], which are based on perturbation techniques and the center manifold theorem [30,31,32,33].
The subsequent sections of this paper are structured as follows: In Section 2, we formulate the complete model and provide a contextual overview of the model. In Section 3, we examine the dynamics of a predator–prey model incorporating alternate food and linear harvesting, while Section 4 extends this analysis in the case of nonlinear harvesting. Lastly, Section 5 summarizes our findings and discusses their implications.
In a heterogeneous patchy environment, the density of the prey in Layer-i (i=1,2) is denoted by ni(t), where Layer-1 serves as a resource layer where the prey gather for feeding despite facing predation, and Layer-2 acts as a refuge for the prey. The predator density in Layer-1 at time t is denoted by p(t), while E represents the harvesting effort, such as number of boats or nets used for harvesting. A schematic representation of the model (2.1) is shown in Figure 1.
It is assumed that the migration of species between two layers occurs much faster than the biological interactions of the system such as the growth of the species population, interactions between the prey and predators, and harvesting effort. To make the system more realistic, the two timescales are introduced: a slow timescale t, which governs the growth, predation, and harvesting processes, and a fast timescale τ, which governs the movement of species between layers. To address the two timescales, a small dimensionless parameter ε is introduced, and the slow timescale is defined as t=ετ, where (0<ε<<1). This scaling allows us to separate the fast dynamics from the slow dynamics, enabling a simplified analysis that focuses on the system's slow evolution while accounting for the rapid migration of species [19,25,26,27].
According to these assumptions, the complete system (slow–fast system) of the predator–prey system, considering an alternate food source for the predators and harvesting at both layers is presented below:
dn1dτ=(kn2−ˆk(p)n1)+ε(r1n1−Aan1p−H1(n1,n2,E))dn2dτ=(ˆk(p)n1−kn2)+ε(−r2n2−H2(n1,n2,E))dpdτ=ε(−μp+Abn1p+β(1−A)p)dEdτ=ε(−cE+p0H1(n1,n2,E)+p0H2(n1,n2,E))} | (2.1) |
n1(0)>0,n2(0)>0,p(0)>0,E(0)>0. |
The constant r1>0 indicates the prey's intrinsic growth rate in the surface layer. The term r1n1 represents the growth rate of the prey in Layer-1, which is assumed to be follow exponential growth, while the term −r2n2 represents the mortality rate of the prey due to the lack of food in Layer 2. Additionally, there is a negative term Aan1p accounting for prey loss due to predation in Layer 1, where the predation rate is given by a for the prey. The functions H1(n1,n2,E) and H2(n1,n2,E) represent the harvesting terms in Layer-1 and Layer-2, respectively. The constant μ is the natural mortality rate for the predators, while b signifies the food conversion rate by the predators with respect to the prey n1(t). Thus, the growth of the predator population is directly proportional to the density of captured prey. The dependence of predators on alternate food resources [10] is considered in predator dynamics through the term β(1−A)p, (0<A<1). When A=1, predators rely solely on prey species in Layer 1. When A=0, the predators depend exclusively on alternate food resources. The parameters c and p0 are assumed to be the cost and price per unit of harvesting, respectively. The term −cE represents the natural reduction in harvesting effort over time due to costs. The harvesting functions p0H1(n1,n2,E) and p0H2(n1,n2,E) represent the revenue or benefit from harvesting.
The parameter k indicates the migration rate of the prey from Layer 2 to Layer 1, which is assumed to be constant. The main motivating factors for this movement of the prey are food and light. The movement of prey from the surface layer to the deeper layer is assumed to be dependent on the density of predators [34] and is denoted by ˆk :
ˆk(p)={αpforp>00forp=0. |
To study the model (2.1), the system is reduced to one timescale (the slow timescale) by substituting fast equilibrium, which is investigated as follows.
To study the fast dispersal model, we neglect the slow part of the system (2.1) by taking ε=0:
dn1dτ=kn2−ˆk(p)n1dn2dτ=ˆk(p)n1−kn2dpdτ=0dEdτ=0 | (2.2) |
The fast equilibria for the fast system (2.2) are obtained as follows:
n∗1=kk+ˆk(p)n=f(p)n=v∗1n;f(p)=kk+ˆk(p)n∗2=ˆk(p)k+ˆk(p)n=(1−f(p))n=v∗2n |
The equilibrium frequencies v∗1 and v∗2 take the following form:
v∗1=f(p)=n∗1nandv∗2=(1−f(p))=n∗2n | (2.3) |
From the expression above (2.3), it is observed that the equilibrium densities (n∗1,n∗2) for the fast system are proportional to the total population. The equilibrium frequencies, denoted as v∗1 and v∗2, represent the proportion of prey in each layer under fast equilibrium. The total frequencies across layers always sum to one, with the sum of their rates of change equating to zero. These equilibrium frequencies vary as functions of the slow variable p, which is assumed to be constant at the fast timescale. For every combination of the slow variable values n, p, and E, the fast system tends towards an equilibrium state. The addition of the first two equations of the system (2.2) gives n(t)=n1(t)+n2(t) as a constant for the fast part of (2.1).
The dynamic behavior of the system (2.1) is explored with linear and nonlinear harvesting functions in Sections 3 and 4.
The dynamics of the corresponding complete predator–prey system, incorporating linear harvesting and an alternate food source for the predators is derived in this section. Here, the harvesting rate remains constant, regardless of the prey's population size. This means that the same number of individuals, or a constant proportion of the population, is harvested at each time step. The harvesting functions in the system (2.1) can be taken as H1(n1,n2,E)=q1En1 and H2(n1.n2,E)=q2En2. It follows that
dn1dτ=(kn2−ˆk(p)n1)+ε(r1n1−Aan1p−q1En1)dn2dτ=(ˆk(p)n1−kn2)+ε(−r2n2−q2En2)dpdτ=ε(−μp+Abn1p+β(1−A)p)dEdτ=ε(−cE+p0q1En1+p0q2En2)} | (3.1) |
n1(0)>0,n2(0)>0,p(0)>0,E(0)>0. |
The parameters q1 and q2 represent the catchability coefficient in Layer-1 and Layer-2, respectively, in which q1En1 and q2En2 represent the linear harvesting terms in Layer-1 and Layer-2, respectively. The aggregated model is further elaborated for the complete system (3.1) in the subsequent section.
Applying the aggregation method, the complete system (3.1) is reduced to a system of three ordinary differential equations. Let n(t)=n1(t)+n2(t) be the aggregated variable. Further, we introduce r(p), q(p), and m(p) as given below:
r(p)=r1k−r2αpk+αp,q(p)=q1k+q2αpk+αp |
The aggregated system presented below is obtained by summing the two prey equations and substituting the fast equilibrium into the complete system (3.1)
dndt=n(r(p)−Aaf(p)p−q(p)E)=n⋅F(p,E)dpdt=p(−μ+Abf(p)n+β(1−A))=p⋅G(n,p)dEdt=E(−c+p0q(p)n)=E⋅H(n,p)} | (3.2) |
n(0)>0,p(0)>0,E(0)>0. |
The aggregated model (3.2) provides an approximation of the complete system (3.1) derived using perturbation theory and the center manifold theorem. The aggregated model (3.2) includes new and different terms with respect to the slow part of the slow–fast system. This is due to the density dependence of the equilibrium frequencies. This process is called functional emergence within the approximated aggregated system.
The system (3.2) possesses four positive equilibrium points:
1) The trivial point P0(0,0,0) is always exists.
2) P1(ˆn,0,ˆE)=P1(cp0q1,0,r1q1), the predator-free boundary's fixed point in the positive nE-plane.
3) The harvesting-free boundary equilibrium point
P2(¯n,¯p,0)=(kβ(A−A0)(r1α+r2α+Aak)(Abk)(Aak+r2α),r1kAak+r2α,0);A0=1−μβ |
exists in the positive np-plane and is positive for the critical value of the alternate food resource as follows:
A>A0. | (3.3) |
4) The unique interior fixed point P3(n∗,p∗,E∗) of (3.2) exists, where n∗,p∗, and E∗ are obtained as shown below:
n∗=c(k+αp∗)p0q1k+p0q2αp∗p∗=k(−bc+p0q1β)(A2−A)αβp0q2(A−A0);A2=β−μβ−bcp0q1E∗=r(p∗)−Aaf(p∗)p∗q(p∗)=r1k−(r2α+Aak)p∗q1k+q2αp∗ |
The interior point P3(n∗,p∗,E∗) is feasible for the condition
A1<A<A2 |
where
A1=β−μβ−r2bcp0(r2q1+r1q2)andA2=β−μβ−bcp0q1 |
It is noted that the dependence of the predators on alternate food (A) is critical for the existence of various equilibrium points. Therefore, the following cases can be investigated:
0<A<A0 | (3.4) |
A0<A<A1 | (3.5) |
A1<A<A2 | (3.6) |
A2<A<1 | (3.7) |
It can be observed that the equilibrium points P0 and P1 may exist irrespective of A. However, the boundary equilibrium point P2 exists for the condition (3.5). There is coexistence of all the species for the conditions in (3.6). The predators may or may not survive under the conditions in (3.4) and (3.7). However, the predators can survive for a suitable range of alternate food, as given in (3.5) and (3.6). Accordingly, the prey and predator populations will not go extinct for the cases of (3.5) and (3.6). It can be observed that for the case of (3.5), no harvesting of prey species is possible because of the nonavailability of a sufficient amount of prey.
This subsection explores the local stability criteria for the feasible equilibrium points of the dynamic system (3.2).
The Jacobian matrix at (n,p,E) for the system described in (3.2) is provided as follows:
J(n,p,E)=[Fn(r′(p)−Aa(f′(p)p+f(p))−q′(p)E)−nq(p)Abf(p)pAbf′(p)pn+G0Ep0q(p)Ep0q′(p)nH] |
where
f′(p)=−αk(k+αp)2<0,r′(p)=−αk(r2+r1)(k+αp)2<0,q′(p)=αk(q2−q1)(k+αp)2>0;q2>q1 | (3.8) |
The Jacobian matrix, when calculated at the equilibrium point P0(0,0,0), takes the following form:
J0(0,0,0)=[r1000−μ+β(1−A)000−c] |
The eigenvalues of J0 are
λ01=r1>0,λ02=−μ+β(1−A)andλ03=−c<0 |
Hence, the origin point P0(0,0,0) consistently manifests as a saddle point, featuring an unstable manifold in the n -direction and a stable manifold in the E -direction. Furthermore, it has a stable manifold in the p -direction for the following condition:
A>1−μβ=A0 | (3.9) |
If the condition in (3.9) is violated, then P0 has an unstable manifold in the p direction. Accordingly, the trajectories along p=0, starting in neighborhood of P0, may be attracted to P1 when the condition in (3.9) is violated.
At the point P1(ˆn,0,ˆE), the Jacobian matrix is given by
J1(ˆn,0,ˆE)=[0c12c130c220c31c320]=[0ˆn(r1−Aa−ˆEkα(q2−q1)k2)−ˆnq10−μ+Abˆn+β(1−A)0ˆEp0q1ˆEp0kˆnα(q2−q1)k20] |
The characteristic equation derived from the Jacobian matrix about the point P1(ˆn,0,ˆE) is provided as follows:
λ3+C11λ2+C12λ+C13=0 | (3.10) |
where C11=−c22,C12=−c13c31,C13=c13c31c22, and C11C12−C13=c22c13c31−c13c31c22=0
Since C11C12=C13, Eq (3.10) may be rewritten as
(λ2+C12)(λ+C11)=0 |
This gives
λ11=c22andλ12,13=±i√C12 |
This confirms the periodic solutions around the equilibrium point P1(ˆn,0,ˆE) in the nE-plane. This indicates that in the absence of predators, the prey population and harvesting effort exhibit cyclic fluctuations instead of reaching a stable equilibrium. These oscillations result from the continuous interaction between the prey's growth and harvesting intensity.
The following Jacobian matrix is obtained for the point P2(¯n,¯p,0)
J2(¯n,¯p,0)=[0d12d13d21d22000d33]=[0¯n(r′(¯p)−Aa(f′(¯p)¯p+f(¯p)))−¯nq(¯p)Abf(¯p)¯pAbf′(¯p)¯p¯n000−c+p0¯nq(¯p)] |
The characteristic equation associated with the matrix J2(¯n,¯p,0) is given by
(λ−(−c+p0¯nq(¯p)))(λ2−(Ab¯pf′(¯p)¯n)λ+Ab¯pf(¯p)¯n(r′(¯p)−Aa(f′(¯p)¯p+f(¯p))))=0 |
One eigenvalue is λ21=−c+p0¯nq(¯p), while the other two eigenvalues can be determined from the following characteristic equation,
λ2−(Ab¯pf′(¯p)¯n)λ+Ab¯pf(¯p)¯n(r′(¯p)−Aa(f′(¯p)¯p+f(¯p)))=0 | (3.11) |
The trace(J2) and det(J2) of Eq (3.11) are computed as follows:
trace(J2)=−αr1β(μ−β(1−A))Aak+r1α+r2α<0det(J2)=r1β(μ−β(1−A))(Aak+r2α)Aak+r1α+r2α>0 |
Accordingly, the point P2(¯n,¯p,0) is locally asymptotically stable for λ21<0 which gives the condition
¯n<cp0q(¯p) | (3.12) |
If the condition of (3.12) is violated, P2(¯n,¯p,0) transitions into a saddle point. This indicates the instability along the z-axis.
The Jacobian matrix computed at the interior point P3(n∗,p∗,E∗) is given by
J3(n∗,p∗,E∗)=[e11e12e13e21e22e23e31e32e33]J3(n∗,p∗,E∗)=[0n∗(r′(p∗)−Aa(f′(p∗)p∗+f(p∗))−q′(p∗)E∗)−n∗q(p∗)Abf(p∗)p∗Abf′(p∗)p∗n∗0E∗p0q(p∗)E∗p0q′(p∗)n∗0] |
The following characteristic equation is determined from the Jacobian matrix above about the point (n∗,p∗,E∗):
λ3+D1λ2+D2λ+D3=0 | (3.13) |
with
D1=−e22=−Abf′(p∗)p∗n∗=Abkp∗n∗(k+αp∗)2>0D2=−e13e31−e12e21=−(−ve)(+ve)−(−ve)(+ve)>0D3=e13e31e22−e13e31e21=Abp0n∗2p∗E∗q(p∗)(αk(q2−q1)(k+αp∗)3+αk(q1k+q2αp∗)(k+αp∗)2)>0 |
Moreover,
D1D2−D3=e12e21e22+e13e31e21=Abkn∗2p∗(k+αp∗)4(kα2Ab(r1+r2)p∗+A2abk2αp∗+(q2−q1)α(k2Abp∗−p0E∗(q1k+q2αp∗)))>0 | (3.14) |
By applying the Routh–Hurwitz criterion, the interior point (n∗,p∗,E∗) is locally asymptotically stable (LAS) if and only if the condition in (3.14) is satisfied. If (3.14) is violated, the point (n∗,p∗,E∗) may become unstable. It should be noted that the condition of (3.14) is always satisfied when q1=q2. In this case, the point P3 is always locally asymptotically stable. In general, it is difficult to analyze the condition in (3.14) to determine the stability. Therefore, it is analyzed for a particular set of data (see Subsection 3.5). The numerical results show that for the given data-set, the system exhibits complex and chaotic dynamics.
Theorem 3.1. The system (3.2) undergoes a transcritical bifurcation around the planar equilibrium point P2(¯n,¯p,0) as the bifurcation parameter A crosses a critical value such that
A=Atc1 | (3.15) |
Example 3.1. For a particular dataset k=0.5,r1=3,r2=1,a=2,α=2.5,q1=1,q2=1.5,μ=2,b=4,β=3,c=4,and p0=10, the Jacobian matrix J2 of P2(0.293797,0.521937,0) is confirmed to have a zero eigenvalue at the threshold value of Atc1=0.373911. The Jacobian matrix J∗2 is given below.
J∗2(P2,Atc1)=[0−0.4678−0.40.2163−0.08800000] | (3.16) |
The eigenvectors corresponding to the zero eigenvalue (i.e., d33=0) are Ve1=(−0.2557,−0.5283,0.7348)T and We1=(0,0,1)T. This analysis confirms that all the conditions of Sotomayar's theorem for a transcritical bifurcation at the point P2 are verified.
Theorem 3.2. The system (3.2) undergoes a transcritical bifurcation around the planar equilibrium point P1(ˆn,0,ˆE) as the bifurcation parameter A crosses a critical value such that
A=Atc2 | (3.17) |
Example 3.2. For a particular dataset k=0.5,r1=3,r2=1,a=2,α=2.5,q1=1,q2=1.5,μ=2,b=4,β=3,c=4, and p0=10, the Jacobian matrix J1 of P1(0.4,0,0.3) is confirmed to have a zero eigenvalue at the threshold value of Atc2=0.714286. The Jacobian matrix J∗1 is given below.
J∗1(P1,Atc2)=[0−2.3714−0.400030300] | (3.18) |
The eigenvectors corresponding to the zero eigenvalue (i.e., c22=0) are Ve2=(−0.1641,0.1641,−0.9727)T and We2=(0,1,0)T. This analysis confirms that all the conditions of Sotomayar's theorem for a transcritical bifurcation at the point P1 are verified.
Let us validate the Hopf bifurcation around an interior equilibrium point with respect to the bifurcation parameter A. Due to the complexity of determining the Hopf point analytically, we will solve this with the help of a numerical example by evaluating the characteristic equation of the Jacobian matrix J3 according to Liu's criterion [35,36,37].
Liu's criterion: The characteristic equation of the Jacobian matrix at an interior equilibrium point is
λ3+ξ1(A)λ2+ξ2(A)λ+ξ3=0. |
It should also satisfy the following conditions for some critical value AH.
1) ξ1(AH)>0, ξ3(AH)>0, Φ(AH)=ξ1(AH)∗ξ2(AH)−ξ3(AH)=0
2) dΦ(AH)dAH≠0
Example 3.3. The system (3.2) undergoes a Hopf bifurcation around an interior point P3(0.301964,0.370331,0.227896) for the following particular data: k=0.5,r1=3,r2=1,a=2,α=2.5,q1=1,q2=1.5,μ=2,b=4,β=3,c=4, and p0=10. This Hopf point occurs when the parameter A crosses a threshold value AH=0.388133. It can be obtained by determining the conditions of Liu's criterion as follows:
ξ1(AH)=0.1067>0andξ3(AH)=0.1459>0Φ(AH)=ξ1(AH)∗ξ2(AH)−ξ3(AH)=0.1067∗1.3673−0.1459=0anddΦ(AH)dAH=−1.4190≠0 |
Since the conditions required by Liu's criterion are satisfied, the Hopf bifurcation occurs around an interior point P3.
Consider the following data for the system (3.2) with linear harvesting:
k=0.5,r1=3,r2=1,a=2,α=2.5,q1=1,q2=1.5,μ=2,b=4,β=3,c=4,p0=10 | (3.19) |
The critical values of A are computed as follows:
A0=0.333,A1=0.37,A2=0.714 |
For the aggregated system (3.2), the boundary equilibrium P1=(4,0,0.3) is obtained. If we choose A=0.35 (A>A0), another boundary equilibrium point P2=(0.1297,0.5263,0) exists. Next, for A=0.38 (A>A1), the interior equilibrium point P3=(0.29967804,0.43207403,0.029233241) exists.
The dynamics of the system were additionally examined through the utilization of the MATLAB software package [38,39,40,41] alongside MATCONT. MATCONT is a toolbox implemented in MATLAB, consisting of numerical algorithms aimed at detecting, continuing, and identifying limit cycles, also referred to as periodic orbits.
In the continuation of the coexistence point P3, some bifurcation points of codimension-1 are detected in the Figure 3 with respect to the bifurcation parameter A. There is a branch point BP1 (transcritical bifurcation) at A=Atc1=0.373911 around the equilibrium point P2=(0.293797,0.521937,0). One Hopf point is obtained at A=AH=0.388133 in the interior R3+ around the equilibrium point P3=(0.301964,0.370331,0.227896). For this Hopf point, the corresponding first Lyapunov coefficient is (−2.912374e−002)<0, indicating a supercritical Hopf bifurcation. Thus, a stable limit cycle exists, bifurcating from the equilibrium. Another branch point BP2 occurs at A=Atc2=0.714286 around the equilibrium point P1≈(0.4,0,3). Figure 3 shows that as the value of A increases close to the branch point BP2 around A=0.714, the solutions in the interior R3+ vanish and will appear in the nE-plane. The bifurcation diagram with respect to the parameter A for the system (3.2) is shown in Figure 4 for A in the interval (0.37,0.714). The complex dynamics are evident from this diagram. To confirm the complexity and chaotic dynamics of the system (3.2), the Lyapunov exponents and their corresponding dimensions are computed and presented in Figures 5 and 6.
The Lyapunov exponents for A=0.39 and A=0.4 are computed in Figure 5(a),(b), with their corresponding dimensions DL=2.4305 and DL=2.0952, respectively. Figure 6 shows the Lyapunov exponents for A=0.47 and A=0.6, with the computed dimensions DL=2.25545 and DL=2.4157, respectively. The presence of positive Lyapunov exponents confirms the complex dynamical behavior of the system for certain values of the parameter A.
The different dynamic behaviors observed in Figure 3 are confirmed by drawing phase portraits with respect to different values of A in Figures 7–9 and Figure 2. In particular, Figure 7A shows the local asymptotic stability of the boundary equilibrium point P2(¯n,¯p,0) for A=0.368. Due to the availability of some alternate food, the predator survives and also takes food from the prey. However, the harvesting effort diminishes to zero. If we increase A beyond BP1, at A=0.374, the interior point P3 is locally asymptotically stable (see Figure 7B). Considering A>AH, the system destabilizes, and strange attractors are obtained in Figures 8 and 9. It can be observed that as the value of A increases close to the branch point BP2 around A=0.714, the chaotic solution in the interior R3+ vanishes, and periodic solutions appear in the nE -plane. The analytical result in (3.3.2) regarding the stability of the point (ˆn,0,ˆE) also confirms the periodic solutions in the nE -plane. The change in the behavior of the solution can be seen in Figure 2.
Further, the bifurcation of codimension-1 are also explored with respect to the mortality rate of the predators μ in Figure 10. The branch point BP1 at μ=1.983924 around P2≈(0.293843,0.520833,0) is computed. The Hopf points H1 and H2 are detected, where the Hopf point H1 at μ=2.019724 around P3≈(0.301694,0.374208,0.221869) and the Hopf point H2 at μ=2.468001 around P1≈(0.4,0,3.000001) are obtained. For these Hopf points, the corresponding first Lyapunov coefficients are computed as (−2.964738×10−2)<0 and 8.841215×10−4>0, respectively. Accordingly, H1 is a supercritical and H2 is a subcritical Hopf bifurcation point.
Next, the numerical simulations of the system (3.2) are discussed when the migration rate is not dependent on the predator density (i.e., ˆk(p) = constant). In this case, the functions r(p), q(p), and f(p) will assume constant values. Let these values be chosen as follows:
r=5,q=2,v1=0.4 | (3.20) |
The dynamic behavior of this system using the datasets in (3.19) and (3.20) is illustrated in Figures 11–13. These figures shows the complex dynamics even with constant migration rates. This complexity occurs due to effort dynamics. Figure 14 considers the case for A=0, i.e., when the predators have no interaction with the prey and it depend only on alternative food. This shows that predators can sustain themselves even in the absence of prey because of availability of alternate food resources. Also, for the case A=1, (i.e., when predators have no dependence on alternate food and the predators interact with prey species only). Due to the refuges and harvesting of prey, there is a scarcity of food (prey) for predators. Therefore, this can lead to the extinction of predators, and this can be seen in Figure 15.
The different bifurcation diagrams for this system are drawn in Figure 16 with respect to the bifurcation parameter A in the interval (0.3357,0.35) for different values of the catchability coefficient q. In Figure 16A–D, different kinds of complexity can be observed for different values of q.
Previously, the model has been studied with linear harvesting, where chaotic behavior is observed within the system. To control this, the model will be analyzed with nonlinear harvesting, which may provide a more effective approach to regulating these dynamics. The following analysis will examine how nonlinear harvesting influences the stability of the system.
The dynamics of the corresponding complete predator–prey system, incorporating nonlinear harvesting and an alternate food source for the predators are illustrated in this section. Here, the harvesting rate depends on the prey's population size. Thus, the harvesting functions in (2.1) can be taken as H1(n1,n2,E)=q1En11+m1n1+m2n2 and H2(n1,n2,E)=q2En21+m1n1+m2n2, and the model is described as follows:
dn1dτ=(kn2−ˆk(p)n1)+ε(r1n1−Aan1p−q1En11+m1n1+m2n2)dn2dτ=(ˆk(p)n1−kn2)+ε(−r2n2−q2En21+m1n1+m2n2)dpdτ=ε(−μp+Abn1p+β(1−A)p)dEdτ=ε(−cE+p0q1En11+m1n1+m2n2+p0q2En21+m1n1+m2n2)} | (4.1) |
The nonlinear terms p0q1En11+m1n1+m2n2 and p0q2En21+m1n1+m2n2 represent the revenue or benefit from harvesting in Layer-1 and Layer-2, respectively, with the denominator capturing the effect of ecological saturation, while m1 and m2 are positive constants. The last equation means that if the earnings from harvesting are higher than costs, more effort is put into harvesting, using additional boats or nets. If the costs are too high, the harvesting effort decreases because harvesting is no longer profitable.
The aggregated model for the system mentioned above is as follows:
dndt=n(r(p)−Aaf(p)p−q(p)E1+m(p)n)=n⋅F(n,p,E)dpdt=p(−μ+Abf(p)n+β(1−A))=p⋅G(n,p)dEdt=E(−c+p0q(p)n1+m(p)n)=E⋅H(n,p)} | (4.2) |
n(0)>0,p(0)>0,E(0)>0, |
where m(p)=m1k+m2αpk+αp.
The system described by (4.2) possesses four positive steady states, enumerated as follows:
1) The trivial fixed point ~P0(0,0,0) is exists and the boundary equilibrium point ~P2(¯n,¯p,0) is feasible for the condition (3.3).
2) The predator-free boundary equilibrium point ~P1(ˆn,0,ˆE)=(cp0q1−cm1,0,r1p0p0q1−cm1) exists in the positive nE-plane for the condition p0q1>cm1.
3) The unique interior equilibrium point ~P3(n∗,p∗,E∗) of (4.2) exists, where n∗,p∗, and E∗ are obtained as follows:
n∗=c(k+αp∗)p0(q1k+q2αp∗)−c(m1k+m2αp∗)p∗=k(−bc+β(p0q1−cm1))(~A2−A)αβ(p0q2−cm2)(A−A0);~A2=β−μβ−bcp0q1−cm1;~A2>AE∗=(r1k−p∗(r2α+Aak))p0p0(q1k+q2αp∗)−c(m1k+m2αp∗) |
It can be observed that
1−μβ<β−μβ−bcp0q1−cm1 |
Accordingly, n∗ and p∗ are positive for
A0<A<~A2 | (4.3) |
The value of E∗ is positive for
p∗<r1kr2α+Aak<r1kr2α⇒Abck−p0q1k(μ−β(1−A))+cm1k(μ−β(1−A))α(p0q2−cm2)(μ−β(1−A))<r1kr2α |
which gives
A>β−μβ−T0(=~A1);T0=bcr2r2(p0q1−cm1)+r1(p0q2−cm2) | (4.4) |
Since A0<~A1, the point ~P3(n∗,p∗,E∗) is feasible for the following condition:
β−μβ−T0<A<β−μβ−bcp0q1−cm1 | (4.5) |
That is,
~A1<A<~A2 |
The local stability conditions for feasible equilibrium points of the system described by the equations in (4.2) are achieved by assessing the characteristics of the eigenvalues derived from the Jacobian matrix calculated at the corresponding equilibrium points.
The Jacobian matrix corresponding to the system (4.2) at (n,p,E) is expressed as:
˜J(n,p,E)=[q(p)m(p)En(1+m(p)n)2+Fn(r′(p)−Aa(f(p)+f′(p)p)−q′(p)E1+m(p)n+q(p)m′(p)En(1+m(p)n)2)−nq(p)1+m(p)nAbf(p)pAbf′(p)pn+G0p0q(p)E(1+m(p)n)2p0q′(p)En1+m(p)n−p0q(p)m′(p)En2(1+m(p)n)2H] |
For the value of f′(p),r′(p) and q′(p), refer to (3.8). Moreover,
m′(p)=αk(m2−m1)(k+αp)2>0;m2>m1 |
This is as discussed in Subsubsection (3.3.1).
At the point ~P1(ˆn,0,ˆE), the Jacobian matrix assumes the following configuration:
~J1(ˆn,0,ˆE)=[a11a12a130a220a31a320]=[q1m1ˆEˆn(1+m1ˆn)2ˆn(−α(r2+r1)k−Aa−ˆEα(q2−q1)k(1+m1ˆn)+ˆEˆnq1α(m2−m1)k(1+m1ˆn)2)−ˆnq11+m1ˆn0−μ+Abˆn+β(1−A)0ˆEp0q1(1+m1ˆn)2p0ˆEˆnα(q2−q1)k(1+m1ˆn)−p0ˆEˆn2q1α(m2−m1)k(1+m1ˆn)20]=[m1r1cp0q1cp0q1−cm1(−α(r2+r1)k−Aa−α(q2−q1)r1kq1+α(m2−m1)r1ckp0q1)−cp00−μ+Abcp0q1−cm1+β(1−A)0r1(p0q1−cm1)2q1p0cr1α(q2−q1)−r1α(m2−m1)c2k(p0q1−cm1)q10] |
The characteristic equation of the matrix ~J1 about the point ~P1 is given by
(λ+cp0)(λ2−a32λ−a31a22)=0 |
One eigenvalue is λ11=−cp0<0 and the remaining two eigenvalues can be obtained from the characteristic equation
λ2−a32λ−a31a22=0 | (4.6) |
The stability of the point ~P1 can be obtained by evaluating trace(~J1)=a32<0 and det(~J1)=a31a22>0 from Eq (4.6) for the following condition:
p0<m2−m1q2−q1andAbc+βp0q1+cm1(βA+μ)>βcm1+p0q1(βA+μ) | (4.7) |
Therefore, the point ~P1 is locally asymptotically stable. However, if the condition (4.7) is violated, then the point ~P1 becomes a saddle point. Therefore, bifurcation may occur around ~P1 provided that
p0=m2−m1q2−q1andAbc+βp0q1+cm1(βA+μ)=βcm1+p0q1(βA+μ) | (4.8) |
The Jacobian matrix computed at the point ~P2(¯n,¯p,0) is given by
~J2(¯n,¯p,0)=[0b12b13b21b22000b33]=[0¯n(r′(¯p)−Aa(f(¯p)+f′(¯p)¯p))−¯nq(¯p)1+m(¯p)¯nAbf(¯p)¯pAbf′(¯p)¯p¯n000−c+p0q(¯p)¯n1+m(¯p)¯n] |
The characteristic equation associated with the matrix ~J2(¯n,¯p,0) is given by
(λ−(−c+p0q(¯p)¯n1+m(¯p)¯n))(λ2−(Abf′(¯p)¯p¯n))λ−Abf(¯p)¯p¯n(r′(¯p)−Aa(f(¯p)+f′(¯p)¯p)))=0 |
One eigenvalue is λ21=−c+p0q(¯p)¯n1+m(¯p)¯n and the remaining two eigenvalues can be derived from the following characteristic equation:
λ2−(Abf′(¯p)¯p¯n))λ−Abf(¯p)¯p¯n(r′(¯p)−Aa(f(¯p)+f′(¯p)¯p))=0 | (4.9) |
The trace(~J2) and det(~J2) of Eq (4.9) are same as those mentioned in Subsubsection (3.3.3). Due to existence condition for the point ~P2, it can be observed that trace(~J2)<0 and det(~J2)>0. Accordingly we have the following:
1) The point ~P2(¯n,¯p,0) is locally asymptotically stable for λ21<0, which gives the condition
¯n>−c(k+α¯p)c(m1k+m2α¯p)+p0(q1k+q2α¯p) | (4.10) |
2) Moreover, the point ~P2(¯n,¯p,0) transforms into a saddle point when (4.10) is violated. Therefore, bifurcation may occur around ~P2(¯n,¯p,0) for
¯n=−c(k+α¯p)c(m1k+m2α¯p)+p0(q1k+q2α¯p) | (4.11) |
The Jacobian matrix for the interior point ~P3(n∗,p∗,E∗) is as follows:
~J3(n∗,p∗,E∗)=[q(p∗)m(p∗)E∗n∗(1+m(p∗)n∗)2n∗(r′(p∗)−Aa(f(p∗)+f′(p∗)p∗)−q′(p∗)E∗1+m(p∗)n∗+q(p∗)m′(p∗)E∗n∗(1+m(p∗)n∗)2)−n∗q(p∗)1+m(p∗)n∗Abf(p∗)p∗Abf′(p∗)p∗n∗0p0q(p∗)E∗(1+m(p∗)n∗)2p0q′(p∗)E∗n∗1+m(p∗)n∗−p0q(p∗)m′(p∗)E∗n∗2(1+m(p∗)n∗)20] |
The characteristic equation for the given Jacobian matrix centered at the point (n∗,p∗,E∗) can be expressed as
λ3+B11λ2+B12λ+B13=0 | (4.12) |
with
B11=Abαkn∗p∗(k+αp∗)2−(q1k+q2αp∗)(m1k+m2αp∗)E∗n∗(k+αp∗+(m1k+m2αp∗)n∗)2B12=Abf′(p∗)q(p∗)m(p∗)n∗2p∗E∗−Abf(p∗)q(p∗)m′(p∗)n∗2p∗E∗(1+m(p∗)n∗)2+Abf(p∗)q′(p∗)n∗p∗E∗1+m(p∗)n∗+p0q(p∗)2n∗E∗(1+m(p∗)n∗)3−Abf(p∗)r′(p∗)n∗p∗+A2abf(p∗)n∗p∗(f(p∗)+f′(p∗)p∗)B13=Aαkbp0q(p∗)n∗2p∗E∗(k+αp∗)2(1+m(p∗)n∗)3(q(p∗)+f(p∗)(1+m(p∗)n∗)(q2−q1)−n∗q(p∗)f(p∗)(m2−m1)) |
By applying the Routh–Hurwitz criterion, the interior equilibrium point ~P3(n∗,p∗,E∗) is locally asymptotically stable if the following conditions are met:
B11>0,B13>0andB11B12>B13 |
Example 4.1. In general, to determine the stability of the point ~P3, it is difficult to analyze the signs B1,B2, and B3 to verify all the conditions above. Therefore, the stability analysis of the point ~P3(n∗,p∗,E∗) is conducted for the particular choice of data, namely k=0.8,r1=2,r2=1,a=2,α=3.5,q1=1,q2=1.5,μ=4,b=4,β=5,c=0.03,p0=0.1,m1=0.6,m2=0.8, and A=0.22. It can be observed that the interior point (n∗,p∗,E∗) occurs at (0.2778,0.3301,0.1555). It also satisfies all the conditions of the Routh– Hurwitz criterion such that B11=0.0312>0, B13=0.00028961>0, and B11B12=0.0312∗0.0789=0.0025>B13. Consequently, the interior point ~P3 is locally asymptotically stable.
Theorem 4.1. The system (4.2) undergoes a transcritical bifurcation around the planar equilibrium point ~P1(ˆn,0,ˆE) as the bifurcation parameter A crosses a critical value such that
A=˜Atc1 | (4.13) |
Proof. The Jacobian matrix ~J1 of the system (4.2) at the equilibrium point ~P1 has a zero eigenvalue for the condition (4.13), then the matrix ~J1 becomes ˜J∗1 (which indicates a matrix when the condition (4.13) is satisfied, i.e., a22=0). Therefore, the Jacobian matrix ˜J∗1 and the transpose of the Jacobian matrix ˜J∗T1 have the following eigenvectors corresponding to a zero eigenvalue:
˜Ve1=(v1,v2,v3=1)Tand˜We1=(w1,w2,w3)T=(0,1,0)T, |
where
v1=q1(p0cα(q2−q1)−α(m2−m1)c2)m1r1(p0cα(q2−q1)−α(m2−m1)c2)+(cm1−p0q1)2(α(r2+r1)p0q1+Aakp0q1+α(q2−q1)r1p0−α(m2−m1)r1c)v2=kq1(cm1−p0q1)3m1r1(p0cα(q2−q1)−α(m2−m1)c2)+(cm1−p0q1)2(α(r2+r1)p0q1+Aakp0q1+α(q2−q1)r1p0−α(m2−m1)r1c) |
The computation for the conditions of Sotomayar's theorem is stated as follows:
Δ11=˜WTe1ΨA(~P1,˜Atc1)=0,Ψ=(nF,pG,EH)TΔ12=˜WTe1[DΨA(~P1,˜Atc1)˜Ve1]=(bcp0q1−cm1−β)v2≠0Δ13=˜WTe1[D2ΨA(˜P1,˜Atc1)(˜Ve1,˜Ve1)]=2Abk2(k+αp)2v1v2≠0 |
Since, all the conditions required by Sotomayar's theorem for a transcritical bifurcation are satisfied, the system (4.2) undergoes a transcritical bifurcation around the point ~P1.
Example 4.2. For a particular dataset k=0.8,r1=2,r2=1,a=2,α=3.5,q1=1,q2=1.5,μ=4,b=4,β=5,c=0.03,p0=0.1,m1=0.6, and m2=0.8, the Jacobian matrix ˜J1 of ~P1(0.365854,0,2.439024) is confirmed to have a zero eigenvalue at ˜Atc1=0.282759. The Jacobian matrix ˜J∗1 is given below:
˜J∗1(~P1,˜Atc1)=[0.3600−6.4173−0.30000.01340.14090] | (4.14) |
The eigenvectors corresponding to the zero eigenvalue (i.e., a22=0) are ˜Ve1=(−0.2946,0.0281,−0.9552)T and ˜We1=(0,1,0)T. This analysis confirms that all the conditions of Sotomayar's theorem for a transcritical bifurcation at the point ~P1 are verified such that Δ11=0, Δ12≠0, and Δ13≠0.
Theorem 4.2. The system (4.2) admits a transcrtical bifurcation around the point ~P2(¯n,¯p,0) when the bifurcation parameter A varies such that
A=˜Atc2 | (4.15) |
Proof. The Jacobian matrix ~J2(¯n,¯p,0) of the system (4.2) at the equilibrium point ~P2 has a zero eigenvalue for the condition (4.15), then the matrix ~J2 becomes ˜J∗2 (which indicates a matrix when the condition (4.15) is satisfied, i.e., b33=0). Therefore, the Jacobian matrix ˜J∗2 and the transpose of the Jacobian matrix ˜J∗T2 have the following eigenvectors corresponding to the zero eigenvalue:
˜Ve2=(v4,v5,v6=1)Tand˜We2=(w4,w5,w6)T=(0,0,1)T, |
where,
v4=−α¯n(q1k+q2α¯p)(k+α¯p)(k+α¯p+(m1k+m2α¯p)¯n)(αk(r2+r1)+Aak2)v5=−(q1k+q2α¯p)(k+α¯p)2(k+α¯p+(m1k+m2α¯p)¯n)(αk(r2+r1)+Aak2) |
The computation for the conditions of Sotomayar's theorem is stated as follows:
Δ21=˜WTe2ΨA(~P2,˜Atc2)=0,Ψ=(nF,pG,EH)TΔ22=˜WTe2[DΨA(~P2,˜Atc2)˜Ve2]=db33dAv6≠0Δ23=˜WTe2[D2ΨA(~P2,˜Atc2)(˜Ve2,˜Ve2)]=2((k+α¯p+m1k¯n+m2α¯p¯n)(p0q1k+p0q2α¯p)−¯np0(q1k+q2α¯p)(m1k+m2α¯p)(k+α¯p+m1k¯n+m2α¯p¯n)2)v4v6+2((k+α¯p+m1k¯n+m2α¯p¯n)(q2α¯np0)−¯np0(q1k+q2α¯p)(α+m2α¯n)(k+α¯p+m1k¯n+m2α¯p¯n)2)v5v6≠0 |
Since, all the conditions required by Sotomayar's theorem for a transcritical bifurcation are satisfied, the system (4.2) undergoes a transcritical bifurcation around the point ~P2(¯n,¯p,0).
Example 4.3. For a particular dataset k=0.8,r1=2,r2=1,a=2,α=3.5,q1=1,q2=1.5,μ=4,b=4,β=5,c=0.03,p0=0.1,m1=0.6, and m2=0.8, the Jacobian matrix ~J2 of ~P2(0.271750,0.415937,0) is confirmed to have a zero eigenvalue at ˜Atc2=0.216708. The Jacobian matrix ˜J∗2 is given below:
˜J∗2(~P2,˜Atc2)=[0−0.4634−0.30.1279−0.05390000] | (4.16) |
The eigenvectors corresponding to the zero eigenvalue (i.e., b33=0) are ˜Ve2=(−0.2233,−0.5297,0.8182)T and ˜We2=(0,0,1)T. This analysis confirms that all the conditions of Sotomayar's theorem for a transcritical bifurcation at the point ~P2 are verified such that Δ21=0, Δ22≠0, and Δ23≠0.
Example 4.4. For the parameter values k=0.8,r1=2,r2=1,a=2,α=3.5,q1=1,q2=1.5,μ=4,b=4,β=5,c=0.03,p0=0.1,m1=0.6, and m2=0.8, the system (4.2) exhibits Hopf bifurcation around the interior equilibrium point (0.284223,0.263246,0.321669) when the critical value of A=˜AH=0.223632 are met. The conditions for the Liu's criterion at the bifurcation parameter ˜AH are as follows:
ξ1(˜AH)=0.0064>0andξ3(˜AH)=0.00064>0Φ(˜AH)=ξ1(˜AH)∗ξ2(˜AH)−ξ3(˜AH)=0.0064∗0.0995−0.00064=0anddΦ(˜AH)d˜AH=0.0280≠0 |
Consider the following choice of hypothetical data for the system (4.2) with nonlinear harvesting in appropriate units:
k=0.8,r1=2,r2=1,a=2,α=3.5,q1=1,q2=1.5,μ=4,b=4,β=5,c=0.03,p0=0.1,m1=0.6,m2=0.8 | (4.17) |
The critical value of A are computed as follows:
A0=0.2,˜A1=0.2155,˜A2=0.2828 |
The aggregated system reaches its equilibrium boundary at ~P1≈(0.3659,0,2.4390). A is chosen as A=0.21 (A0<A<˜A1), then another equilibrium point on the boundary, denoted as ˜P2(0.1681,0.4171,0), is identified and presented in Figure 19. Following that, with A=0.22, an interior equilibrium point is obtained ˜P3(0.2778,0.3301,0.1555) as shown in Figure 20. (see the example in (4)).
Following the interior point ˜P3, several codimension-1 bifurcation points are identified in Figure 17 with respect to the bifurcation parameter A. These are as follows:
1) One such branch point BP1 (also known as transcritical bifurcation), occurs at A=˜Atc1=0.216708 around ~P2=(0.271750,0.415937,0).
2) Additionally, a Hopf point emerges at A=˜AH=0.223632 within the interior R3+ around the point ~P3=(0.284223,0.263246,0.321669). This Hopf point has a first Lyapunov coefficient of (−3.835306e−01)<0, indicating a supercritical Hopf bifurcation and implying the presence of a stable limit cycle.
3) Another branch point BP2 is obtained at A=˜Atc2=0.282759 around the point ~P1≈(0.365854,0,2.439024).
Figure 18 illustrates the bifurcation diagram with respect to the parameter A for the system (4.2) in the interval (0.21,0.3). Figures 19–23 represent the phase portraits and time series plots of the system (4.2) and are generated with different values of A to validate the diverse dynamic behaviors of the system, which have been observed in Figure 17. For instance, Figure 19 highlights the locally asymptotically stable boundary equilibrium point around ~P2(0.1681,0.4171,0) at A=0.21. Meanwhile, for A=0.22, Figure 20 showcases the locally asymptotically stable behavior around the point ~P3(0.2778,0.3301,0.1555). When A>˜AH=0.223632, the system becomes destabilized, resulting in the appearance of strange attractors, as illustrated in Figures 21–23.
This article examined the dynamic interactions between predators and prey in an ecosystem where the prey have access to a refuge and there is an alternate food source for the predators. The model consists of two layers: Layer 1 for predator dynamics and Layer 2 for prey population dynamics. While the prey are capable of migrating between both layers, the predators remain confined to Layer 1. The migration rate of the prey is assumed to be density-dependent on the predator population. Additionally, the model incorporated the presence of alternate food sources for the predators, which support their survival when the availability of prey is reduced due to the refuge or harvesting. Both linear and nonlinear harvesting strategies were applied to the prey population in both layers, with effort dynamics considered in each case. Two distinct timescales are introduced: a fast scale describing the inter-layer migration and a slow scale governing the growth of the population, the predator–prey interactions, and harvesting efforts. By exploiting these timescales, aggregation techniques are employed to reduce the dimensionality of the system. This study highlights the significant impact of both linear and nonlinear harvesting effort on the dynamic behavior of the predator–prey system.
The initial analysis focuses on the linear harvesting dynamics in the predator–prey system described by (3.2). Four distinct feasible equilibrium points are computed: the trivial equilibrium point (P0), the predator-free boundary equilibrium (P1), the boundary equilibrium in the absence of harvesting (P2), and the interior equilibrium point (P3). The dynamic behavior of the aggregated system is then discussed in detail. It is found that with A acting as the bifurcation parameter, the system displays two transcritical bifurcation points and one Hopf bifurcation point. The system shows local asymptotic stability at the equilibrium points P2 and P3 when A<Atc1 and Atc1<A<AH, respectively. For values of A within the range AH<A<Atc2, the system exhibits chaotic behavior. Beyond Atc2, periodic solutions are observed. Additionally, the dynamics of the Lyapunov exponents and their corresponding dimensions are calculated, confirming the chaotic dynamics of the system. The presence of positive Lyapunov exponents indicate the complex dynamical behavior of the system for certain ranges of the parameter A.
The study of the dynamic behavior of predator–prey system with linear harvesting can lead to chaotic behavior in the system. To address this, the dynamics of the system with nonlinear harvesting (4.2) are examined. Theoretical analysis reveals that by choosing the additional food parameter A as the bifurcation parameter, the model undergoes a Hopf bifurcation near ~P3 (the interior equilibrium point), as well as two transcritical bifurcations around the equilibrium points ~P1 (the predator-free equilibrium point) and ~P2 (the harvesting-free equilibrium point). Numerical simulations are conducted for various values of A to validate the theoretical results. The simulations indicated that the system reaches the point ~P2 point when A0<A<˜A1, which implies no harvesting of prey species is possible due to an insufficient amount of prey, and the system stabilizes at ~P3 for the range ˜A1<A<˜A2. The coexistence of prey and predator populations occurs within this suitable range of additional food (˜A1<A<˜A2). Moreover, the study demonstrated that the system exhibits local asymptotic stability at ~P2 when A<˜Atc1 and at ~P3 when ˜Atc1<A<˜AH.
The study demonstrated that linear harvesting causes chaotic behavior in the system. This means that population levels become unpredictable, fluctuating in an irregular way, which can lead to instability or even extinction in some cases. However, nonlinear harvesting maintains stability. In biological terms, this is similar to natural population control mechanisms, such as predators consuming more prey when prey populations are high and consuming less when prey populations are low. This kind of adaptive response helps regulate the ecosystem, preventing extreme fluctuations. Therefore, the study concludes that nonlinear harvesting strategies can effectively regulate chaos, ensuring long-term population stability. Ecologically, these findings highlight the necessity of integrating nonlinear harvesting approaches into sustainable harvesting practices to prevent population declines and preserve ecosystem balance.
The authors declare they have not used artificial intelligence (AI) tools in the creation of this article.
We sincerely thank all the reviewers for their insightful comments, which have greatly contributed to improving our research and manuscript. The authors would like to thank the management of Vellore Institute of Technology (VIT) for providing the facility to carry out this study.
The authors declare there is no conflict of interest.
[1] |
P. K. Santra, G. S. Mahapatra, G. R. Phaijoo, Bifurcation analysis and chaos control of discrete prey–predator model incorporating novel prey–refuge concept, Comput. Math. Methods, 3 (2021), e1185. https://doi.org/10.3934/mbe.2022313 doi: 10.3934/mbe.2022313
![]() |
[2] |
H. Molla, Sabiar Rahman, S. Sarwardi, Dynamics of a predator–prey model with Holling type Ⅱ functional response incorporating a prey refuge depending on both the species, Int. J. Nonlinear Sci. Numer. Simul., 20 (2019), 89–104. https://doi.org/10.1515/ijnsns-2017-0224 doi: 10.1515/ijnsns-2017-0224
![]() |
[3] |
T. K. Kar, Stability analysis of a prey–predator model incorporating a prey refuge, Commun. Nonlinear Sci. Numer. Simul., 10 (2005), 681–691. https://doi.org/10.1016/j.cnsns.2003.08.006 doi: 10.1016/j.cnsns.2003.08.006
![]() |
[4] |
S. Pandey, U. Ghosh, D. Das, S. Chakraborty, A. Sarkar, Rich dynamics of a delay-induced stage-structure prey-predator model with cooperative behaviour in both species and the impact of prey refuge, Math. Comput. Simul., 216 (2024), 49–76. https://doi.org/10.1016/j.matcom.2023.09.002 doi: 10.1016/j.matcom.2023.09.002
![]() |
[5] |
W. Lu, Y. Xia, Y. Bai, Periodic solution of a stage-structured predator-prey model incorporating prey refuge, Math. Biosci. Eng., 17 (2020), 3160–3174. https://doi.org/10.3934/mbe.2020179 doi: 10.3934/mbe.2020179
![]() |
[6] |
J. Ghosh, B. Sahoo, S.Poria, Prey-predator dynamics with prey refuge providing additional food to predator, Chaos, Solitons Fractals, 96 (2017), 110–119. https://doi.org/10.1016/j.chaos.2017.01.010 doi: 10.1016/j.chaos.2017.01.010
![]() |
[7] |
P. D. N. Srinivasu, B. S. R. V. Prasad, M. Venkatesulu, Biological control through provision of additional food to predators: a theoretical study, Theor. Popul Biol., 72 (2007), 111–120. https://doi.org/10.1016/j.tpb.2007.03.011 doi: 10.1016/j.tpb.2007.03.011
![]() |
[8] |
S. Samaddar, M. Dhar, P. Bhattacharya, Effect of fear on prey–predator dynamics: Exploring the role of prey refuge and additional food, Chaos: Interdiscip. J. Nonlinear Sci., 30 (2020). https://doi.org/10.1063/5.0006968 doi: 10.1063/5.0006968
![]() |
[9] |
R. Rani, S. Gakkhar, The impact of provision of additional food to predator in predator–prey model with combined harvesting in the presence of toxicity, J. Appl. Math. Comput., 60 (2019), 673–701. https://doi.org/10.1007/s12190-018-01232-z doi: 10.1007/s12190-018-01232-z
![]() |
[10] |
D. K. Jana, P. Panja, Effects of supplying additional food for a scavenger species in a prey-predator-scavenger model with quadratic harvesting, Int. J. Modell. Simul., 43 (2022), 250–264. https://doi.org/10.1080/02286203.2022.2065658 doi: 10.1080/02286203.2022.2065658
![]() |
[11] |
P. Panja, S. Jana, S. Kumar Mondal, Effects of additional food on the dynamics of a three species food chain model incorporating refuge and harvesting, Int. J. Nonlinear Sci. Numer. Simul., 20 (2019), 787–801. https://doi.org/10.1515/ijnsns-2018-0313 doi: 10.1515/ijnsns-2018-0313
![]() |
[12] |
M. Kaur, R. Rani, R. Bhatia, G. N. Verma, S. Ahiwar, Dynamical study of quadrating harvesting of a predator–prey model with Monod–Haldane functional response, J. Appl. Math. Comput., 66 (2021), 397–422. https://doi.org/10.1007/s12190-020-01438-0 doi: 10.1007/s12190-020-01438-0
![]() |
[13] |
S. R. Jawad, M. Al Nuaimi, Persistence and bifurcation analysis among four species interactions with the influence of competition, predation and harvesting, Iraqi J. Sci., (2023), 1369–1390. https://doi.org/10.24996/ijs.2023.64.3.30 doi: 10.24996/ijs.2023.64.3.30
![]() |
[14] |
S. Tolcha, B. K. Bole, P. R. Koya, Population dynamics of two mutuality preys and one predator with harvesting of one prey and allowing alternative food source to predator, Math. Model. Appl., 5 (2020), 55. https://doi.org/10.11648/j.mma.20200502.12 doi: 10.11648/j.mma.20200502.12
![]() |
[15] | D. S. Al-Jaf, The role of linear type of harvesting on two competitive species interaction, Commun. Math. Biol. Neurosci., (2024). lhttps://doi.org/10.28919/cmbn/8426 |
[16] |
Y. Wang, X. Zhou, W. Jiang, Bifurcations in a diffusive predator-prey system with linear harvesting, Chaos, Solitons Fractals, 169 (2023), 113286. https://doi.org/10.1016/j.chaos.2023.113286 doi: 10.1016/j.chaos.2023.113286
![]() |
[17] | M. Chen, W. Yang, Effect of non-linear harvesting and delay on a predator-prey model with Beddington-DeAngelis functional response, Commun. Math. Biol. Neurosci., 2024. https://doi.org/10.28919/cmbn/8382 |
[18] |
S. Pal, S. Majhi, S. Mandal, N. Pal, Role of fear in a predator-prey model with Beddington–DeAngelis functional response, Z. Naturforsch. A, 74 (2019), 581–595. https://doi.org/10.1515/zna-2018-0449 doi: 10.1515/zna-2018-0449
![]() |
[19] |
R. Rani, S. Gakkhar, A. Moussaoui, Dynamics of a fishery system in a patchy environment with nonlinear harvesting, Math. Methods Appl. Sci., 42 (2019), 7192–7209. https://doi.org/10.1002/mma.5826 doi: 10.1002/mma.5826
![]() |
[20] | R. Rani, S. Gakkhar, Non-linear effort dynamics for harvesting in a Predator- Prey system, J. Nat. Sci. Res., 5 (2015). |
[21] |
M. Li, B. Chen, H. Ye, A bioeconomic differential algebraic predator–prey model with nonlinear prey harvesting, Appl. Math. Modell., 42 (2017), 17–28. https://doi.org/10.1016/j.apm.2016.09.029 doi: 10.1016/j.apm.2016.09.029
![]() |
[22] |
J. C. Poggiale, P. Auger, Impact of spatial heterogeneity on a predator–prey system dynamics, C.R. Biol., 327 (2004), 1058–1063. https://doi.org/10.1016/j.crvi.2004.06.006 doi: 10.1016/j.crvi.2004.06.006
![]() |
[23] | S. A. Levin, Population dynamic models in heterogeneous environments, Ann. Rev. Ecol. Syst., 7 (1976), 287–310. |
[24] | K. Dao Duc, P. Augera, T. Nguyen-Huu, Predator density-dependent prey dispersal in a patchy environment with a refuge for the prey: biological modelling, S. Afr. J. Sci., 104 (2008), 180–184. |
[25] |
R. Mchich, P. Auger, J. C. Poggiale, Effect of predator density dependent dispersal of prey on stability of a predator-prey system, Math. Biosci., 206 (2007), 343–356. https://doi.org/10.1016/j.mbs.2005.11.005 doi: 10.1016/j.mbs.2005.11.005
![]() |
[26] |
A. El Abdllaoui, P. Auger, B.W. Kooi, R. B. De la Parra, R. Mchich, Effects of density-dependent migrations on stability of a two-patch predator–prey model, Math. Biosci., 210 (2007), 335–354. https://doi.org/10.1016/j.mbs.2007.03.002 doi: 10.1016/j.mbs.2007.03.002
![]() |
[27] |
A. Moussaoui, P. Auger, Simple fishery and marine reserve models to study the SLOSS problem, ESAIM: Proc. Surv., 49 (2015), 78–90. https://doi.org/10.1051/proc/201549007 doi: 10.1051/proc/201549007
![]() |
[28] |
M. Bensenane, A. Moussaoui, P. Auger, On the optimal size of marine reserves, Acta Biotheor., 61 (2013), 109–118. https://doi.org/10.1007/s10441-013-9173-9 doi: 10.1007/s10441-013-9173-9
![]() |
[29] | A. Moussaoui, M. Bensenane, P. Auger, A. Bah, On the optimal size and number of reserves in a multi-site fishery model, J. Biol. Syst., 22 (2014), 1–17. |
[30] | J. Carr, Applications of Centre Manifold Theory, Springer, 1981. https://doi.org/10.1007/978-1-4612-5929-9 |
[31] | P. Auger, R. B. De la Parra, J. C. Poggiale, E. Sánchez, T. Nguyen-Huu, Aggregation of variables and applications to population dynamics, Struct. Popul. Models Biol. Epidemiol., (2008), 209–263. https://doi.org/10.1007/978-3-540-78273-55 |
[32] |
T. Grozdanovski, J. J. Shepherd, A. Stacey, Multi-scaling analysis of a logistic model with slowly varying coefficients, Appl. Math. Lett., 22 (2009), 1091–1095. https://doi.org/10.1016/j.aml.2008.10.002 doi: 10.1016/j.aml.2008.10.002
![]() |
[33] | M. H. Holmes, Introduction to Perturbation Methods, Springer, 1995. https://doi.org/10.1002/zamm.19960760709 |
[34] |
P. Auger, S. Charles, M. Viala, J. C. Poggiale, Aggregation and emergence in ecological modelling: integration of ecological levels, Ecol. Modell., 127 (2000), 11–20. https://doi.org/10.1016/S0304-3800(99)00201-X doi: 10.1016/S0304-3800(99)00201-X
![]() |
[35] |
W. M. Liu, Criterion of Hopf bifurcations without using eigenvalues, J. Math. Anal. Appl., 182 (1994), 250–256. https://doi.org/10.1006/jmaa.1994.1079 doi: 10.1006/jmaa.1994.1079
![]() |
[36] |
D. Sen, S. Ghorai, M. Banerjee, Complex dynamics of a three species prey-predator model with intraguild predation, Ecol. Complexity, 34 (2018), 9–22. https://doi.org/10.1016/j.ecocom.2018.02.002 doi: 10.1016/j.ecocom.2018.02.002
![]() |
[37] |
R. P. Gupta, M. Banerjee, P. Chandra, Period doubling cascades of prey-predator model with nonlinear harvesting and control of over exploitation through taxation, Commun. Nonlinear Sci. Numer. Simul., 19 (2014), 2382–2405. https://doi.org/10.1016/j.cnsns.2013.10.033 doi: 10.1016/j.cnsns.2013.10.033
![]() |
[38] |
W. Govaerts, Y. A. Kuznetsov, A. Dhooge, Numerical continuation of bifurcations of limit cycles in MATLAB, SIAM J. Sci. Comput., 27 (2005), 231–252. https://doi.org/10.1137/030600746 doi: 10.1137/030600746
![]() |
[39] | A. Riet, A Continuation Toolbox in Matlab, Master's thesis, Mathematical Institute, Utrecht University, Utrecht, The Netherlands, 2000. |
[40] | W. Mestrom, Continuation of Limit Cycles in MATLAB, Master's thesis, Mathematical Institute, Utrecht University, Utrecht, The Netherlands, 2002. |
[41] |
A. Dhooge, W. Govaerts, Y. A. Kuznetsov, MATCONT: A MATLAB package for numerical bifurcation analysis of ODEs, ACM Trans. Math. Software (TOMS), 90 (2003), 141–164. https://doi.org/10.1145/779359.779362 doi: 10.1145/779359.779362
![]() |