Next Article in Journal
Solution Bounds and Numerical Methods of the Unified Algebraic Lyapunov Equation
Previous Article in Journal
The Stochastic Team Orienteering Problem with Position-Dependent Rewards
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modeling and Analysis of the Influence of Fear on the Harvested Modified Leslie–Gower Model Involving Nonlinear Prey Refuge

by
Abdul Rahman Mahmoud Jamil
and
Raid Kamel Naji
*
Department of Mathematics, College of Science, University of Baghdad, Baghdad 10071, Iraq
*
Author to whom correspondence should be addressed.
Mathematics 2022, 10(16), 2857; https://doi.org/10.3390/math10162857
Submission received: 19 July 2022 / Revised: 3 August 2022 / Accepted: 4 August 2022 / Published: 10 August 2022
(This article belongs to the Section Mathematical Biology)

Abstract

:
Understanding the effects of fear, quadratic fixed effort harvesting, and predator-dependent refuge are essential topics in ecology. Accordingly, a modified Leslie–Gower prey–predator model incorporating these biological factors is mathematically modeled using the Beddington–DeAngelis type of functional response to describe the predation processes. The model’s qualitative features are investigated, including local equilibria stability, permanence, and global stability. Bifurcation analysis is carried out on the temporal model to identify local bifurcations such as transcritical, saddle-node, and Hopf bifurcation. A comprehensive numerical inquiry is carried out using MATLAB to verify the obtained theoretical findings and understand the effects of varying the system’s parameters on their dynamical behavior. It is observed that the existence of these factors makes the system’s dynamic behavior richer, so that it involves bi-stable behavior.

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:
d N d T = N [ r 1 ( m + n ( 1 m ) n + P ) d b N a P c + N + e P q 1 E N ] = N f ( N , P ) , d P d T = P [ r 2 ( 1 e P K + N ) q 2 E P ] = P g ( N , P ) ,                                                                                           .
where N ( T ) and P ( T ) 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 δ N P [19], with δ being the refuge coefficient. The remaining ( N δ N P ) prey species are vulnerable to predators. Throughout the manuscript, it is also assumed that 0 < δ < 1 , which is due to the fact that only a portion of the prey will have a refuge, and 0 1 δ P 1 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:
d N d T = N [ r 1 ( m + n ( 1 m ) n + P ) d b N a ( 1 δ P ) P c + N ( 1 δ P ) + e P q 1 E N ] = N f ( N , P ) , d P d T = P [ r 2 ( 1 e P K + N ( 1 δ P ) ) q 2 E P ] = P g ( N , P ) ,                                                                                            
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 f ( N , P ) and g ( N , P ) are continuously differentiable functions representing the growth rates of the prey and predator, respectively. Hence, using the positive conditions ( N ( 0 ) , P ( 0 ) ) , we can solve (2) to obtain:
N ( T ) = N ( 0 )   e 0 t [ r 1 ( m + n ( 1 m ) n + P ( s ) ) d b N ( s ) a ( 1 δ P ( s ) ) P ( s ) c + N ( s ) ( 1 δ P ( s ) ) + e P ( s ) q 1 E N ( s ) ] d s P ( T ) = P ( 0 ) e 0 t [ r 2 ( 1 e P ( s ) K + N ( s ) ( 1 δ P ( s ) ) ) q 2 E P ( s ) ] d s
Consequently, from the definition of the exponential function, any solution in the interior of + 2 = { ( N , P ) 2 : N ( T ) 0 , P ( T ) 0 } that starts with positive initial conditions ( N ( 0 ) , P ( 0 ) ) stays there indefinitely, according to the previous two equations. □
Theorem 2.
All the solutions of System (2) are uniformly bounded in the region.
Υ = { ( N , P ) + 2 : N [ 0 , r 1 d b ] ,   P [ 0 , r 2 δ ] } ,
where  r 1 ,   r 2 ,   d , and b  are positive constants that satisfy  r 1 d > 0 , 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
d N d T ( r 1 d ) N b N 2 .
As a result of solving the differential inequality above:
N ( T ) r 1 d b [ 1 e ( r 1 d ) T ] + ( r 1 d ) N ( 0 ) e ( r 1 d ) T
Therefore, when T , it is obtained that N ( T ) r 1 d b = ϵ > 0 , 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:
d P d T r 2 P [ r 2 e K + ϵ + q 2 E ] P 2 = r 2 P δ P 2 .
Similarly, solving the last differential inequality gives:
P ( T ) r 2 δ   [ 1 e r 2 T ] + r 2 N ( 0 ) e r 2 T .
Therefore, as T , we obtain P ( T ) r 2 δ . Consequently, the proof is complete. □
In light of the preceding, System (2) comprises continuous partial derivative interaction functions in the + 2 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 W 1 = ( 0 , 0 ) , always exists.
The predator-free equilibrium point (PDFEP) is denoted by W 2 = ( N ¯ , 0 ) , where N ¯ = r 1 d b + q 1 E , which can be obtained from solving f ( N , 0 ) = 0 , exists on the N axis if and only if
d < r 1 .
The prey-free equilibrium point (PYFEP), which is written as W 3 = ( 0 , P ˜ ) , where P ˜ = r 2 K r 2 e + K q 2 E , which can be obtained from solving g ( 0 , P ) = 0 , 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 W 4 = ( N , P ) , where
N = e r 2 P K ( r 2 E q 2 P ) ( 1 δ P ) ( r 2 E q 2 P )
However, P is the positive root of the next sixth-order polynomial equation.
B 6 P 6 + B 5 P 5 + B 4 P 4 + B 3 P 3 + B 2 P 2 + B 1 P + B 0 = 0 ,
where:
B 6 = a δ 2 E 2 q 2 2 < 0 , B 5 = δ E q 2 [ E q 2 ( 2 a a δ n + d e e m r 1 ) + 2 a δ r 2 ] , B 4 = E 2 q 2 2 [ a + c d δ d e d δ K + b e K + 2 a δ n + d δ e n + e E K q 1   c δ m r 1 + e m r 1 + δ K m r 1 δ e n r 1 ] + E q 2 [ 4 a δ r 2 3 d δ e r 2 , + b e 2 r 2 + 2 a δ 2 n r 2 + e 2 E q 1 r 2 + 3 δ e m r 1 r 2 ] a δ 2 r 2 2 , B 3 = E 2 q 2 2 [ c d + b c K + d K b K 2 a n + c d δ n d e n d δ K n + b e K n   + c E K q 1 E K 2 q 1 + e E K n q 1 + c m r 1 K m r 1 c δ n r 1 + e n r 1 + δ K n r 1 ] + E q 2 [ 2 a r 2 2 c d δ r 2 + b c e r 2 + 3 d e r 2 + 2 d δ K r 2 4 b e K r 2 4 a δ n r 2 3 d δ e n r 2 + b e 2 n r 2 + c e E q 1 r 2 4 e E K q 1 r 2 + e 2 E n q 1 r 2 + 2 c δ m r 1 r 2 3 e m r 1 r 2 2 δ K m r 1 r 2 + 3 δ e n r 1 r 2 ] + 2 a δ r 2 2 + 2 d δ e r 2 2 2 b e 2 r 2 2 a δ 2 n r 2 2 2 e 2 E q 1 r 2 2 2 δ e m r 1 r 2 2 . B 2 = E 2 q 2 2 [ c d n + b c K n + d K n b K 2 n + c E K n q 1 E K 2 n q 1 + c n r 1 K n r 1 ] + E q 2 r 2 [ 2 c d 2 b c K 2 d K + 2 b K 2 + 2 a n 2 c d δ n + b c e n + 3 d e n + 2 d δ K n 4 b e K n ] + E 2 q 1 q 2 r 2 [ 2 c K + 2 K 2 + c e n 4 e K n ] + E q 2 r 1 r 2 [ 2 c m + 2 K m + 2 c δ n 3 e n 2 δ K n ] + r 2 2 [ a + c d δ b c e 2 d e d δ K + 3 b e K + 2 a δ n + 2 d δ e n 2 b e 2 n c e E q 1 + 3 e E K q 1 2 e 2 E n q 1 c δ m r 1 + 2 e m r 1 + δ K m r 1 2 δ e n r 1 ] . , B 1 = E n q 2 r 2 [ 2 c ( d b K ) 2 K ( d b K ) 2 E K q 1 ( c K ) 2 r 1 ( c K ) ] + r 2 2 [ c d + b c K + d K r 2 2 b K 2 a n + c d δ n b c e n 2 d e n d δ K n + 3 b e K n + c E K q 1 r 2 2 E K 2 q 1 c e E n q 1 + 3 e E K n q 1 + c m r 1 K m r 1 c δ n r 1 + 2 e n r 1 + δ K n r 1 ] , , B 0 = n r 2 2 ( c K ) [ E K q 1 + r 1 ( d b K ) ] .
Therefore, the SEP exists in the interior of a positive quadrant of the N P plane uniquely if the following conditions are met.
0 < K ( r 2 E q 2 P ) < e r 2 P ,
with one set of the following sets of sufficient conditions.
B 0 > 0 ,   B 1 > 0 ,   B 2 > 0 ,   B 4 < 0 ,   B 5 < 0 B 0 > 0 ,   B 1 > 0 ,   B 2 > 0 ,   B 3 > 0 ,   B 4 > 0 ,   B 5 < 0 B 0 > 0 ,   B 1 > 0 ,   B 2 < 0 ,   B 3 < 0 ,   B 4 < 0 ,   B 5 < 0 B 0 > 0 ,   B 1 > 0 ,   B 2 > 0 ,   B 3 > 0 ,   B 4 > 0 ,   B 5 > 0 B 0 > 0 ,   B 1 < 0 ,   B 2 < 0 ,   B 3 < 0 ,   B 4 < 0 ,   B 5 < 0 } .
In order to study the local asymptotic stability (LAS) of System (2), the Jacobian matrix (JM) of System (2) is computed at the point ( N , P ) by:
J ( N , P ) = ( N f N + f ( N , P ) N f P P g N P g P + g ( N , P ) ) ,
where:
f N = b + a ( 1 δ P ) 2 P Ω 1 2 q 1 E , f P = r 1 n ( 1 m ) ( n + P ) 2 [ c + N ( 1 δ P ) + e P ] ( a 2 δ P ) a ( 1 δ P ) P ( N δ + e ) Ω 1 2 , g N = r 2 e P ( 1 δ P ) Ω 2 2 , g P = q 2 E r 2 e ( K + N ) Ω 2 2 ,
with Ω 1 = c + N ( 1 δ P ) + e P , Ω 2 = K + N ( 1 δ P ) .
Accordingly, for W 1 , the JM becomes
J ( W 1 ) = ( r 1 d 0 0 r 2 ) .
Hence, J ( W 1 ) has the eigenvalues λ 11 = r 1 d and λ 12 = r 2 > 0 . Thus, W 1 is an unstable node or saddle point depending on r 1 > d , or r 1 < d , respectively.
The JM at W 2 can be written as:
J ( W 2 ) = ( ( r 1 d ) N ¯ ( r 1 ( 1 m ) n + a c + N ¯ ) 0 r 2 ) .
Thus, J ( W 2 ) has the following eigenvalues λ 21 = ( r 1 d ) < 0 due to the existing condition (3), and λ 22 = r 2 > 0 . Thus, W 2 is a saddle point.
The JM at W 3 can be written as:
J ( W 3 ) = ( r 1 ( m + n ( 1 m ) n + P ˜ ) d a ( 1 δ P ˜ ) P ˜ c + e P ˜ 0 r 2 e P ˜ 2 ( 1 δ P ˜ ) K 2 r 2 ) .
Therefore, the following roots λ 31 = r 1 ( m + n ( 1 m ) n + P ˜ ) d a ( 1 δ P ˜ ) P ˜ c + e P ˜ and λ 31 = r 2 < 0 are the eigenvalues of J ( W 3 ) . Therefore, the PYFEP is LAS provided the following requirement is met.
r 1 ( m + n ( 1 m ) n + P ˜ ) < d + a ( 1 δ P ˜ ) P ˜ c + e P ˜ .
Finally, System (2) JM at W 4 can be written as:
J ( W 4 ) = ( w 11 w 12 w 21 w 22 )
where:
w 11 = N ( b + a ( 1 δ P ) 2 P Ω 1 2 q 1 E ) , w 12 = N ( r 1 n ( 1 m ) ( n + P ) 2 + [ c + N ( 1 δ P ) + e P ] ( a 2 δ P ) a P ( 1 δ P ) ( e N δ ) Ω 1 2 ) , w 21 = r 2 e P 2 ( 1 δ P ) Ω 2 2 ,   w 22 = P ( q 2 E + r 2 e ( K + N ) Ω 2 2 ) ,
where Ω 1 = c + N ( 1 δ P ) + e P and Ω 2 = K + N ( 1 δ P ) . Direct computation reveals that J ( W 4 ) has two negative real-part eigenvalues, implying that W 4 is LAS if the following necessary condition is satisfied.
a ( 1 δ P ) 2 P Ω 2 < b + q 1 E .
2 δ P Ω 1 + e a P ( 1 δ P ) Ω 1 2 < r 1 n ( 1 m ) ( n + P ) 2 + a Ω 1 + δ a N P ( 1 δ P ) Ω 1 2 .
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.
r 1 ( m + n ( 1 m ) n + P ˜ ) > d + a ( 1 δ P ˜ ) P ˜ c + e P ˜ .
Proof. 
Consider the function U ( N , P ) = N α 1 P α 2 , where α 1 and α 2 are positive constants.
Obviously, U ( N , P ) > 0 for all ( N , P ) + 2 , and U ( N , P ) = 0 for all ( N , P ) + 2 , where + 2 represents the boundary of + 2 . Hence, U ( N , P ) 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 φ ( N , P ) = U ( N , P ) U ( N , P ) is positive for all points ( N , P ) that belong to the ω limit sets of System (2) in + 2 . The following is the result of direct computation:
φ ( N , P ) = U ( N , P ) U ( N , P ) = α 1 N d N d T + α 2 P d P d T .
Accordingly, W 1 , W 2 , and W 3 belong to the ω limit sets of System (2) in + 2 . Moreover, the direct calculation gives that:
φ ( W 1 ) = α 1 ( r 1 d ) + α 2 r 2 . φ ( W 2 ) = α 2 r 2 . φ ( W 3 ) = α 1 ( r 1 ( m + n ( 1 m ) n + P ˜ ) d a ( 1 δ P ˜ ) P ˜ c + e P ˜ ) .
Clearly, for a suitable choice of α 1 and α 2 , so that α 1 is sufficiently larger than α 2 , it is obtained that φ ( W 1 ) > 0 , while φ ( W 2 ) > 0 always. However, φ ( W 3 ) > 0 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.
r 1 < d .
( r 2 e P ˜ ( 1 δ P ) K Ω 2 ) 2 < 4 ( b + q 2 E ) ( r 2 e Ω 2 + q 2 E ) .
Proof. 
Define the function V 1 ( N , P ) = N + [ P P ˜ P ˜ l n ( P P ˜ ) ] , which is a positive definite real-valued function on the region D 1 = { ( N , P ) + 2 : N 0 ,   P > 0 } . Then, after performing some direct calculations, it is deduced that
d V 1 d T = d N d T + ( P P ˜ P ) d P d T ( r 1 d ) N ( b + q 1 E ) N 2 a ( 1 δ P ) N P Ω 1 ( r 2 e Ω 2 + q 2 E ) ( P P ˜ ) 2 + ( r 2 e P ˜ ( 1 δ P ) K Ω 2 ) N ( P P ˜ ) .
Further simplification leads to:
d V 1 d T ( r 1 d ) N [ b + q 1 E N r 2 e Ω 2 + q 2 E ( P P ˜ ) ] 2 .
Obviously, d V 1 d T 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 D ( N , P ) (called the Dulac function) such that the expression Δ = N ( D d N d T ) + P ( D d P d T ) has the same sign ( 0 ) 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 D ( N , P ) = 1 N P , which is a C 1 in a simply connected region of the domain of System (2), then direct computation with respect to System (2) shows that
Δ = b q 1 E P + a ( 1 + δ P ) 2 Ω 1 2 r 2 e N ( K + N ) Ω 2 2 q 2 E N .
Clearly, Δ < 0 due to Condition (12) of the local stability of the SEP.
Thus, System (2) has no periodic dynamics in the interior of + 2 , 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.
Rewrite System (2) as:
d Z d T = F ( Z ) ,   with   Z = ( N P ) ,   and   F = ( N f ( N , P ) P g ( N , P ) ) .
Hence, the second directional derivative of F , where V = ( v 1 , v 2 ) T is any vector with μ + , can be written using direct computation as:
D 2 F ( Z , μ ) . ( V , V ) = ( c 11 c 21 ) ,
where:
c 11 = c a c + N Ω 1 3 v 1 v 2 + 2 a c + N c δ + e N Ω 1 3 v 2 2 + a 2 c 2 δ c e + c δ N 2 e N P Ω 1 3 v 1 v 2 + 3 a δ e c + N P 2 Ω 1 3 v 1 v 2 + a δ e e δ N P 3 Ω 1 3 v 1 v 2 + 2 1 m n r 1 N n + P 3 v 2 2 + 1 m n r 1 n + P 2 v 1 v 2 + a N 1 δ P c 3 c δ P + x 1 + δ P 2 e P 1 + δ P Ω 1 3 v 1 v 2 2 n + P 2 a N P 1 δ P 3 a P 1 + δ P 2 Ω 1 + b + E q 1 Ω 1 3 n + P 2 Ω 1 3 v 1 2 Ω 1 1 m n r 1 Ω 1 2 + a n + P 2 c + x 2 δ c + N P + δ e + δ N P 2 n + y 2 Ω 1 3 v 1 v 2
c 21 = 2 E q 2 v 2 2 2 e r 2 K + N 2 Ω 2 3 v 2 2 + 4 e r 2 K + N P Ω 2 3 v 1 v 2 2 e r 2 P 2 Ω 2 3 v 1 2   + 4 e r 2 δ P 3 Ω 2 3 v 1 2 6 e r 2 δ v 1 v 2 K + N P 2 Ω 2 3 + 2 e r 2 δ 2 N P 3 Ω 2 3 v 1 v 2 2 e r 2 δ 2 P 4 Ω 2 3 v 1 2
Theorem 6.
System (2) undergoes a transcritical bifurcation (TB) at the PYFEP when the parameter d  passes through the value d = r 1 ( m + n ( 1 m ) n + P ˜ ) a ( 1 δ P ˜ ) P ˜ c + e P ˜ , provided that the following condition holds:
c ˜ 11 0 ,
 where the symbol  c ˜ 11   is determined in the proof.
