Predator-prey dynamics with refuge, alternate food, and harvesting strategies in a patchy habitat

1.   Introduction

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.

2.   The complete (slow–fast) system

In a heterogeneous patchy environment, the density of the prey in Layer-$i$ $(i = 1, 2)$ is denoted by $n_i(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 $\tau$, which governs the movement of species between layers. To address the two timescales, a small dimensionless parameter $\varepsilon$ is introduced, and the slow timescale is defined as $t = \varepsilon \tau$, where $(0 < \varepsilon < < 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:

The constant $r_1 > 0$ indicates the prey's intrinsic growth rate in the surface layer. The term $r_1n_1$ represents the growth rate of the prey in Layer-1, which is assumed to be follow exponential growth, while the term $-r_2n_2$ represents the mortality rate of the prey due to the lack of food in Layer 2. Additionally, there is a negative term $Aan_1p$ accounting for prey loss due to predation in Layer 1, where the predation rate is given by $a$ for the prey. The functions $H_1(n_1, n_2, E)$ and $H_2(n_1, n_2, E)$ represent the harvesting terms in Layer-1 and Layer-2, respectively. The constant $\mu$ is the natural mortality rate for the predators, while $b$ signifies the food conversion rate by the predators with respect to the prey $n_1(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 $\beta(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 $p_0$ 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 $p_0 H_1(n_1, n_2, E)$ and $p_0H_2(n_1, n_2, 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 $\hat{k}$ :

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.

2.1.   The fast system

To study the fast dispersal model, we neglect the slow part of the system (2.1) by taking $\varepsilon = 0$:

The fast equilibria for the fast system (2.2) are obtained as follows:

The equilibrium frequencies $v_1^{*}$ and $v_2^{*}$ take the following form:

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) = n_1(t) + n_2(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.

3.   Mathematical model 1 with linear harvesting

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 $H_1(n_1, n_2, E) = q_1En_1$ and $H_2(n_1.n_2, E) = q_2En_2$. It follows that

The parameters $q_1$ and $q_2$ represent the catchability coefficient in Layer-1 and Layer-2, respectively, in which $q_1En_1$ and $q_2En_2$ 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.

3.1.   The aggregated model

Applying the aggregation method, the complete system (3.1) is reduced to a system of three ordinary differential equations. Let $n(t) = n_1(t) + n_2(t)$ be the aggregated variable. Further, we introduce $r(p)$, $q(p)$, and $m(p)$ as given below:

The aggregated system presented below is obtained by summing the two prey equations and substituting the fast equilibrium into the complete system (3.1)

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.

3.2.   Existence of equilibrium points in the aggregated model (3.2)

The system (3.2) possesses four positive equilibrium points:

1) The trivial point $P_0(0, 0, 0)$ is always exists.

2) $P_1(\hat{n}, 0, \hat{E}) = P_1\bigg(\dfrac{c}{p_0q_1}, 0, \dfrac{r_1}{q_1} \bigg)$, the predator-free boundary's fixed point in the positive $nE$-plane.

3) The harvesting-free boundary equilibrium point

exists in the positive $np$-plane and is positive for the critical value of the alternate food resource as follows:

4) The unique interior fixed point $P_3(n^*, p^*, E^*)$ of (3.2) exists, where $n^*, p^*$, and $E^*$ are obtained as shown below:

The interior point $P_3(n^*, p^*, E^*)$ is feasible for the condition

where

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:

It can be observed that the equilibrium points $P_0$ and $P_1$ may exist irrespective of $A$. However, the boundary equilibrium point $P_2$ 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.

3.3.   Stability analysis of the equilibrium points of the model (3.2)

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:

where

3.3.1.   Stability of $P_0(0, 0, 0)$

The Jacobian matrix, when calculated at the equilibrium point $P_0(0, 0, 0)$, takes the following form:

The eigenvalues of $J_0$ are

