1. Introduction
The exploitation of natural resources has increased as a result of the rising need for food and energy. Therefore, in theoretical ecology and evolutionary biology throughout the past few decades, many scientists have examined prey–predator interactions, and mathematical models have tremendously benefited in better understanding these challenging scenarios [
1,
2,
3,
4]. In the early Nineteenth Century, Malthus seems to have utilized mathematical models to explain the patterns of prey–predator relationships. In order to simulate prey–predator interactions with realistic accuracy, the well-known Lotka–Volterra model was modified to include a prey logistic growth factor, as well as a number of population-dependent reaction functions [
1,
3].
In order to improve the birth rate of prey animals, fear effects were also exploited. Numerous studies in basic ecology and environmental biology have examined the impact of fear [
5,
6,
7,
8,
9,
10,
11,
12,
13,
14,
15]. In addition to preventing direct ingestion, fear of predators has been found to increase vigilance and reduce the amount of time spent searching for wild animals living in social groups when population sizes decline. When averaged across numerous trials, the impact of population-level anxiety on prey survival may be comparable to that of direct predator eating. The cost of fear in prey reproduction, which was recently developed and tested by Sarkar and Khajanchi [
10], is utilized in this paper. They demonstrated how strong anti-predator reflexes could maintain relationships between prey and predator. Regarding fear, Maghool and Naji [
12] built and researched a tri-trophic Leslie–Gower food-web system. A time-delayed predator–prey model with Holling type II functional response was put forth by Liu et al. [
13] and includes the gestation period and the cost of fear on prey reproduction. In their study of the dynamics of a Leslie–Gower ratio-dependent predator–prey model that took into consideration the prey population’s fear, Vinoth et al. [
14] also took into account the Allee effect in the predator growth from both a biological and mathematical standpoint. Due to the prey group’s protective skills, the Sokol–Howell kind of function response is employed to describe the predation process. The effects of predation anxiety on the dynamics of the three-species food chain system at the first and second levels were postulated and studied by Maghool and Naji [
15]. An investigation of the dynamics of a prey–predator model took into account the predator stage structure, and the fear effect was performed recently in [
16].
A refuge is a term used in ecology to describe a location where an organism might hide from predators in order to avoid being discovered. Due to population dynamics, populations of both predators and prey are much larger, and a region can support a significant number of additional species when shelters are present. Therefore, any method that lowers the risk of predation, such as prey aggregations, geographic or temporal refuges, or reduced prey search activity, can be broadly referred to as a refuge. There is a natural occurrence that serves as a haven for it to survive. The prey’s use of regional shelters is one of the more significant behavioral traits that influences the dynamics of predator–prey systems. By using shelters, some of the prey population is partially protected from predators. The coexistence of the predator–prey system is significantly impacted by the presence of shelters. The most widely used definitions of refuge in the literature at this time are a constant refuge and refuge proportionate to prey; see for example [
17,
18] and the references therein. Researchers are now being drawn to novel, diverse refuge concepts [
19,
20,
21]. Due to the danger that predators cause, this study explores a refuge notion proportional to the predator. Prey refuges grow as predators do because as fears grow, so do the predators’ numbers.
On the other hand, harvesting is a substantial and regular occurrence. Fisheries frequently use harvesting since natural systems are largely regenerative. In an exploited fishing system with two interacting species, researchers are looking into capturing either prey or predator species or both prey and predator species. Many different harvesting techniques have been used. Some people use continuous threshold harvesting, fixed effort harvesting (or proportional harvesting), and constant harvesting [
20,
22,
23,
24], while others researched nonlinear harvesting [
18,
25,
26].
The Leslie–Gower model of prey and predators places a strong emphasis on the fact that both prey and predator growth rates have upper bounds and that the quantity of available prey affects the predator’s carrying capacity. This is not taken into consideration by the Lotka–Volterra model. Both predators and prey can reach these upper limits under ideal circumstances: for predators, when there is much prey per predator, and for prey, when there are few predators. Numerous scholars have become interested in Leslie–Gower prey–predator systems [
4,
12,
14,
25,
26,
27]. The Leslie–Gower model of prey and predator in the presence of fear, quadratic fixed effort harvesting, and predator-dependent refuge is therefore suggested and explored in this study.
2. Model Construction
The modified Leslie–Gower prey–predator model with the Sarkar and Khajanchi fear function [
10] that influences the prey’s birth rate and the quadratic fixed effort harvesting with the Beddington–DeAngelis type of functional response can be stated in the following form:
where
and
are the prey and predator densities at time
T, respectively.
Due to the importance of the prey’s refuge, it is observed that prey refugees are thought to reduce predator–prey oscillations and avoid prey extinction. An examination of the empirical evidence reveals that refuges can fulfill the former function. Therefore, the Beddington–DeAngelis response function is used to create a harvested in the above dynamical model so that the quantity of prey refuge is dependent on both species. Assume that the amount of prey refuge is
[
19], with
being the refuge coefficient. The remaining
prey species are vulnerable to predators. Throughout the manuscript, it is also assumed that
, which is due to the fact that only a portion of the prey will have a refuge, and
ensures an appropriate range of refuge for a realistic environment.
Accordingly, the dynamics of the above-described model can be represented in the following form:
where all the parameters are nonnegative and described in
Table 1.
Theorem 1. System (2) is a positively invariant.
Proof. According to the form of System (2), it is clear that the system is a Kolmogorov system in which
and
are continuously differentiable functions representing the growth rates of the prey and predator, respectively. Hence, using the positive conditions
, we can solve (2) to obtain:
Consequently, from the definition of the exponential function, any solution in the interior of that starts with positive initial conditions stays there indefinitely, according to the previous two equations. □
Theorem 2. All the solutions of System (2) are uniformly bounded in the region.where
, and are positive constants that satisfy
, which represents the survival condition of the prey species in the absence of the predator, while all the new symbols are given in the proof. Proof. According to the first equation of System (2), it is observed that
As a result of solving the differential inequality above:
Therefore, when
, it is obtained that
, since the birth rate for the survivor species is biologically greater than the death rate. Using the prey’s bound in the second equation of System (2), it is deduced that:
Similarly, solving the last differential inequality gives:
Therefore, as , we obtain . Consequently, the proof is complete. □
In light of the preceding, System (2) comprises continuous partial derivative interaction functions in the domain. As a result, for any given beginning condition, System (2) has a unique solution.
3. Equilibria and Their Stability
There are only four non-negative equilibrium points (EP) in System (2), and the existing conditions and their forms are as follows:
The trivial equilibrium point (TEP), which is defined by , always exists.
The predator-free equilibrium point (PDFEP) is denoted by
, where
, which can be obtained from solving
, exists on the
axis if and only if
The prey-free equilibrium point (PYFEP), which is written as , where , which can be obtained from solving , always exists.
Note that all the above EPs are the same as those of System (1).
Finally, there is the survival equilibrium point (SEP), which is expressed as
, where
However,
is the positive root of the next sixth-order polynomial equation.
where:
Therefore, the SEP exists in the interior of a positive quadrant of the
plane uniquely if the following conditions are met.
with one set of the following sets of sufficient conditions.
In order to study the local asymptotic stability (LAS) of System (2), the Jacobian matrix (JM) of System (2) is computed at the point
by:
where:
with
,
.
Accordingly, for
, the JM becomes
Hence, has the eigenvalues and . Thus, is an unstable node or saddle point depending on , or , respectively.
The JM at
can be written as:
Thus, has the following eigenvalues due to the existing condition (3), and . Thus, is a saddle point.
The JM at
can be written as:
Therefore, the following roots
and
are the eigenvalues of
. Therefore, the PYFEP is LAS provided the following requirement is met.
Finally, System (2) JM at
can be written as:
where:
where
and
. Direct computation reveals that
has two negative real-part eigenvalues, implying that
is LAS if the following necessary condition is satisfied.
A fundamental challenge in mathematical biology is the long-term persistence of each component of a system of interacting components, which is often a population in an ecological environment. Long-term survival is determined by a variety of factors. As a result, studying each species’ long-term survival in a prey–predator relationship is a fascinating topic. In this section, the Gard technique is utilized, which is based on constructing Lyapunov-like persistence functions [
28]. The criteria for the system’s survival are specified in the following theorem.
Theorem 3. If the following condition holds, then System (2) is uniformly persistent.
Proof. Consider the function , where and are positive constants.
Obviously,
for all
, and
for all
, where
represents the boundary of
. Hence,
is the so-called persistence function or named the average Lyapunov function in the sense of the Gard approach. Therefore, according to Gard, the proof follows if and only if
is positive for all points
that belong to the
limit sets of System (2) in
. The following is the result of direct computation:
Accordingly,
,
, and
belong to the
limit sets of System (2) in
. Moreover, the direct calculation gives that:
Clearly, for a suitable choice of and , so that is sufficiently larger than , it is obtained that , while always. However, under Condition (14). Hence, the proof is complete. □
4. Global Stability Analysis
The following theorems are concerned with the system’s global dynamics (2). The attractive basin of trajectories of a dynamical System (2), according to global stability, is either the state-space or the interior of the state-space that determines the system’s state variables. In other words, global stability implies that, regardless of the initial conditions, all paths eventually drift to the system’s attractor. Global stability is required by most biological systems, such as gene regulatory systems.
Theorem 4. The PYFEP is globally asymptotically stable (GAS) if and only if the following requirement holds.
Proof. Define the function , which is a positive definite real-valued function on the region . Then, after performing some direct calculations, it is deduced that
Further simplification leads to:
Obviously, is negative definite under Conditions (15) and (16). Hence, the PYFEP is GAS. □
Theorem 5. Assume that System (2) has a unique SEP that is LAS, then it is GAS.
Proof. If there exists a continuously differentiable function
(called the Dulac function) such that the expression
has the same sign
almost everywhere in a simply connected region of the plane, the plane autonomous System (2) has no non-constant periodic solutions lying entirely within the region, according to the Bendixson–Dulac theorem [
29].
Hence, define
, which is a
in a simply connected region of the domain of System (2), then direct computation with respect to System (2) shows that
Clearly, due to Condition (12) of the local stability of the SEP.
Thus, System (2) has no periodic dynamics in the interior of , which is assumed to have a unique SEP. Consequently, the Poincare–Bendixson theorem that states the solution guarantees that the SEP is GAS. □
5. Bifurcation Analysis
Bifurcation theory investigates changes in the qualitative structure of a collection of curves, such as the integral curves of a set of vector fields or the solutions of a set of differential equations. A bifurcation occurs when a small smooth change in a system’s parameter values causes a significant qualitative change in its behavior. It is most commonly used in the mathematical study of dynamical systems. Bifurcation can take two different forms: local bifurcations, which can be seen when parameters cross critical thresholds by looking for changes in the local stability properties of equilibria, periodic orbits or other invariant sets; global bifurcations, which happen when the system’s larger invariant sets interfere with each other or with the system’s equilibria. They cannot be found solely by examining the stability of the equilibria. In this section, the detection of the possibility of having local bifurcation is carried out.
Hence, the second directional derivative of
, where
is any vector with
, can be written using direct computation as:
where:
Theorem 6. System (2) undergoes a transcritical bifurcation (TB) at the PYFEP when the parameter passes through the value, provided that the following condition holds: where the symbol
is determined in the proof. Proof. From the JM that is given by Equation (9), it is observed that for
, it becomes
Therefore,
has eigenvalues given by
and
. Thus,
becomes a non-hyperbolic point. Let
and
be the eigenvectors corresponding to the
of
and their transpose, respectively. Then, the direct calculation gives that:
Moreover, simple computation gives that:
. Hence, it is obtained that
.
Moreover, according to Equation (17), the following result is obtained
where:
Hence, it is obtained that
, which is not equal to zero under Condition (18). As a result of the Sotomayor theorem of local bifurcation [
29], System (2) undergoes a TB at
, and this completes the proof. □
Theorem 7. If Condition (12) is met along with the following condition:
then, as the parameter, System (2) presents a saddle-node bifurcation (SNB) near the SEP provided that the following condition holds. where: All the new symbols are defined in the proof.
Proof. Recall the JM of System (2) at the SEP that is given by Equation (11) with
, then it can be written as:
where
, and
are the elements of
, with
. Direct calculation shows that the determinant of
is zero.
Then has a zero eigenvalue () with the second eigenvalue under Condition (12). Thus, the SEP is a non-hyperbolic point when .
Let
and
be the eigenvectors associated with
of
and their transpose, respectively. Then, the direct calculation gives that:
Note that, according to the Jacobian elements, we obtain that
and
due to Conditions (12) and (19). Moreover, simple computation gives that:
Hence, It is concluded that
. Moreover, according to Equation (17), the following is gained:
Therefore, the expression , is not identical to zero due to Condition (20). Thus, System (2) presents an SNB near the SEP, which completes the proof. □
Theorem 8. Assume that Condition (13) is satisfied along with the following conditions:where
, andare the JM elements that are given in Equation (11). Then, System (2) presents a Hopf bifurcation around the SEP when the parameter passes through the value with Proof. Direct computation shows that the characteristic equation of the JM of System (2) at SEP can be determined as:
Note that at
, it is easy to verify that
, due to Condition (21). Consequently, in this case, the roots of Equation (23) are written as:
where
and
. Clearly, the above eigenvalues are complex conjugate pure imaginary under Condition (22). However, in the neighborhood of
, Equation (23) has two complex conjugate eigenvalues given by:
Now, since
therefore, according to the Hopf bifurcation theorem, System (2) has a Hopf bifurcation at
, and this completes the proof. □
6. Numerical Simulation
In this section, to illustrate several dynamical scenarios in accordance with the theoretical findings in the preceding sections, numerical simulations are performed. The following hypothetical set of parameter values, which are biologically plausible, are utilized to numerically solve System (2).
The trajectories of System (2) are shown in the figure after being acquired using Dataset (24) starting from various initial values (1).
According to
Figure 1, System (2) has a unique SEP in its state-space
, which is GAS. Now, a numerical investigation of the effects of changing the parameters on the dynamics of System (2) is conducted. It is noted that System (2) has multiple SEPs in its state-space
with various stability types for the range
; for instance, look at
Figure 2 below at
. System (2) does not have a SEP; however, for the range
, the solution approaches the PYFEP asymptotically, as shown in
Figure 3. Otherwise, System (2) has a singular SEP that is GAS in its state-space.
According to
Figure 2, System (2) has a bi-stable case in which System (2) approaches asymptotically two different SEPs depending on the initial points using the same dataset.
The influence of varying the parameter
on the dynamics of System (2) is presented in
Figure 4, for the values
, keeping the rest of the parameters as in Dataset (24).
As shown in
Figure 4, the population of
decreases while that of
increases with the increase of
. The influence of the fear parameters on the dynamics of System (2) is studied in
Figure 5 and
Figure 6. It is observed that, for the range
, System (2) has no SEP and the system approaches the PYFEP; see
Figure 5; however, varying the parameter
, such that
with the other parameters as in (24), is investigated in
Figure 6.
As shown in
Figure 6, the populations of
and
increase with the increase of
. Now, varying the parameter
, such that
, with the other parameters as in (24) is investigated in
Figure 7. As d increases, the SEP converges to the PYFEP and, then, coincides with it. However, for the range
. System (2) has multiple SEPs; otherwise, the SEP is still unique in the interior of state-space
; see
Figure 8, for instance.
From
Figure 7, it is deduced that increasing
causes the decrease in both species
and
, with
decreasing more and faster. Moreover, it is observed that the increase in the value of
so that
has a similar influence as that with
on the dynamics of System (2). On the other hand, System (2) has a unique SEP that is GAS in its state-space for
.
Now that the parameter
has been varied in the ranges of
,
, and
, it is clear from
Figure 9 that System (2) approaches a PYFEP globally, two separate EPs locally, and a unique SEP globally, respectively.
Now that the parameter
has been varied in the ranges of
,
, and
, it is clear from
Figure 10 that System (2) approaches a PYFEP globally, two separate EPs locally, and a unique SEP globally, respectively.
Obviously, as shown in
Figure 10, increasing
leads to an increase in the population of
and a decrease in
. The influences of changing the harvesting parameters
, and
are investigated and presented in
Figure 11.
It is clear from
Figure 11 that
and
have the opposite influence on the populations of System (2), while both populations decrease as
increases.
For different values of
with the rest of the parameters as in Dataset (24), the dynamics of System (2) is drawn in
Figure 12.
Finally, the influence of the predator-dependent refuge parameter
is shown in
Figure 13.
7. Discussion
This study modified the Leslie–Gower fear and harvesting prey–predator model to include a predator-dependent refuge function. The goal was to comprehend how this kind of refuge affected the dynamics of the prey–predator system. The system underwent theoretical and numerical analysis. The system was shown to have a maximum of three boundary equilibrium points with no, a single, or many SEPs belonging to the state-space defined by . According to the theoretical results, two of the border EPs are unstable, while the third one can either be an unstable saddle or a stable node (locally or globally) depending on whether the SEP is present. The uniquely existing SEP, on the other hand, is GAS or there is a periodic dynamic surrounding it (Hopf bifurcation). However, the numerical simulation that was achieved for the chosen hypothetical set of data revealed a rich dynamical behavior that may be summed up in the next section.
The existence of periodic dynamics due to the Hopf bifurcation was proven theoretically; however, the set of selected hypothetical data (24) does not support their existence, and it is still possible to have different sets of data. The system has a GAS point attractor that may transfer to being LAS for multiple points with a disjoint basin of attractions in the case of having more than one SEP in the interior of System (2)’s state-space.
8. Conclusions
It is concluded that the birth rate of the prey population has two bifurcation points at which the stability of System (2) is transferred from the PYFEP to a unique SEP and, then, to two different SEPs. However, the birth rate of the predator population has no bifurcation point; instead, it causes a quantitative influence on the size of the species’ populations. Regarding the fear function parameters, decreasing the minimum cost of fear leads to extinction of the prey species and the system approaches asymptotically the PYFEP globally. However, increasing the level of fear keeps the persistence of the system and increases the population’s size. It was obtained also that an increase in the natural prey’s death rate causes extinction of the prey, and the system is GAS at the PYFEP. For small values of the prey’s intraspecific competition, the system undergoes a bi-stable behavior due to the existence of two SEPs, which are LAS. However, increasing the prey’s intraspecific competition above a specific value gives a unique SEP in the state-space, which is GAS. The attack rate has a similar influence on the dynamics of System (2) as that of the prey’s natural death rate. For the lower value of a half-saturation constant, there is no SEP and System (2) approaches asymptotically the PYFEP. Once the half-saturation constant crosses that lower value, System (2) has bi-stable behavior between the PYFEP and SEP with a disjoint basin of attraction for them. Further, with increases in this parameter, the basin of attraction of the SEP becomes the whole of the state-space, which makes the SEP GAS. A similar behavior is shown as that of the half-saturation constant when varying the food quantity conversion rate to a predator. The influence of varying harvesting parameters on the dynamics of System (2) is quantitative, so that the catchability coefficients of the prey and predator have the opposite impact on each other on the population size. However, increasing the effort level leads to decreasing in both populations. On the other hand, increasing the carrying capacity of the predator leads to an increase in the predator population and a decay in the prey population. Moreover, varying the refuge parameter has a great impact on the dynamics of System (2), so that as the parameter increases, the system transfers from having a GAS point at a unique SEP to the bi-stable behavior between two SEPs in the state-space, then returns back again to the GAS at a unique SEP.
Finally, for future work, it is advised that the prey–predator model be expanded to become a food-web by including a top predator or else embedding one or more biological factors, such as stage structure, infectious diseases, delay, diffusion, and time-dependent parameters.