Proof. 
From the JM that is given by Equation (9), it is observed that for d = d , it becomes
J 1 = J ( W 3 ,   d ) = ( 0 0 r 2 e P ˜ 2 ( 1 δ P ˜ ) K 2 r 2 ) = ( d i j ) .
Therefore, J 1 has eigenvalues given by λ 31 = 0 and λ 32 = r 2 . Thus, W 3 becomes a non-hyperbolic point. Let V 1 = ( v 11 v 21 ) and U 1 = ( u 11 u 21 ) be the eigenvectors corresponding to the λ 31 = 0 of J 1 and their transpose, respectively. Then, the direct calculation gives that:
V 1 = ( 1 e P ˜ 2 ( 1 δ P ˜ ) K 2 ) = ( 1 ϵ 1 ) ,   and   U 1 = ( 1 0 ) .
Moreover, simple computation gives that:
F d ( Z , d ) = ( N 0 ) F d ( W 3 , d ) = ( 0 0 ) . Hence, it is obtained that U 1 Τ F d ( W 3 , d ) = 0 .
U 1 Τ [ D F d ( W 3 , d ) V 1 ] = 1 0 .
Moreover, according to Equation (17), the following result is obtained
D 2 F ( W 3 , d ) ( V 1 , V 1 ) = ( c ˜ 11 c ˜ 21 ) ,
where:
c ˜ 11 = ( 1 m ) n r 1 ( n + P ˜ ) 2 ϵ 1 a c 2 ( c + e P ˜ ) 3 ϵ 1 + a ( 2 c 2 δ c e ) P ˜ ( c + e P ˜ ) 3 ϵ 1 + 3 a c δ e P ˜ 2 ( c + e P ˜ ) 3 ϵ 1 + a δ e 2 P ˜ 3 ( c + e P ˜ ) 3 ϵ 1 2 [ a P ˜ ( 1 δ P ˜ ) 2 + ( b + E q 1 ) ( c + e P ˜ ) 2 ] ( c + e P ˜ ) 2 [ ( 1 m ) n r 1 ( c + e P ˜ ) 2 + a ( n + P ˜ ) 2 ( c 2 c δ P ˜ δ e P ˜ 2 ) ] ( n + P ˜ ) 2   ( c + e P ˜ ) 2 ϵ 1 . c ˜ 21 = 2 E q 2 ϵ 1 2 2 e r 2 K ϵ 1 2 + 4 e r 2 P ˜ K 2 ϵ 1 2 e r 2 P ˜ 2 K 3 6 e r 2 δ P ˜ 2 K 2 ϵ 1 + 4 e r 2 δ P ˜ 3 K 3 2 e r 2 δ 2 P ˜ 4 K 3 .
Hence, it is obtained that U 1 Τ [ D 2 F ( W 3 , d ) ( V 1 , V 1 ) ] = c ˜ 11 , 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 W 3 , and this completes the proof. □
Theorem 7.
If Condition (12) is met along with the following condition:
r 1 n ( 1 m ) ( n + P ) 2 + a Ω 1 + δ a N P ( 1 δ P ) Ω 1 2 < 2 δ P Ω 1 + e a P ( 1 δ P ) Ω 1 2 .
 then, as the parameter r 1 = r 1 , System (2) presents a saddle-node bifurcation (SNB) near the SEP provided that the following condition holds.