Hence, the origin point $P_0(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:

If the condition in (3.9) is violated, then $P_0$ has an unstable manifold in the $p$ direction. Accordingly, the trajectories along $p = 0$, starting in neighborhood of $P_0$, may be attracted to $P_1$ when the condition in (3.9) is violated.

3.3.2.   Stability of $P_1(\hat{n}, 0, \hat{E})$

At the point $P_1(\hat{n}, 0, \hat{E})$, the Jacobian matrix is given by

The characteristic equation derived from the Jacobian matrix about the point $P_1(\hat{n}, 0, \hat{E})$ is provided as follows:

where $\quad C_{11} = -c_{22}, \quad C_{12} = -c_{13}c_{31}, \quad C_{13} = c_{13}c_{31}c_{22},$ and $C_{11}C_{12}-C_{13} = c_{22}c_{13}c_{31}-c_{13}c_{31}c_{22} = 0$

Since $C_{11}C_{12} = C_{13}$, Eq (3.10) may be rewritten as

This gives

This confirms the periodic solutions around the equilibrium point $P_1(\hat{n}, 0, \hat{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.

3.3.3.   Stability of $P_2(\overline{n}, \overline{p}, 0)$

The following Jacobian matrix is obtained for the point $P_2(\overline{n}, \overline{p}, 0)$

The characteristic equation associated with the matrix $J_2(\overline{n}, \overline{p}, 0)$ is given by

One eigenvalue is $\lambda_{21} = -c+p_0\overline{n}q(\overline{p})$, while the other two eigenvalues can be determined from the following characteristic equation,

The $trace(J_2)$ and $det(J_2)$ of Eq (3.11) are computed as follows:

Accordingly, the point $P_2(\overline{n}, \overline{p}, 0)$ is locally asymptotically stable for $\lambda_{21} < 0$ which gives the condition

If the condition of (3.12) is violated, $P_2(\overline{n}, \overline{p}, 0)$ transitions into a saddle point. This indicates the instability along the $z$-axis.

3.3.4.   Stability of $P_3(n^*, p^*, E^*)$

The Jacobian matrix computed at the interior point $P_3(n^*, p^*, E^*)$ is given by

The following characteristic equation is determined from the Jacobian matrix above about the point $(n^*, p^*, E^*)$:

with

Moreover,

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 $q_1 = q_2$. In this case, the point $P_3$ 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.

3.4.   Bifurcation analysis

3.4.1.   Transcritical analysis

Theorem 3.1. The system (3.2) undergoes a transcritical bifurcation around the planar equilibrium point $P_2(\overline{n}, \overline{p}, 0)$ as the bifurcation parameter $A$ crosses a critical value such that

Example 3.1. For a particular dataset $k = 0.5, r_1 = 3, r_2 = 1, a = 2, \alpha = 2.5, q_1 = 1, q_2 = 1.5, \mu = 2, b = 4, \beta = 3, c = 4, \text{and}$ $p_0 = 10$, the Jacobian matrix $J_2$ of $P_2(0.293797, 0.521937, 0)$ is confirmed to have a zero eigenvalue at the threshold value of $A_{tc_1} = 0.373911$. The Jacobian matrix $J_2^{*}$ is given below.

The eigenvectors corresponding to the zero eigenvalue (i.e., $d_{33} = 0$) are $V_{e_{1}} = (-0.2557, -0.5283, 0.7348)^T$ and $W_{e_{1}} = (0, 0, 1)^T$. This analysis confirms that all the conditions of Sotomayar's theorem for a transcritical bifurcation at the point $P_2$ are verified.

Theorem 3.2. The system (3.2) undergoes a transcritical bifurcation around the planar equilibrium point $P_1(\hat{n}, 0, \hat{E})$ as the bifurcation parameter $A$ crosses a critical value such that

Example 3.2. For a particular dataset $k = 0.5, r_1 = 3, r_2 = 1, a = 2, \alpha = 2.5, q_1 = 1, q_2 = 1.5, \mu = 2, b = 4, \beta = 3, c = 4,$ and $p_0 = 10$, the Jacobian matrix $J_1$ of $P_1(0.4, 0, 0.3)$ is confirmed to have a zero eigenvalue at the threshold value of $A_{tc_2} = 0.714286$. The Jacobian matrix $J_1^{*}$ is given below.

The eigenvectors corresponding to the zero eigenvalue (i.e., $c_{22} = 0$) are $V_{e_2} = (-0.1641, 0.1641, -0.9727)^T$ and $W_{e_2} = (0, 1, 0)^T$. This analysis confirms that all the conditions of Sotomayar's theorem for a transcritical bifurcation at the point $P_1$ are verified.

3.4.2.   Hopf bifurcation analysis

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 $J_3$ according to Liu's criterion ^[35,36,37].

Liu's criterion: The characteristic equation of the Jacobian matrix at an interior equilibrium point is

It should also satisfy the following conditions for some critical value $A_H$.

1) $\xi_1(A_H) > 0$, $\quad \xi_3(A_H) > 0$, $\quad \Phi(A_H) = \xi_1(A_H)*\xi_2(A_H)-\xi_3(A_H) = 0$

2) $\dfrac{d\Phi(A_H)}{dA_H}\neq0$

Example 3.3. The system (3.2) undergoes a Hopf bifurcation around an interior point $P_3(0.301964, 0.370331, 0.227896)$ for the following particular data: $k = 0.5, r_1 = 3, r_2 = 1, a = 2, \alpha = 2.5, q_1 = 1, q_2 = 1.5, \mu = 2, b = 4, \beta = 3, c = 4,$ and $p_0 = 10$. This Hopf point occurs when the parameter $A$ crosses a threshold value $A_H = 0.388133$. It can be obtained by determining the conditions of Liu's criterion as follows:

Since the conditions required by Liu's criterion are satisfied, the Hopf bifurcation occurs around an interior point $P_3$.

3.5.   Numerical simulation of the model (3.2)

Consider the following data for the system (3.2) with linear harvesting:

The critical values of $A$ are computed as follows:

For the aggregated system (3.2), the boundary equilibrium $P_1 = (4, 0, 0.3)$ is obtained. If we choose $A = 0.35$ $(A > A_0)$, another boundary equilibrium point $P_2 = (0.1297, 0.5263, 0)$ exists. Next, for $A = 0.38$ $(A > A_1)$, the interior equilibrium point $P_3 = (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 $P_3$, some bifurcation points of codimension-1 are detected in the Figure 3 with respect to the bifurcation parameter $A$. There is a branch point $BP_1$ (transcritical bifurcation) at $A = A_{tc_1} = 0.373911$ around the equilibrium point $P_2 = (0.293797, 0.521937, 0)$. One Hopf point is obtained at $A = A_H = 0.388133$ in the interior $\mathbb{R}_+^3$ around the equilibrium point $P_3 = (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 $BP_2$ occurs at $A = A_{tc_2} = 0.714286$ around the equilibrium point $P_1\approx (0.4, 0, 3)$. Figure 3 shows that as the value of $A$ increases close to the branch point $BP_2$ around $A = 0.714$, the solutions in the interior $\mathbb{R}_+^3$ 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 $D_L = 2.4305$ and $D_L = 2.0952$, respectively. Figure 6 shows the Lyapunov exponents for $A = 0.47$ and $A = 0.6$, with the computed dimensions $D_L = 2.25545$ and $D_L = 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 $P_2(\overline{n}, \overline{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 $BP_1$, at $A = 0.374$, the interior point $P_3$ is locally asymptotically stable (see Figure 7B). Considering $A > A_H$, 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 $BP_2$ around $A = 0.714$, the chaotic solution in the interior $\mathbb{R}_+^3$ vanishes, and periodic solutions appear in the $nE$ -plane. The analytical result in (3.3.2) regarding the stability of the point $(\hat{n}, 0, \hat{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 $\mu$ in Figure 10. The branch point $BP_1$ at $\mu = 1.983924$ around $P_2 \approx (0.293843, 0.520833, 0)$ is computed. The Hopf points $H_1$ and $H_2$ are detected, where the Hopf point $H1$ at $\mu = 2.019724$ around $P_3 \approx (0.301694, 0.374208, 0.221869)$ and the Hopf point $H_2$ at $\mu = 2.468001$ around $P_1 \approx (0.4, 0, 3.000001)$ are obtained. For these Hopf points, the corresponding first Lyapunov coefficients are computed as $(-2.964738 \times 10^{-2}) < 0$ and $8.841215 \times 10^{-4} > 0$, respectively. Accordingly, $H_1$ is a supercritical and $H_2$ 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., $\hat{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:

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.

4.   Mathematical model 2 with nonlinear harvesting

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 $H_1(n_1, n_2, E) = \dfrac{q_1En_1}{1+m_1n_1+m_2n_2}$ and $H_2(n_1, n_2, E) = \dfrac{q_2En_2}{1+m_1n_1+m_2n_2}$, and the model is described as follows:

The nonlinear terms $\dfrac{p_0 q_1En_1}{1+m_1n_1+m_2n_2}$ and $\dfrac{p_0q_2En_2}{1+m_1n_1+m_2n_2}$ represent the revenue or benefit from harvesting in Layer-1 and Layer-2, respectively, with the denominator capturing the effect of ecological saturation, while $m_1$ and $m_2$ 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.

4.1.   The aggregated model

The aggregated model for the system mentioned above is as follows:

where $\quad m(p) = \dfrac{m_1 k + m_2 \alpha p}{k + \alpha p}$.

4.2.   Existence of equilibrium points of aggregated model (4.2)

The system described by (4.2) possesses four positive steady states, enumerated as follows:

1) The trivial fixed point $\tilde{P_0}(0, 0, 0)$ is exists and the boundary equilibrium point $\tilde{P_2}(\overline{n}, \overline{p}, 0)$ is feasible for the condition (3.3).

2) The predator-free boundary equilibrium point $\tilde{P_1}(\hat{n}, 0, \hat{E}) = \left(\dfrac{c}{p_0q_1-cm_1}, 0, \dfrac{r_1p_0}{p_0q_1-cm_1}\right)$ exists in the positive $nE$-plane for the condition $p_0q_1 > cm_1$.

3) The unique interior equilibrium point $\tilde{P_3}(n^*, p^*, E^*)$ of (4.2) exists, where $n^*, p^*$, and $E^*$ are obtained as follows:

It can be observed that

Accordingly, $n^*$ and $p^*$ are positive for

The value of $E^*$ is positive for

which gives

Since $A_0 < \tilde{A_1}$, the point $\tilde{P_3}(n^*, p^*, E^*)$ is feasible for the following condition:

That is,

4.3.   Stability analysis of the equilibrium points of the model (4.2)

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:

For the value of $f'(p), r'(p)$ and $q'(p)$, refer to (3.8). Moreover,

4.3.1.   Stability of $\tilde{P_0}(0, 0, 0)$

This is as discussed in Subsubsection (3.3.1).

4.3.2.   Stability of $\tilde{P_1}(\hat{n}, 0, \hat{E})$

At the point $\tilde{P_1}(\hat{n}, 0, \hat{E})$, the Jacobian matrix assumes the following configuration:

The characteristic equation of the matrix $\tilde{J_1}$ about the point $\tilde{P_1}$ is given by

One eigenvalue is $\lambda_{11} = \dfrac{-c}{p_0} < 0$ and the remaining two eigenvalues can be obtained from the characteristic equation

The stability of the point $\tilde{P_1}$ can be obtained by evaluating $trace(\tilde{J_1}) = a_{32} < 0$ and $det(\tilde{J_1}) = a_{31}a_{22} > 0$ from Eq (4.6) for the following condition:

Therefore, the point $\tilde{P_1}$ is locally asymptotically stable. However, if the condition (4.7) is violated, then the point $\tilde{P_1}$ becomes a saddle point. Therefore, bifurcation may occur around $\tilde{P_1}$ provided that

4.3.3.   Stability of $\tilde{P_2}(\overline{n}, \overline{p}, 0)$

The Jacobian matrix computed at the point $\tilde{P_2}(\overline{n}, \overline{p}, 0)$ is given by

The characteristic equation associated with the matrix $\tilde{J_2}(\overline{n}, \overline{p}, 0)$ is given by

One eigenvalue is $\lambda_{21} = -c+\dfrac{p_0q(\overline{p})\overline{n}}{1+m(\overline{p})\overline{n}}$ and the remaining two eigenvalues can be derived from the following characteristic equation:

The $trace(\tilde{J_2})$ and $det(\tilde{J_2})$ of Eq (4.9) are same as those mentioned in Subsubsection (3.3.3). Due to existence condition for the point $\tilde{P_2}$, it can be observed that $trace(\tilde{J_2}) < 0$ and $det(\tilde{J_2}) > 0$. Accordingly we have the following:

1) The point $\tilde{P_2}(\overline{n}, \overline{p}, 0)$ is locally asymptotically stable for $\lambda_{21} < 0$, which gives the condition

2) Moreover, the point $\tilde{P_2}(\overline{n}, \overline{p}, 0)$ transforms into a saddle point when (4.10) is violated. Therefore, bifurcation may occur around $\tilde{P_2}(\overline{n}, \overline{p}, 0)$ for

4.3.4.   Stability of $\tilde{P_3}(n^*, p^*, E^*)$

The Jacobian matrix for the interior point $\tilde{P_3}(n^*, p^*, E^*)$ is as follows:

The characteristic equation for the given Jacobian matrix centered at the point $(n^*, p^*, E^*)$ can be expressed as

with

By applying the Routh–Hurwitz criterion, the interior equilibrium point $\tilde{P_3}(n^*, p^*, E^*)$ is locally asymptotically stable if the following conditions are met:

Example 4.1. In general, to determine the stability of the point $\tilde{P_3}$, it is difficult to analyze the signs $B_1, B_2$, and $B_3$ to verify all the conditions above. Therefore, the stability analysis of the point $\tilde{P_3}(n^*, p^*, E^*)$ is conducted for the particular choice of data, namely $k = 0.8, r_1 = 2, r_2 = 1, a = 2, \alpha = 3.5, q_1 = 1, q_2 = 1.5, \mu = 4, b = 4, \beta = 5, c = 0.03, p_0 = 0.1, m_1 = 0.6, m_2 = 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 $B_{11} = 0.0312 > 0$, $B_{13} = 0.00028961 > 0$, and $B_{11}B_{12} = 0.0312*0.0789 = 0.0025 > B_{13}$. Consequently, the interior point $\tilde{P_3}$ is locally asymptotically stable.

4.4.   Bifurcation analysis

4.4.1.   Transcritical analysis

Theorem 4.1. The system (4.2) undergoes a transcritical bifurcation around the planar equilibrium point $\tilde{P_1}(\hat{n}, 0, \hat{E})$ as the bifurcation parameter $A$ crosses a critical value such that

Proof. The Jacobian matrix $\tilde{J_1}$ of the system (4.2) at the equilibrium point $\tilde{P_1}$ has a zero eigenvalue for the condition (4.13), then the matrix $\tilde{J_1}$ becomes $\tilde{J}_1^*$ (which indicates a matrix when the condition (4.13) is satisfied, i.e., $a_{22} = 0$). Therefore, the Jacobian matrix $\tilde{J}_1^*$ and the transpose of the Jacobian matrix $\tilde{J}_1^{*T}$ have the following eigenvectors corresponding to a zero eigenvalue:

where

The computation for the conditions of Sotomayar's theorem is stated as follows:

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 $\tilde{P_1}$. □

Example 4.2. For a particular dataset $k = 0.8, r_1 = 2, r_2 = 1, a = 2, \alpha = 3.5, q_1 = 1, q_2 = 1.5, \mu = 4, b = 4, \beta = 5, c = 0.03, p_0 = 0.1, m_1 = 0.6,$ and $m_2 = 0.8$, the Jacobian matrix $\tilde{J}_1$ of $\tilde{P_1}(0.365854, 0, 2.439024)$ is confirmed to have a zero eigenvalue at $\tilde{A}_{tc_1} = 0.282759$. The Jacobian matrix $\tilde{J}_1^*$ is given below:

The eigenvectors corresponding to the zero eigenvalue (i.e., $a_{22} = 0$) are $\tilde{V}_{e_1} = (-0.2946, 0.0281, -0.9552)^T$ and $\tilde{W}_{e_1} = (0, 1, 0)^T$. This analysis confirms that all the conditions of Sotomayar's theorem for a transcritical bifurcation at the point $\tilde{P_1}$ are verified such that $\Delta_{11} = 0$, $\Delta_{12}\neq0$, and $\Delta_{13}\neq0$.

Theorem 4.2. The system (4.2) admits a transcrtical bifurcation around the point $\tilde{P_2}(\overline{n}, \overline{p}, 0)$ when the bifurcation parameter $A$ varies such that

Proof. The Jacobian matrix $\tilde{J_2}(\overline{n}, \overline{p}, 0)$ of the system (4.2) at the equilibrium point $\tilde{P_2}$ has a zero eigenvalue for the condition (4.15), then the matrix $\tilde{J_2}$ becomes $\tilde{J}_2^*$ (which indicates a matrix when the condition (4.15) is satisfied, i.e., $b_{33} = 0$). Therefore, the Jacobian matrix $\tilde{J}_2^*$ and the transpose of the Jacobian matrix $\tilde{J}_2^{*T}$ have the following eigenvectors corresponding to the zero eigenvalue:

where,

The computation for the conditions of Sotomayar's theorem is stated as follows:

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 $\tilde{P_2}(\overline{n}, \overline{p}, 0)$. □

Example 4.3. For a particular dataset $k = 0.8, r_1 = 2, r_2 = 1, a = 2, \alpha = 3.5, q_1 = 1, q_2 = 1.5, \mu = 4, b = 4, \beta = 5, c = 0.03, p_0 = 0.1, m_1 = 0.6,$ and $m_2 = 0.8$, the Jacobian matrix $\tilde{J_2}$ of $\tilde{P_2}(0.271750, 0.415937, 0)$ is confirmed to have a zero eigenvalue at $\tilde{A}_{tc_2} = 0.216708$. The Jacobian matrix $\tilde{J}_2^*$ is given below:

The eigenvectors corresponding to the zero eigenvalue (i.e., $b_{33} = 0$) are $\tilde{V}_{e_2} = (-0.2233, -0.5297, 0.8182)^T$ and $\tilde{W}_{e_2} = (0, 0, 1)^T$. This analysis confirms that all the conditions of Sotomayar's theorem for a transcritical bifurcation at the point $\tilde{P_2}$ are verified such that $\Delta_{21} = 0$, $\Delta_{22}\neq0$, and $\Delta_{23}\neq0$.

4.4.2.   Hopf bifurcation analysis

Example 4.4. For the parameter values $k = 0.8, r_1 = 2, r_2 = 1, a = 2, \alpha = 3.5, q_1 = 1, q_2 = 1.5, \mu = 4, b = 4, \beta = 5, c = 0.03, p_0 = 0.1, m_1 = 0.6,$ and $m_2 = 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 = \tilde{A}_H = 0.223632$ are met. The conditions for the Liu's criterion at the bifurcation parameter $\tilde{A}_H$ are as follows:

4.5.   Numerical solutions for model 2

Consider the following choice of hypothetical data for the system (4.2) with nonlinear harvesting in appropriate units:

The critical value of A are computed as follows:

The aggregated system reaches its equilibrium boundary at $\tilde{P_1} \approx (0.3659, 0, 2.4390)$. $A$ is chosen as $A = 0.21$ ($A_0 < A < \tilde{A}_1$), then another equilibrium point on the boundary, denoted as $\tilde{P}_2(0.1681, 0.4171, 0)$, is identified and presented in Figure 19. Following that, with $A = 0.22$, an interior equilibrium point is obtained $\tilde{P}_3(0.2778, 0.3301, 0.1555)$ as shown in Figure 20. (see the example in (4)).

Following the interior point $\tilde{P}_3$, 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 $BP_1$ (also known as transcritical bifurcation), occurs at $A = \tilde{A}_{tc_1} = 0.216708$ around $\tilde{P_2} = (0.271750, 0.415937, 0)$.

2) Additionally, a Hopf point emerges at $A = \tilde{A}_H = 0.223632$ within the interior $\mathbb{R}_{+}^3$ around the point $\tilde{P_3} = (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 $BP_2$ is obtained at $A = \tilde{A}_{tc_2} = 0.282759$ around the point $\tilde{P_1}\approx(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 $\tilde{P_2}(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 $\tilde{P_3}(0.2778, 0.3301, 0.1555)$. When $A > \tilde{A}_H = 0.223632$, the system becomes destabilized, resulting in the appearance of strange attractors, as illustrated in Figures 21–23.

5.   Conclusions and discussion

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 ($P_0$), the predator-free boundary equilibrium ($P_1$), the boundary equilibrium in the absence of harvesting ($P_2$), and the interior equilibrium point ($P_3$). 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 $P_2$ and $P_3$ when $A < A_{tc_1}$ and $A_{tc_1} < A < A_H$, respectively. For values of $A$ within the range $A_H < A < A_{tc_2}$, the system exhibits chaotic behavior. Beyond $A_{tc_2}$, 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 $\tilde{P_3}$ (the interior equilibrium point), as well as two transcritical bifurcations around the equilibrium points $\tilde{P_1}$ (the predator-free equilibrium point) and $\tilde{P_2}$ (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 $\tilde{P_2}$ point when $A_0 < A < \tilde{A}_1$, which implies no harvesting of prey species is possible due to an insufficient amount of prey, and the system stabilizes at $\tilde{P_3}$ for the range $\tilde{A}_1 < A < \tilde{A}_2$. The coexistence of prey and predator populations occurs within this suitable range of additional food ($\tilde{A}_1 < A < \tilde{A}_2$). Moreover, the study demonstrated that the system exhibits local asymptotic stability at $\tilde{P_2}$ when $A < \tilde{A}_{tc_1}$ and at $\tilde{P_3}$ when $\tilde{A}_{tc_1} < A < \tilde{A}_H$.

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.

Use of AI tools declaration

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

Acknowledgments

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.

Conflict of interest

The authors declare there is no conflict of interest.