ϵ 3 c 11 + c 21 0 ,
 where:
r 1 = ( n + P ) 2 n ( 1 m ) ( 2 δ P Ω 1 + e a P ( 1 δ P ) Ω 1 2 a Ω 1 + δ a N P ( 1 δ P ) Ω 1 2 w 11 w 22 N w 21 ) .
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 r 1 = r 1 , then it can be written as:
J 2 = J ( W 4 , r 1   ) = ( w 11 w 12 w 21 w 22 ) ,
where w 11 ,   w 12 , w 21 , and w 22 are the elements of J ( W 4 ) , with w 12 = w 12 ( r 1 ) . Direct calculation shows that the determinant of J 2 is zero.
Then J 2 has a zero eigenvalue ( λ 41 = 0 ) with the second eigenvalue λ 42 = w 11 + w 22 < 0 under Condition (12). Thus, the SEP is a non-hyperbolic point when r 1 = r 1 .
Let V 2 = ( v 12 v 22 ) and U 2 = ( u 12 u 22 ) be the eigenvectors associated with λ 41 = 0 of J 2 and their transpose, respectively. Then, the direct calculation gives that:
V 2 = ( w 12 w 11 1 ) = ( ϵ 2 1 ) ,   and   U 2 = ( w 21 w 11 1 ) = ( ϵ 3 1 ) .
Note that, according to the Jacobian elements, we obtain that ϵ 2 > 0 and ϵ 3 > 0 due to Conditions (12) and (19). Moreover, simple computation gives that:
F r 1 ( Z , r 1 ) = ( m N + n ( 1 m ) N n + P 0 ) F r 1 ( W 4 , r 1 ) = ( m N + n ( 1 m ) N n + P 0 )
Hence, It is concluded that U 2 Τ F r 1 ( Z 4 , r 1 ) = [ m + n ( 1 m ) n + P ] N ϵ 3 0 . Moreover, according to Equation (17), the following is gained:
D 2 F ( W 4 , r 1 ) ( V 2 , V 2 ) = ( c 11 c 21 ) = ( c 11 ( W 4 , r 1 , V 2 ) c 21 ( W 4 , r 1 , V 2 ) ) .
Therefore, the expression U 2 Τ [ D 2 F ( W 4 , r 1 ) ( V 2 , V 2 ) ] = ϵ 3 c 11 + c 21 , 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:
b + q 1 E < a ( 1 δ P ) 2 P Ω 2 ,
w 11 w 22 w 12 w 21 > 0 ,
where  w 11 , w 12 ,   w 21 , and w 22 are the JM elements that are given in Equation (11). Then, System (2) presents a Hopf bifurcation around the SEP when the parameter a passes through the value a with
a = Ω 1 2 ( 1 δ P ) 2 N P [ q 2 E P + r 2 e ( K + N ) P Ω 2 2 + ( b + q 1 E ) N ] .
Proof. 
Direct computation shows that the characteristic equation of the JM of System (2) at SEP can be determined as:
λ 2 ( w 11 + w 22 ) λ + ( w 11 w 22 w 12 w 21 ) = 0 .
Note that at a = a , it is easy to verify that w 11 + w 22 = 0 , due to Condition (21). Consequently, in this case, the roots of Equation (23) are written as:
λ 1 , 2 = ± i w 11 w 22 w 12 w 21 ,
where w 11 = w 11 ( a ) and w 12 = w 12 ( a ) . Clearly, the above eigenvalues are complex conjugate pure imaginary under Condition (22). However, in the neighborhood of a = a , Equation (23) has two complex conjugate eigenvalues given by:
λ 1 , 2 = w 11 + w 22 2 ± ( w 11 + w 22 ) 2 4 ( w 11 w 22 w 12 w 21 ) 2 .
Now, since
d d a [ R e ( λ 1 , 2 ) ] a = a = d d a [ w 11 + w 22 2 ] a = a = ( 1 δ P ) 2 N P Ω 1 2 0 ,
therefore, according to the Hopf bifurcation theorem, System (2) has a Hopf bifurcation at a = a , 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).
r 1 = 2 ,   m = 0.5 ,   n = 1 ,   d = 0.05 ,   b = 0.1 ,   a = 0.75 ,   c = 2 ,   e = 0.2 ,   δ = 0.05 ,   q 1 = 0.1 ,   q 2 = 0.2 ,   E = 0.75 ,   r 2 = 1 ,   K = 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 r 1 2.24 ; for instance, look at Figure 2 below at r 1 = 3 . System (2) does not have a SEP; however, for the range r 1 1.5 , 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 r 2 on the dynamics of System (2) is presented in Figure 4, for the values r 2 = 0.5 ,   1 ,   1.5 ,   2 ,   2.5 ,   3 , keeping the rest of the parameters as in Dataset (24).
As shown in Figure 4, the population of N decreases while that of P increases with the increase of r 2 . 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 m 0.31 , System (2) has no SEP and the system approaches the PYFEP; see Figure 5; however, varying the parameter n , such that n = 0.5 ,   1.5 ,   2.5 ,   3.5 ,   4.5 ,   5.5 with the other parameters as in (24), is investigated in Figure 6.
As shown in Figure 6, the populations of N and P increase with the increase of n . Now, varying the parameter d , such that d = 0.01 ,   0.11 ,   0.21 ,   0.31 ,   0.41 , 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 b 0.08 . 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 d causes the decrease in both species N and P , with N decreasing more and faster. Moreover, it is observed that the increase in the value of a so that a 1.03 has a similar influence as that with d on the dynamics of System (2). On the other hand, System (2) has a unique SEP that is GAS in its state-space for 0 < a < 1.03 .
Now that the parameter c has been varied in the ranges of c 1.02 , 1.03 c 1.19 , and 1.20 c , 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 e has been varied in the ranges of e 0.08 , c = 0.09 , and 0.1 c , 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 e leads to an increase in the population of N and a decrease in P . The influences of changing the harvesting parameters q 1 ,   q 2 , and E are investigated and presented in Figure 11.
It is clear from Figure 11 that q 1 and q 2 have the opposite influence on the populations of System (2), while both populations decrease as E increases.
For different values of K 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.

Author Contributions

The authors A.R.M.J. and R.K.N. contributed to the conception, design of the model, acquisition of the data, analysis, interpretation, drafting of the MS, revision, and proofreading. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors are thankful to the anonymous Reviewers for their careful reading, useful comments, and constructive suggestions for the improvement of the present research work. We are also thankful to the Editor for his speed and cooperation.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Holling, C.S. The components of predation as revealed by a study of small mammal predation of the european pine sawfly. Can. Entomol. 1959, 91, 293–320. [Google Scholar] [CrossRef]
  2. Berryman, A.A. The origins and evolution of predator-prey theory. Ecology 1992, 75, 1530–1535. [Google Scholar] [CrossRef] [Green Version]
  3. Ruan, S.; Xiao, D. Global analysis in a predator-prey system with nonmonotonic functional response. SIAM J. Appl. Math. 2001, 61, 1445–1472. [Google Scholar]
  4. Aziz-Alaoui, M.A.; Okiye, M.D. Boundedness and global stability for a predator-prey model with modified Leslie-Gower and Holling-type II schemes. Appl. Math. Lett. 2003, 16, 1069–1075. [Google Scholar] [CrossRef] [Green Version]
  5. Wang, X.; Zanette, L.; Zou, X. Modelling the fear effect in predator-prey interactions. J. Math. Biol. 2016, 73, 1179–1204. [Google Scholar] [CrossRef]
  6. Wang, X.; Zou, X. Modeling the Fear Effect in Predator-Prey Interactions with Adaptive Avoidance of Predators. Bull. Math. Biol. 2017, 79, 1325–1359. [Google Scholar] [CrossRef]
  7. Elliott, K.H.; Betini, G.S.; Norris, D.R. Fear creates an Allee effect: Experimental evidence from seasonal populations. Proc. R. Soc. B Biol. Sci. 2017, 284, 20170878. [Google Scholar] [CrossRef]
  8. Pal, S.; Majhi, S.; Mandal, S.; Pal, N. Role of Fear in a Predator-Prey Model with Beddington–DeAngelis Functional Response. Z. Für Nat. A 2019, 74, 581–595. [Google Scholar] [CrossRef]
  9. Fakhry, N.H.; Naji, R.K. The Dynamics of a Square Root Prey-Predator Model with Fear. Iraqi J. Sci. 2020, 61, 139–146. [Google Scholar] [CrossRef]
  10. Sarkar, K.; Khajanchi, S. Impact of fear effect on the growth of prey in a predator-prey interaction model. Ecol. Complex. 2020, 42, 100826. [Google Scholar] [CrossRef]
  11. Sasmal, S.K.; Takeuchi, Y. Dynamics of a predator-prey system with fear and group defense. J. Math. Anal. Appl. 2020, 481, 123471. [Google Scholar] [CrossRef]
  12. Maghool, F.H.; Naji, R.K. The dynamics of a tritrophic Leslie-Gower food-web system with the effect of fear. J. Appl. Math. 2021, 2021, 2112814. [Google Scholar] [CrossRef]
  13. Liu, J.; Lv, P.; Liu, B.; Zhang, T. Dynamics of a Predator-Prey Model with Fear Effect and Time Delay. Complexity 2021, 2021, 9184193. [Google Scholar] [CrossRef]
  14. Vinoth, S.; Sivasamy, R.; Sathiyanathan, K.; Unyong, B.; Rajchakit, G.; Vadivel, R.; Gunasekaran, N. The dynamics of a Leslie type predator–prey model with fear and Allee effect. Adv. Differ. Equ. 2021, 2021, 338. [Google Scholar] [CrossRef]
  15. Maghool, F.H.; Naji, R.K. Chaos in the three-species Sokol-Howell food chain system with fear. Commun. Math. Biol. Neurosci. 2022, 2022, 14. [Google Scholar] [CrossRef]
  16. Rahi, S.A.; Kurnaz, S.; Naji, R.K. The impact of fear on a stage structure prey-predator system with anti-predator behavior. Appl. Nanosci. 2022. [Google Scholar] [CrossRef]
  17. Sarwardi, S.; Mandal, P.K.; Ray, S. Analysis of a competitive prey-predator system with a prey refuge. Biosystems 2012, 110, 133–148. [Google Scholar] [CrossRef]
  18. Abdulghafour, A.S.; Naji, R.K. A Study of a Diseased Prey-Predator Model with Refuge in Prey and Harvesting from Predator. J. Appl. Math. 2018, 2018, 2952791. [Google Scholar] [CrossRef]
  19. González-Olivares, E.; González-Yañez, B.; Becerra-Klix, R.; Ramos-Jiliberto, R. Multiple stable states in a model based on predator-induced defenses. Ecol. Complex. 2017, 32, 111–120. [Google Scholar] [CrossRef]
  20. Manarul Haque, M.; Sarwardi, S. Dynamics of a Harvested Prey–Predator Model with Prey Refuge Dependent on Both Species. Int. J. Bifurc. Chaos Appl. Sci. Eng. 2018, 28, 1830040. [Google Scholar] [CrossRef] [Green Version]
  21. Molla, H.; Rahman, M.S.; Sarwardi, S. Dynamics of a predator-prey model with Holling Type II Functional Response Incorporating a Prey Refuge Depending on Both the Species. Int. J. Nonlinear Sci. Numer. Simul. 2019, 20, 89–104. [Google Scholar] [CrossRef]
  22. Ji, L.; Wu, C. Qualitative analysis of a predator-prey model with constant-rate prey harvesting incorporating a constant prey refuge. Nonlinear Anal. Real World Appl. 2010, 11, 2285–2295. [Google Scholar] [CrossRef]
  23. Pei, Y.; Lv, Y.; Li, C. Evolutionary consequences of harvesting for a two-zooplankton one-phytoplankton system. Appl. Math. Model. 2012, 36, 1752–1765. [Google Scholar] [CrossRef] [Green Version]
  24. Lv, Y.; Yuan, R.; Pei, Y. Dynamics in two nonsmooth predator-prey models with threshold harvesting. Nonlinear Dyn. 2013, 74, 107–132. [Google Scholar] [CrossRef]
  25. Gupta, R.P.; Chandra, P. Bifurcation analysis of modified Leslie–Gower predator-prey model with Michaelis–Menten type prey harvesting. J. Math. Anal. Appl. 2013, 398, 278–295. [Google Scholar] [CrossRef]
  26. Al-Momen, S.; Naji, R.K. The Dynamics of Modified Leslie-Gower Predator-Prey Model under the Influence of Nonlinear Harvesting and Fear Effect. Iraqi J. Sci. 2022, 63, 259–282. [Google Scholar] [CrossRef]
  27. Yu, S. Global stability of a modified Leslie-Gower model with Beddington-DeAngelis functional response. Adv. Differ. Equ. 2014, 2014, 84. [Google Scholar] [CrossRef] [Green Version]
  28. Gard, T.C. Uniform persistence in multispecies population models Author links open overlay panel. Math. Biosci. 1987, 85, 93–104. [Google Scholar] [CrossRef]
  29. Perko, L. Differential Equations and Dynamical Systems, 3rd ed.; Springer: New York, NY, USA, 2001. [Google Scholar]
Figure 1. The phase portrait of System (2) using Dataset (24) and starting from different points. (a) The trajectories approach asymptotically the SEP given by W 4 = ( 3.41 ,   5.15 ) . (b) The time series of the trajectories given by the phase portrait.
Figure 1. The phase portrait of System (2) using Dataset (24) and starting from different points. (a) The trajectories approach asymptotically the SEP given by W 4 = ( 3.41 ,   5.15 ) . (b) The time series of the trajectories given by the phase portrait.
Mathematics 10 02857 g001
Figure 2. The phase portrait of System (2) using Dataset (24) with r 1 = 3 and starting from different points. (a) The trajectories approach asymptotically two different SEPs depending on their initial points. (b) The time series of the trajectories given by the phase portrait.
Figure 2. The phase portrait of System (2) using Dataset (24) with r 1 = 3 and starting from different points. (a) The trajectories approach asymptotically two different SEPs depending on their initial points. (b) The time series of the trajectories given by the phase portrait.
Mathematics 10 02857 g002
Figure 3. The phase portrait of System (2) using Dataset (24) with r 1 = 1.4 and starting from different points. (a) The trajectories approach asymptotically the PYFEP, W 3 = ( 0.4 ) . (b) The time series of the trajectories given by the phase portrait.
Figure 3. The phase portrait of System (2) using Dataset (24) with r 1 = 1.4 and starting from different points. (a) The trajectories approach asymptotically the PYFEP, W 3 = ( 0.4 ) . (b) The time series of the trajectories given by the phase portrait.
Mathematics 10 02857 g003
Figure 4. The phase portrait of System (2) using Dataset (24) with r 2 = 0.5 ,   1 ,   1.5 ,   2 ,   2.5 ,   3  (a) The trajectories approach asymptotically the different SEPs depending on the value of r 2 . (b) Time series of N ‘s trajectories as given by the phase portrait. (c) Time series of P ’s trajectories as given by the phase portrait.
Figure 4. The phase portrait of System (2) using Dataset (24) with r 2 = 0.5 ,   1 ,   1.5 ,   2 ,   2.5 ,   3  (a) The trajectories approach asymptotically the different SEPs depending on the value of r 2 . (b) Time series of N ‘s trajectories as given by the phase portrait. (c) Time series of P ’s trajectories as given by the phase portrait.
Mathematics 10 02857 g004
Figure 5. The phase portrait of System (2) using Dataset (24) with m = 0.2  and starting from different points. (a) The trajectories approach asymptotically the PYFEP, W 3 = ( 0, 4 ) . (b) The time series of the trajectories given by the phase portrait.
Figure 5. The phase portrait of System (2) using Dataset (24) with m = 0.2  and starting from different points. (a) The trajectories approach asymptotically the PYFEP, W 3 = ( 0, 4 ) . (b) The time series of the trajectories given by the phase portrait.
Mathematics 10 02857 g005
Figure 6. The phase portrait of System (2) using Dataset (24) with n = 0.5 ,   1.5 ,   2.5 ,   3.5 ,   4.5 ,   5.5 . (a) The trajectories approach asymptotically the different SEPs depending on the value of n. (b) Time series of N‘s trajectories as given by the phase portrait. (c) Time series of P’s trajectories as given by the phase portrait.
Figure 6. The phase portrait of System (2) using Dataset (24) with n = 0.5 ,   1.5 ,   2.5 ,   3.5 ,   4.5 ,   5.5 . (a) The trajectories approach asymptotically the different SEPs depending on the value of n. (b) Time series of N‘s trajectories as given by the phase portrait. (c) Time series of P’s trajectories as given by the phase portrait.
Mathematics 10 02857 g006
Figure 7. The phase portrait of System (2) using Dataset (24) with d = 0.01 ,   0.11 ,   0.21 ,   0.31 ,   0.41 . (a) The trajectories approach asymptotically the different SEPs depending on the value of d. (b) Time series of N’s trajectories as given by the phase portrait. (c) Time series of P’s trajectories as given by the phase portrait.
Figure 7. The phase portrait of System (2) using Dataset (24) with d = 0.01 ,   0.11 ,   0.21 ,   0.31 ,   0.41 . (a) The trajectories approach asymptotically the different SEPs depending on the value of d. (b) Time series of N’s trajectories as given by the phase portrait. (c) Time series of P’s trajectories as given by the phase portrait.
Mathematics 10 02857 g007
Figure 8. The phase portrait of System (2) using Dataset (24) with b = 0.01 and starting from different points. (a) The trajectories approach asymptotically two different SEPs depending on their initial points. (b) The time series of the trajectories given by the phase portrait.
Figure 8. The phase portrait of System (2) using Dataset (24) with b = 0.01 and starting from different points. (a) The trajectories approach asymptotically two different SEPs depending on their initial points. (b) The time series of the trajectories given by the phase portrait.
Mathematics 10 02857 g008
Figure 9. The phase portrait of System (2) using Dataset (24) with c and starting from different points. (a) The trajectories approach asymptotically the PYFEP globally. (b) The trajectories approach asymptotically the PYFEP and SEP locally (bi-stable). (c) The trajectories approach asymptotically the SEP globally.
Figure 9. The phase portrait of System (2) using Dataset (24) with c and starting from different points. (a) The trajectories approach asymptotically the PYFEP globally. (b) The trajectories approach asymptotically the PYFEP and SEP locally (bi-stable). (c) The trajectories approach asymptotically the SEP globally.
Mathematics 10 02857 g009
Figure 10. The phase portrait of System (2) using Dataset (24) with varying e and starting from different points. (a) The trajectories approach asymptotically the PYFEP globally for e = 0.08 . (b) The trajectories approach asymptotically the PYFEP and SEP locally (bi-stable) for e = 0.09 . (c) The trajectories approach asymptotically the SEP globally for e = 0.1 . (d) The trajectories approach asymptotically another SEP globally for e = 0.05 .
Figure 10. The phase portrait of System (2) using Dataset (24) with varying e and starting from different points. (a) The trajectories approach asymptotically the PYFEP globally for e = 0.08 . (b) The trajectories approach asymptotically the PYFEP and SEP locally (bi-stable) for e = 0.09 . (c) The trajectories approach asymptotically the SEP globally for e = 0.1 . (d) The trajectories approach asymptotically another SEP globally for e = 0.05 .
Mathematics 10 02857 g010
Figure 11. The phase portrait of System (2) using Dataset (24). (a) The trajectories approach asymptotically the different SEPs for q 1 = 0.01 ,   0.21 ,   0.41 ,   0.61 ,   0.81 . (b) The trajectories approach asymptotically the different SEPs for q 2 = 0.01 ,   0.21 ,   0.41 ,   0.61 ,   0.81 . (c) The trajectories approach asymptotically the different SEPs for E = 0.5 ,   1.5 ,   2.5 ,   3.5 ,   4.5 ,   5.5 .
Figure 11. The phase portrait of System (2) using Dataset (24). (a) The trajectories approach asymptotically the different SEPs for q 1 = 0.01 ,   0.21 ,   0.41 ,   0.61 ,   0.81 . (b) The trajectories approach asymptotically the different SEPs for q 2 = 0.01 ,   0.21 ,   0.41 ,   0.61 ,   0.81 . (c) The trajectories approach asymptotically the different SEPs for E = 0.5 ,   1.5 ,   2.5 ,   3.5 ,   4.5 ,   5.5 .
Mathematics 10 02857 g011aMathematics 10 02857 g011b
Figure 12. The phase portrait of System (2) using Dataset (24) with K = 0.5 ,   1.5 ,   2.5 ,   3.5 ,   4.5 ,   5.5 . (a) The trajectories approach asymptotically the different SEPs depending on the value of K . (b) The time series of the trajectories given by the phase portrait, which indicates decreasing in N and increasing in P as K increases.
Figure 12. The phase portrait of System (2) using Dataset (24) with K = 0.5 ,   1.5 ,   2.5 ,   3.5 ,   4.5 ,   5.5 . (a) The trajectories approach asymptotically the different SEPs depending on the value of K . (b) The time series of the trajectories given by the phase portrait, which indicates decreasing in N and increasing in P as K increases.
Mathematics 10 02857 g012
Figure 13. The phase portrait of System (2) using Dataset (24) with varying δ and starting from different points. The plotted arrows show the direction of the trajectories. (a) The trajectories approach asymptotically a SEP globally for δ = 0.01 . (b) The trajectories approach asymptotically another SEP globally for δ = 0.03 . (c) The trajectories approach asymptotically two different SEP locally (bi-stable) for δ = 0.06 . (d) The trajectories approach asymptotically two different SEP locally (bi-stable) for δ = 0.07 .
Figure 13. The phase portrait of System (2) using Dataset (24) with varying δ and starting from different points. The plotted arrows show the direction of the trajectories. (a) The trajectories approach asymptotically a SEP globally for δ = 0.01 . (b) The trajectories approach asymptotically another SEP globally for δ = 0.03 . (c) The trajectories approach asymptotically two different SEP locally (bi-stable) for δ = 0.06 . (d) The trajectories approach asymptotically two different SEP locally (bi-stable) for δ = 0.07 .
Mathematics 10 02857 g013
Table 1. Parameters’ description.
Table 1. Parameters’ description.
ParameterDescription
r 1 ,   r 2 The birth rate of prey population and predator population, respectively.
m The minimum cost of fear with m [ 0 , 1 ] .
n The level of fear.
d The natural death rate of the prey.
b Decay rate due to intraspecific competition.
a Attack rate.
c Half saturation constant.
e A measure of the food quantity that the prey provides converted to predator birth.
q 1 ,   q 2 The catchability coefficients of the prey and predator, respectively.
E The effort level for harvesting the prey and predator.
K The carrying capacity of the predator in the absence of its prey.
δ The refuge coefficient.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Jamil, A.R.M.; Naji, R.K. Modeling and Analysis of the Influence of Fear on the Harvested Modified Leslie–Gower Model Involving Nonlinear Prey Refuge. Mathematics 2022, 10, 2857. https://doi.org/10.3390/math10162857

AMA Style

Jamil ARM, Naji RK. Modeling and Analysis of the Influence of Fear on the Harvested Modified Leslie–Gower Model Involving Nonlinear Prey Refuge. Mathematics. 2022; 10(16):2857. https://doi.org/10.3390/math10162857

Chicago/Turabian Style

Jamil, Abdul Rahman Mahmoud, and Raid Kamel Naji. 2022. "Modeling and Analysis of the Influence of Fear on the Harvested Modified Leslie–Gower Model Involving Nonlinear Prey Refuge" Mathematics 10, no. 16: 2857. https://doi.org/10.3390/math10162857

APA Style

Jamil, A. R. M., & Naji, R. K. (2022). Modeling and Analysis of the Influence of Fear on the Harvested Modified Leslie–Gower Model Involving Nonlinear Prey Refuge. Mathematics, 10(16), 2857. https://doi.org/10.3390/math10162857

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop