Next Article in Journal
Functional Subspace Variational Autoencoder for Domain-Adaptive Fault Diagnosis
Previous Article in Journal
Multimodal Prompt Learning in Emotion Recognition Using Context and Audio Information
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The Impact of Fear on a Harvested Prey–Predator System with Disease in a Prey

by
Hiba Abdullah Ibrahim
and
Raid Kamel Naji
*
Department of Mathematics, College of Science, University of Baghdad, Baghdad 10071, Iraq
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(13), 2909; https://doi.org/10.3390/math11132909
Submission received: 3 June 2023 / Revised: 22 June 2023 / Accepted: 26 June 2023 / Published: 28 June 2023
(This article belongs to the Section Mathematical Biology)

Abstract

:
A mathematical eco-epidemiological model consisting of harvested prey–predator system involving fear and disease in the prey population is formulated and studied. The prey population is supposed to be separated into two groups: susceptible and infected. The susceptible prey grows logistically, whereas the infected prey cannot reproduce and instead competes for the environment’s carrying capacity. Furthermore, the disease is transferred through contact from infected to susceptible individuals, and there is no inherited transmission. The existence, positivity, and boundedness of the model’s solution are discussed. The local stability analysis is carried out. The persistence requirements are established. The global behavior of the system is investigated with the use of the Lyapunov method. An application to the Sotomoyar theorem of local bifurcation is performed around the equilibrium points. In the end, the system is numerically simulated to confirm our obtained analytical results and specify the control set of parameters. Bifurcation diagrams are used to show the dynamical behavior as a function of some parameters. It is obtained that the prey’s fear stabilizes the system, while the disease and harvest cause extinction in one or more species.
MSC:
92D40; 92D30; 34D20; 34C23

1. Introduction

Due to its widespread existence and significance, the dynamic interaction between prey and predator has been the subject of research and will continue to be a crucial topic; see [1,2,3,4,5,6,7,8,9] and the references cited therein. The well-known Lotka–Volterra type prey–predator system has been the focus of numerous good articles since being raised by Lotka and Volterra [10,11]. A Canadian scientist by the name of Holling [12] introduced the comparable functional response functions for diverse species in 1959, including type I, type II, and type III, based on his experimental results, which comprised three fundamental types of predation rates. Since weak (or slow) prey is relevant to Holling type I functional response function, the rate of prey consumption rises linearly. In contrast, type II functional response shows that a predator’s rate of prey intake grows as prey density increases but eventually levels off at a particular point, remaining constant regardless of increases in prey density. The type III functional response, on the other hand, is comparable to the type II response in that saturation happens at high prey density levels. However, at low prey density levels, the graph between the number of prey eaten and the size of the prey population is now a strictly increasing function of prey consumed by predators, making it suitable for the prey with a group defense capability.
The consistent environment is unusual in natural existence. Fear, stage structure, delay, harvesting, toxicity, cannibalism, anti-predator capabilities, refuge, infective disease, and other important factors affecting the population all have an impact on the natural environment. The principal biotic threats to the prey–predator relationship are infectious disease, which causes both species to suffer and has economic consequences, harvesting, and fear, which lowers prey growth. As a result, knowing the impact of biotic and abiotic forces on the ecosystem is critical. Consequently, it is worth looking into the prey–predator system in the context of fear, harvesting, and illness. For many years, mathematical biologists have attempted to integrate the two primary research areas, ecology and epidemiology. Eco-epidemiological refers to the integration of infection into the ecological model [13]. Eco-epidemiological models, which are described as mathematical representations of ecological systems where the disease impacts prey and/or predators, have been the subject of many investigations. When analyzing the effects of disease on an environment, disease in the predator, disease in the prey, and disease in both groups can all be taken into consideration. Predators may eat susceptible and ill animals if the disease is present in the prey. Studies on disease, infection, prey–predator dynamics, and possible ecological and biological implications are mentioned in [14,15,16,17,18,19,20,21,22,23,24]. On the other hand, the impact of harvesting a natural population is one of the most important issues in population ecology. According to Ganguli and Kar [25] and the references therein, harvesting can signify a population decline because of hunting or capturing individuals and removing them from the population. Numerous studies have looked at the dynamic behaviors of harvesting populations because there is so much interest in using bio-economic models; see, for example, [19,21,25,26,27]. This component is also useful for population control [28]. Most of the time, achieving a sizable harvest gain from the population is the aim in addition to population management.
Additionally, fear has a remarkable impact on ecosystems. The impact of terror on prey species was a topic of study for many scholars. The effects of fear can be either direct or indirect; in the former, the predator directly kills the prey [29]. In contrast, in the indirect example, a predator’s presence drastically alters the prey’s behavior due to their dread of the threat of predation [30]. As a result, numerous studies have been conducted to comprehend how prey fear affects the dynamics of the prey–predator system, see [9,31,32,33,34,35,36,37,38,39] and the references therein. Despite the aforementioned, this paper combines the effects of fear, disease, and harvest in a single prey–predator system. Moreover, this paper is structured as follows: In Section 2, we develop a harvested prey–predator model with fear and sickness in prey species and then investigate the existence, boundedness, and non-negativity of solutions. In Section 3, we addressed the equilibria’s stability and persistence. Section 4 focuses on global dynamics. Section 5 investigates the local bifurcation. Section 6 has some numerical simulations, while Section 5 contains concluding observations.

2. Model Construction

In this section, a real-world eco-epidemiological problem is simulated mathematically depending on the following assumptions. First, let X ( t ) and Y ( t ) be the density of the prey and predator species at time t, respectively. Consider now the assumptions:
  • The prey species X ( t ) is divided into two compartments due to the existence of the disease in their population: Susceptible prey is denoted by S ( t ) and Infected prey is represented by I ( t ) .
  • In the absence of predation, the susceptible prey grows logistically, while the infected prey cannot reproduce; rather, it competes for the carrying capacity of the environment. According to the mass action principle, the disease is transmitted from the infected to susceptible individuals by contact. Finally, the existence of disease causes a specific death rate in the infected population.
  • The susceptible prey growth reduces due to the influence of the fear of the predation process by the predator.
  • The predator species does not distinguish between the susceptible and infected prey and attack the available prey according to the Holling type-II and Lotka–Volttera functional responses for the susceptible and infected prey, respectively.
  • There is an external force that imposes harvesting on the system.
Now, according to the above assumptions, the eco-epidemiological model can be represented as follows:
d S d t = r ( θ + μ ( 1 θ ) μ + Y ) S α S ( S + I ) β S I a 1 S Y b + S h 1 S , d I d t = β S I a 2 I Y ( d 1 + h 2 ) I ,   d Y d t = e 1 a 1 S Y b + S + e 2 a 2 I Y ( d 2 + h 3 ) Y ,  
with initial condition S ( 0 ) > 0 ,   I ( 0 ) > 0 ,   Y ( 0 ) > 0 . While the parameters are positive constants and are described the Table 1.
Moreover, the solution of system (1) is bounded, as seen in the next theorem. It is easy to confirm that the necessary requirement for the survival of all species in the system (1) is given by:
h 1 < r
Theorem 1. 
The solutions of system (1) starting in  + 3  remain positive for all time.
Proof. 
The solution ( S ( t ) , I ( t ) , Y ( t ) ) of the system (1) exists and is unique on [ 0 , τ ] , where 0 < τ < + [40], since the right-hand side of the system (1) is continuous and locally Lipschitzian on C (space of continuous functions). Now, from the first equation of Equation (1), it is obtained that
S ( t ) = S ( 0 ) exp [ 0 t { r ( θ + μ ( 1 θ ) μ + Y ( u ) ) α ( S ( u ) + I ( u ) ) β I ( u ) a 1 Y ( u ) b + S ( u ) h 1 } d u ] > 0 .
Similarly,
I ( t ) = I ( 0 ) exp [ 0 t { β S ( u ) a 2 Y ( u ) ( d 1 + h 2 ) } d u ] > 0 ,
and
Y ( t ) = Y ( 0 ) exp [ 0 t { e 1 a 1 S ( u ) b + S ( u ) + e 2 a 2 I ( u ) ( d 2 + h 3 ) } d u ] > 0 .
Hence, for all starting points in + 3 , the solution remains positive for all time. □
Theorem 2. 
All solutions of system (1) are uniformly bounded.
Proof: 
From the first equation, it is obtained that
d S d t + α S 2 ( r h 1 ) S
Hence,
lim t sup S ( t ) r h 1 α .
Let M ( t ) = S ( t ) + I ( t ) + Y ( t ) , then applying some calculations gives that
d M d t ( r h 1 ) S ( d 1 + h 2 ) I ( d 2 + h 3 ) Y   = ( r h 1 + 1 ) S S ( d 1 + h 2 ) I ( d 2 + h 3 ) Y ( r h 1 + 1 ) r h 1 α L M = N L M  
where L = m i n { 1 ,   d 1 + h 2 ,   d 2 + h 3 } , and N = ( r h 1 ) 2 α + r h 1 α . Direct computation gives that
M ( t ) N L ( 1 e L t ) + M ( 0 ) e L t
Therefore, for t , it is obtained that M ( t ) N L . Therefore, all of the system (1)’s solutions will reach the region:
D = { ( S , I , Y ) + 3 : S ( t ) r h 1 α ,   0 < M ( t ) N L }
Thus, all the solutions are uniformly bounded. □

3. Stability Analysis with Persistence

In this part, the existence of the equilibrium points ( E . P ) , stability analysis, and persistence of the system (1) are studied. It is spotted that the next equilibrium points exist as follows:
The vanishing equilibrium point ( V . E . P ) , which is denoted by m 0 = ( 0 , 0 , 0 ) always exists.
The axial equilibrium point ( A . E . P ) is denoted by m 1 = ( S ^ , 0 , 0 ) where S ^ = r h 1 α , exists under Condition (2).
The predator-free equilibrium point ( P . F . E . P ) that is denoted by m 2 = ( S ¯ , I ¯ , 0 ) , where
S ¯ = d 1 + h 2 β ,   I ¯ = ( r h 1 ) β α ( d 1 + h 2 ) β ( α + β ) ,
exists if the following condition holds:
α ( d 1 + h 2 ) < ( r h 1 ) β
The disease-free equilibrium point ( D . F . E . P ) that is denoted by m 3 = ( S ˜ , 0 , Y ˜ ) , where
S ˜ = b ( d 2 + h 3 ) e 1 a 1 ( d 2 + h 3 ) ,
while Y ˜ is a positive root of the following second-order polynomial equation
B 1 Y 2 + B 2 Y + B 3 = 0 ,
where
B 1 = [ a 1 e 1 d 2 h 3 ] 2 , B 2 = μ [ d 2 + h 3 a 1 e 1 ] 2 b 2 α e 1 ( d 2 + h 3 ) b e 1 ( r θ h 1 ) [ d 2 + h 3 a 1 e 1 ] . B 3 = b μ e 1 [ b α ( d 2 + h 3 ) ( r h 1 ) ( a 1 e 1 d 2 h 3 ) ] .
Obviously, m 3 exists uniquely in the positive quadrant of the S Y -plane provided that the following requirements are met.
( d 2 + h 3 ) < e 1 a 1 .  
b α ( d 2 + h 3 ) < ( r h 1 ) ( a 1 e 1 d 2 h 3 ) .  
The positive equilibrium point ( P . E . P ) that is denoted by m 4 = ( S , I , Y ) , where
I = ( b + S ) ( d 2 + h 3 ) S a 1 e 1 ( b + S ) a 2 e 2 Y = S β ( d 1 + h 2 ) a 2 } .
While S is a positive root of the following third-order polynomial equation
E 1 S 3 + E 2 S 2 + E 3 S + E 4 = 0 ,
where
  E 1 = α β a 2 e 2 . E 2 = β ( α + β ) ( d 2 a 1 e 1 + h 3 ) + e 2 [ β 2 a 1 + β a 2 ( b α r θ + h 1 ) + α a 2 ( μ a 2 d 1 1 ) ] . E 3 = ( α + β ) [ ( d 2 + h 3 ) ( b β + μ a 2 ) μ a 1 a 2 e 1 + ( d 1 + h 2 ) ( a 1 e 1 d 2 h 3 ) ]   e 2 ( d 1 + h 2 ) [ 2 β a 1 + a 2 ( b α r θ + h 1 ) ] + a 2 e 2 h 1 ( b β + μ a 2 ) + β a 2 e 2 ( μ a 1 b r θ ) + μ a 2 2 e 2 ( b α r ) . E 4 = ( α + β ) ( d 2 + h 3 ) [ b μ a 2 b ( d 1 + h 2 ) ] a 2 e 2 [ b r μ a 2 + ( μ a 1 b r θ ) d 1 ] + a 1 e 2 ( d 1 + h 2 ) 2 + b a 2 e 2 h 1 ( μ a 2 d 1 ) + a 2 e 2 h 2 ( b r θ μ a 1 b h 1 ) .
So by Descartes’ rule of sign, Equation (10) has a unique positive root, and hence, P.E.P exists uniquely in the interior of + 3 provided that the following condition is met
( d 1 + h 2 ) β < S < ( b + S ) ( d 2 + h 3 ) a 1 e 1 ,  
with one set of the following sets of conditions.
E 2 > 0 ,   E 3 0 ,   E 4 < 0 E 2 0 ,   E 3 < 0 ,   E 4 < 0 } .
Now, the local stability of the above equilibrium points of system (1) is studied. The Jacobian matrix at m 0 = ( 0 , 0 , 0 ) can be written as
J ( m 0 ) = [ r h 1 0 0 0 ( d 1 + h 2 ) 0 0 0 ( d 2 + h 3 ) ] .
Obviously, the eigenvalues of J ( m 0 ) are given by:
λ 01 = r h 1 ,   λ 02 = ( d 1 + h 2 ) < 0 ,   λ 03 = ( d 2 + h 3 ) < 0 .
Hence, the V.E.P is locally asymptotically stable (L.A.S) if and only if the survival Condition (2) is reflected ( r < h 1 ) , and otherwise, it is a saddle point.
The Jacobian matrix at m 1 = ( r h 1 α , 0 , 0 ) can be written as
J ( m 1 ) = [ r + h 1 ( α + β ) ( r h 1 ) α ( r ( 1 θ ) ( b α + r h 1 ) + a 1 α μ ) ( r h 1 ) μ α ( b α + r h 1 ) 0 β ( r h 1 ) α ( d 1 + h 2 ) α 0 0 0 a 1 e 1 ( r h 1 ) ( b α + r h 1 ) ( d 2 + h 3 ) ] .  
So, the eigenvalues of J ( m 1 ) are given by
λ 11 = ( r h 1 ) < 0 ,   λ 12 = β ( r h 1 ) α ( d 1 + h 2 ) α ,   λ 13 = a 1 e 1 ( r h 1 ) ( b α + r h 1 ) ( d 2 + h 3 ) .
Therefore, the A.E.P. is a L.A.S provided that the following condition holds:
β ( r h 1 ) < α ( d 1 + h 2 ) a 1 e 1 ( r h 1 ) ( b α + r h 1 ) < ( d 2 + h 3 ) } .  
The Jacobian matrix at m 2 = ( S ¯ , I ¯ , 0 ) can be written as
J ( m 2 ) = [ α S ¯ ( α + β ) S ¯ ( r ( 1 θ ) μ + a 1 b + S ¯ ) S ¯ β I ¯ 0 a 2 I ¯ 0 0 e 1 a 1 S ¯ b + S ¯ + e 2 a 2 I ¯ ( d 2 + h 3 ) ] .  
The characteristic equation of J ( m 2 ) can be determined as follows:
( λ 2 T 2 λ + D 2 )   ( e 1 a 1 S ¯ b + S ¯ + e 2 a 2 I ¯ ( d 2 + h 3 ) λ ) = 0 ,  
where T 2 = α S ¯ < 0 , and D 2 = β ( α + β ) S ¯ I ¯ > 0 .
Since, T 2 < 0 and D 2 > 0 , then the two eigenvalues obtain from the first factor of Equation (19) have negative real parts. While the third eigenvalue is given by λ 23 = e 1 a 1 S ¯ b + S ¯ + e 2 a 2 I ¯ ( d 2 + h 3 ) , will be negative, provided that the following condition holds.
e 1 a 1 S ¯ b + S ¯ + e 2 a 2 I ¯ < ( d 2 + h 3 ) .  
Accordingly, the P.F.E.P. is a L.A.S. provided that Condition (20) is met.
The Jacobian matrix at m 3 = ( S ˜ , 0 , Y ˜ ) can be written as
J ( m 3 ) = [ α S ˜ + a 1 Y ˜ ( b + S ˜ ) 2 ( α + β ) S ˜ ( r μ ( 1 θ ) ( μ + Y ˜ ) 2 + a 1 b + S ˜ ) S ˜ 0 β S ˜ a 2 Y ˜ ( d 1 + h 2 ) 0 b e 1 a 1 ( b + S ˜ ) 2 Y ˜ e 2 a 2 Y ˜ 0 ] .  
The characteristic equation of J ( m 3 ) can be determined as follows:
( λ 2 T 3 λ + D 3 )   ( β S ˜ a 2 Y ˜ ( d 1 + h 2 ) λ ) = 0 ,  
where T 3 = α S ˜ + a 1 Y ˜ ( b + S ˜ ) 2 , and D 3 = b e 1 a 1 ( b + S ˜ ) 2 ( r μ ( 1 θ ) ( μ + Y ˜ ) 2 + a 1 b + S ˜ ) S ˜ Y ˜ > 0 .
It is easy to verify that all the eigenvalues have negative real parts and D.F.E.P. is L.A.S. if the following condition holds.
β S ˜ < a 2 Y ˜ + ( d 1 + h 2 ) .  
a 1 Y ˜ < α S ˜ ( b + S ˜ ) 2 .  
Finally, The Jacobian matrix at m 4 = ( S , I , Y ) can be written as
J ( m 4 ) = [ α + a 1 Y ( b + S ) 2 ( α + β ) S ( r μ ( 1 θ ) ( μ + Y ) 2 + a 1 b + S ) S β I 0 a 2 I b e 1 a 1 Y ( b + S ) 2 e 2 a 2 Y 0 ] = [ a i j ] 3 × 3 .  
Then, the characteristic equation of J ( m 4 ) is
λ 3 + A 1 λ 2 + A 2 λ + A 3 = 0 ,  
where
A 1 = a 11 , A 2 = ( a 23 a 32 + a 12 a 21 + a 13 a 31 ) , A 3 = a 23 ( a 11 a 32 a 12 a 31 ) a 21 a 32 a 13 ,
While
Δ = A 1   A 2 A 3 = a 11 a 12 a 21 + a 13 ( a 11 a 31 + a 21 a 32 ) + a 12 a 23 a 31 .
Consequently, the L.A.S requirements of P.E.P. can be constructed in the following theorem.
Theorem 3. 
A P.E.P. of the system (1) is L.A.S. if and only if the following condition is met.
max { a 1 Y ( b + S ) 2 + b e 1 a 1 ( α + β ) S e 2 a 2 ( b + S ) 2 , a 1 Y ( b + S ) 2 + e 2 a 2 β I ( b + S ) 2 b e 1 a 1   } < α .  
Proof. 
According to the Routh–Hurwitz criterion, the roots of the Jacobian matrix J ( m 4 ) have negative real parts provided that A 1 > 0 ,   A 3 > 0 ,   and   Δ > 0 . Direct computation shows that Condition (27) guarantees the satisfaction of Routh–Hurwitz criterion requirements. Hence the proof is done.
In the following, the persistence requirements of the system (1) are investigated. It is clear that system (1) has two subsystems lying in the non-negative quadrant of S I and S Y planes. These subsystems can be written as
S [ r α S ( α + β ) I h 1 ] = ϕ 1 ( S   ,   I ) I [ β S ( d 1 + h 2 ) ] = ϕ 2 ( S   ,   I )   .  
S [ r ( θ + μ ( 1 θ ) μ + Y ) α S a 1 Y b + S h 1 ] = ρ 1 ( S   ,   I ) Y [ e 1 a 1 S b + S ( d 2 + h 3 ) ] = ρ 2 ( S   ,   I )   .  
Clearly, each of these subsystems has three nonnegative equilibrium points that represent the projection of the corresponding boundary equilibrium points of the system (1).
Define the Dulac functions as B 1 ( S , I ) = 1 S I and B 2 ( S , Y ) = 1 S Y that are obviously B 1 ( S , I ) > 0 , B 2 ( S , Y ) > 0 , and C 1 functions in the I n t .   + 2 of the S I and S Y planes. Now, direct computation gives that
Δ ( S , I ) = ( B 1 ϕ 1 ) S + ( B 1 ϕ 2 ) I = α I < 0 .
Δ ( S , Y ) = ( B 2 ρ 1 ) S + ( B 2 ρ 21 ) Y = α Y + a 1 ( b + S ) 2 .
Then, the expression Δ ( S , I ) is always negative, and hence, the first subsystem (28) does not contain a periodic dynamic fall in the positive quadrant of the S I -plane. However, the expression Δ ( S , Y ) does not identically zero in the I n t .   + 2 of the S Y -plane and does not change the sign under the condition
a 1 Y ( b + S ) 2 < α .
Therefore, by using the Dulac–Bendixson criterion [41], there is no closed curve lying in the I n t .   + 2 of the S Y -plane for all the trajectories satisfying Condition (30). Hence, according to the Poincare–Bendixon theorem [41], the interior equilibrium points of the subsystems (28) and (29) that are given by m 2 and m 3 respectively, will be a globally asymptotically stable G.A.S whenever they are L.A.S. □
Theorem 4. 
Assume that Condition (30) holds, then system (1) is uniformly persistent if the following conditions are met.
α ( d 1 + h 2 ) < β ( r h 1 ) .
( d 2 + h 3 ) < e 1 a 1 S ¯ b + S ¯ + e 2 a 2 I ¯ .
  a 2 Y ˜ + ( d 1 + h 2 ) < β S ˜ .
Proof. 
Consider the following function: ( S , I , Y ) = S p 1   I p 2   Y p 3 , where p j , j = 1 , 2 , 3 are positive constants. Clearly, ξ ( S , I , Y ) > 0   for   all   ( S , I , Y ) I n t . + 3 and ξ ( S , I , Y ) 0 , when S 0 ,   I 0 , or Y 0 . Consequently, it is obtained that
Ω ( S , I , Y ) = ξ ( S , I , Y ) ξ ( S , I , Y ) = p 1 [ r ( θ + μ ( 1 θ ) μ + Y ) α ( S + I ) β I a 1 Y b + S h 1 ] + p 2 [ β S a 2 Y ( d 1 + h 2 ) ] + p 3 [ e 1 a 1 S b + S + e 2 a 2 I ( d 2 + h 3 ) ] .
Now the proof follows if Ω ( P ) > 0 for any boundary equilibrium point P , with a suitable choice of constants p 1 > 0 , p 2 > 0 , and p 3 > 0 .
Ω ( m 0 ) = p 1 ( r h 1 ) p 2 ( d 1 + h 2 ) p 3 ( d 2 + h 3 ) . Ω ( m 1 ) = p 2 ( β S ^ ( d 1 + h 2 ) ) + p 3 ( e 1 a 1 S ^ b + S ^ ( d 2 + h 3 ) ) . Ω ( m 2 ) = p 3 ( e 1 a 1 S ¯ b + S ¯ + e 2 a 2 I ¯ ( d 2 + h 3 ) ) . Ω ( m 3 ) = p 2 ( β S ˜ a 2 Y ˜ ( d 1 + h 2 ) ) .
Clearly, Condition (31) guarantees that Ω ( m 0 ) > 0 , under Condition (31) with a suitable choice of positive constants p j , so that p 1 is sufficiently larger than p 2 and p 3 . Moreover, Ω ( m 1 ) > 0 , for a suitable choice of positive constants p j , so that p 2 is sufficiently larger than p 3 . Finally, Ω ( m j ) > 0 ; j = 2 , 3 due to Conditions (32) and (33), respectively. Thus the proof is performed. □

4. Global Dynamics

In this section, an investigation of the global dynamics of the system (1) around each L.A.S. point is carried out utilizing a suitable Lyapunov function.
Theorem 5. 
The V.E.P. is a G.A.S. whenever it is L.A.S. provided that the following condition holds
e 2 < e 1 .
Proof. 
Define the function W 0 ( S , I , Y ) = S + I + 1 e 1 Y . Clearly, the function W 0 is a positive definite so that W 0 ( 0 , 0 , 0 ) = 0 and W 0 ( S , I , Y ) > 0 for all ( S , I , Y ) + 3 with ( S , I , Y ) ( 0 , 0 , 0 ) . Then, the derivative d W 0 d t can be determined as
d W 0 d t = r ( θ + μ ( 1 θ ) μ + Y ) S α S ( S + I ) h 1 S a 2 ( 1 e 2 e 1 ) I Y ( d 1 + h 2 ) I ( d 2 + h 3 ) e 1 Y .
Now, using Condition (34) and after some algebraic manipulation, it is obtained that
d W 0 d t ( r h 1 ) S ( d 1 + h 2 ) I ( d 2 + h 3 ) e 1 Y .
Therefore, d W 0 d t is negative definite due to the L.A.S. condition of m 0 . Hence, the V.E.P. is a G.A.S. point whenever it is L.A.S. under Condition (34). □
Theorem 6. 
The A.E.P. is a G.A.S. whenever it is a L.A.S. provided that the following conditions are met along with the Condition (34).
( α + β ) S ^ < ( d 1 + h 2 ) .
a 1 S ^ b +   r ( 1 θ ) S ^ μ < ( d 2 + h 3 ) e 1 .
Proof. 
Define the function W 1 ( S , I , Y ) = S ^ S u S ^ u d u + I + 1 e 1 Y . Clearly, the function W 1 is positive definite so that W 1 ( S ^ , 0 , 0 ) = 0 and W 1 ( S , I , Y ) > 0 for all ( S , I , Y ) + 3 with ( S , I , Y ) ( S ^ , 0 , 0 ) and S > 0 . Then, the derivative d W 1 d t can be determined as
d W 1 d t α ( S S ^ ) 2 [ ( d 1 + h 2 ) ( α + β ) S ^ ] I a 2 ( 1 e 2 e 1 ) I Y [ ( d 2 + h 3 ) e 1 a 1 S ^ b + S   r ( 1 θ ) S ^ μ + Y ] Y .
Now, using Conditions (34)–(36) and after some algebraic manipulation it is obtained that
d W 1 d t α ( S S ^ ) 2 [ ( d 1 + h 2 ) ( α + β ) S ^ ] I [ ( d 2 + h 3 ) e 1 a 1 S ^ b + S   r ( 1 θ ) S ^ μ + Y ] Y .
Therefore, d W 1 d t is negative definite due to the given conditions. Hence, the A.E.P. is a G.A.S., and the proof is complete. □
Theorem 7. 
The basin of attraction of the P.F.E.P is given by the region  B ( m 2 ) = { ( S , I , Y ) + 3 : S > S ¯ , I > I ¯ , Y 0 }  provided that the following condition is met along with Condition (34).
r ( 1 θ ) S ¯ μ + a 1 S ¯ b + a 2 I ¯ < ( d 2 + h 3 ) e 1   .
Proof. 
Consider the following function: W 2 ( S , I , Y ) = S ¯ S u S ¯ u d u + I ¯ I v I ¯ v d v + 1 e 1 Y . Clearly, the function W 2 ( S , I , Y ) > 0 is a continuously differentiable real-valued function for all ( S , I , Y ) + 3 with ( S , I , Y ) ( S ¯ , I ¯ , 0 ) and W 2 ( S ¯ , I ¯ , 0 ) = 0 . Then, the derivative d W 2 d t can be determined as
d W 2 d t α ( S S ¯ ) 2 α ( S S ¯ ) ( I I ¯ ) a 2 ( 1 e 2 e 1 ) I Y [ ( d 2 + h 3 ) e 1 r ( 1 θ ) S ¯ μ + Y a 1 S ¯ b + S a 2 I ¯ ] Y .
So, by using Conditions (34) and (37), it is then obtained after some algebraic manipulation that
d W 2 d t α ( S S ¯ ) 2 α ( S S ¯ ) ( I I ¯ ) [ ( d 2 + h 3 ) e 1 r ( 1 θ ) S ¯ μ + Y a 1 S ¯ b + S a 2 I ¯ ] Y .
Therefore, d W 2 d t is negative definite due to the given conditions in the region B ( m 2 ) . Hence, the proof is complete. □
Theorem 8. 
The basin of attraction of the D.F.E.P is given by the region  B ( m 3 ) = { ( S , I , Y ) + 3 : S > S ¯ ,   I 0 , Y > Y ˜ }  provided that the following conditions are met along with Condition (34).
a 1 Y ˜ ( b + S ) ( b + S ˜ ) < α
( α + β ) S ˜ < ( d 1 + h 2 ) + e 2 a 2 e 1 Y ˜
Proof. 
Consider the following function: W 3 ( S , I , Y ) = S ˜ S u S ˜ u d u + I + 1 e 1   Y ˜ Y w Y ˜ w d w . Clearly, the function W 3 ( S , I , Y ) > 0 is a continuously differentiable real-valued function for all ( S , I , Y ) + 3 with ( S , I , Y ) ( S ˜ , 0 , Y ˜ ) and W 3 ( S ˜ , 0 , Y ˜ ) = 0 . Then, the derivative d W 3 d t can be determined as
d W 3 d t [ α a 1 Y ˜ ( b + S ) ( b + S ˜ ) ] ( S S ˜ ) 2 [ μ r ( 1 θ ) ( μ + Y ) ( μ + Y ˜ ) + a 1 S ˜ ( b + S ) ( b + S ˜ ) ] ( S S ˜ ) ( Y Y ˜ ) a 2 ( 1 e 2 e 1 ) I Y [ ( d 1 + h 2 ) + e 2 a 2 e 1 Y ˜ ( α + β ) S ˜ ] I .
So, by using Conditions (34), (38), and (39), it is obtained after some algebraic manipulation that
d W 3 d t [ α a 1 Y ˜ ( b + S ) ( b + S ˜ ) ] ( S S ˜ ) 2 [ μ r ( 1 θ ) ( μ + Y ) ( μ + Y ˜ ) + a 1 S ˜ ( b + S ) ( b + S ˜ ) ] ( S S ˜ ) ( Y Y ˜ ) [ ( d 1 + h 2 ) + e 2 a 2 e 1 Y ˜ ( α + β ) S ˜ ] I .
Therefore, d W 3 d t is negative definite due to the given conditions in the region B ( m 3 ) . Hence, the proof is complete. □
Theorem 9. 
The basin of attraction of the P.E.P is given by the region  B ( m 4 ) = { ( S , I , Y ) + 3 : S > S ,   I I , Y > Y }  provided that the following conditions are met along with Condition (34). 
a 1 Y ( b + S ) ( b + S ) < α .
a 1 S ( b + S ) ( b + S ) < r   μ ( 1 θ ) ( μ + Y ) ( μ + Y ) .  
Proof. 
Consider the function: W 4 ( S , I , Y ) = S S u S u d u + I I v I v d v + 1 e 1 Y Y w Y w d w . Clearly, the function W 4 ( S , I , Y ) > 0 is a continuously differentiable real-valued function for all ( S , I , Y ) + 3 with ( S , I , Y ) ( S , I , Y ) and W 4 ( S , I , Y ) = 0 . Then, the derivative d W 4 d t can be determined as
d W 4 d t = ( α a 1 Y ( b + S ) ( b + S ) ) ( S S ) 2 α ( S S ) ( I I ) ( r   μ ( 1 θ ) ( μ + Y ) ( μ + Y ) a 1 S ( b + S ) ( b + S ) ) ( S S ) ( Y Y ) a 2 ( 1 e 2 e 1 ) ( I I ) ( Y Y ) .  
So, by using Conditions (34), (40) and (41), it is obtained that d W 4 d t is negative definite due to the given conditions in the region B ( m 4 ) . Hence, the proof is complete. □

5. Local Bifurcation

By using Sotomayor’s theorem [41], the local bifurcation (L.B.) around the equilibrium points of system (1) is explored. It is common knowledge that a non-hyperbolic equilibrium point is a requirement, but not a sufficient one, for bifurcation to take place. In order to ensure that the equilibrium point is non-hyperbolic at a given value of the parameter, the candidate bifurcation parameter is chosen. Now, rewrite system (1) in the form:
d X d t = F ( X )
where X = ( S , I , Y ) T and F = ( S f 1 , I f 2 , Y f 3 ) T with f i ; i = 1 , 2 , 3 represent the interaction functions on the right-hand side of the system (1). Then, using the Jacobian matrix of system (1), simple computation demonstrates that we have the following second directional derivative for every non-zero vector V = ( v 1 , v 2 , v 3 ) T , and ϑ + .
D 2 F ( X , ϑ ) ( V , V ) = [ k i 1 ] 3 × 1 ,
where
k 11 = 2 ( α a 1 b   Y ( b + S ) 3 ) v 1 2 2 ( r μ ( 1 θ ) ( μ + Y ) 2 + a 1 b   ( b + S ) 2 ) v 1 v 3 2 ( α + β ) v 1 v 2 + 2 ( r μ ( 1 θ ) S ( μ + Y ) 3 ) v 3 2 , k 21 = 2 β v 1 v 2 2 a 1 v 2 v 3 , k 31 = 2 b e 1 a 1   Y ( b + S ) 3 v 1 2 + 2 b e 1 a 1   ( b + S ) 2 v 1 v 3 + 2 e 2 a 2 v 2 v 3 .
While the third directional derivative is determined as
D 3 F ( X , ϑ ) ( V , V , V ) = [ h i 1 ] 3 × 1 ,
where
h 11 = 6 r ( 1 θ ) μ v 1 v 3 2 ( μ + Y ) 3 6 r ( 1 θ ) μ S v 3 3 ( μ + Y ) 4 6 b a 1 Y v 1 3 ( b + S ) 4 + 6 b a 1 v 1 2 v 3 ( b + S ) 3 , h 21 = 0 , h 31 = 6 b a 1 e 1 Y v 1 3 ( b + S ) 4 6 b a 1 e 1 v 1 2 v 3 ( b + S ) 3 .
Theorem 10. 
A transcritical bifurcation takes place near the V.E.P. when the parameter  r  moves through the point  r = h 1 ( r ) .
Proof. 
As r = r the Jacobian matrix at V.E.P. becomes
J 0 = J ( m 0 , r ) = [ 0 0 0 0 ( d 1 + h 2 ) 0 0 0 ( d 2 + h 3 ) ] .
Hence, the eigenvalues of J 0 are given in Equation (14) with λ 01 = 0 . Accordingly, the eigenvectors corresponding to λ 01 = 0 of J 0 and J 0 T can be determined as follows, respectively.
V 0 = ( v 10 , v 20 , v 30 ) T = ( 1 , 0 , 0 ) T ;   U 0 = ( u 10 , u 20 , u 30 ) T = ( 1 , 0 , 0 ) T .
Moreover, direct computation shows that
F r ( X ,   r ) = ( ( θ + μ ( 1 θ ) μ + Y ) S 0 0 ) F r ( m 0 ,   r ) = ( 0 0 0 ) U 0 T F r ( m 0 ,   r ) = 0 . D F r ( X ,   r ) . V 0 = ( θ + μ ( 1 θ ) μ + Y 0 0 ) D F r ( m 0 ,   r ) . V 0 = ( 1 0 0 ) U 0 T [ D F r ( m 0 ,   r ) . V 0 ] = 1 .
Using Equation (43) gives that
D 2 F ( m 0 ,   r ) ( V 0 , V 0 ) = ( 2 α 0 0 ) U 0 T [ D 2 F ( m 0 ,   r ) ( V 0 , V 0 ) ] = 2 α .
Hence, according to Sotomayor’s theorem, the transcritical bifurcation takes place, and the proof is complete. □
Theorem 11. 
Assume that the second part of Condition (17) is satisfied, then a transcritical bifurcation takes place near the A.E.P when the parameter  β  moves through the point  β = α ( d 1 + h 2 ) r h 1 ( β ) .
Proof. 
As β = β the Jacobian matrix at A.E.P. becomes
J 1 = J ( m 1 , β ) = [ r + h 1 ( α + β ) ( r h 1 ) α ( r ( 1 θ ) ( b α + r h 1 ) + a 1 α μ ) ( r h 1 ) μ α ( b α + r h 1 ) 0 0 0 0 0 a 1 e 1 ( r h 1 ) ( b α + r h 1 ) ( d 2 + h 3 ) ] .
Hence, the eigenvalues of J 1 are given in Equation (16) with λ 12 = 0 . Accordingly, under the second part of Condition (17), the eigenvectors corresponding to λ 12 = 0 of J 1 and J 1 T can be determined as follows, respectively.
V 1 = ( v 11 , v 21 , v 31 ) T = ( ( α + β ) α , 1 , 0 ) T ;   U 1 = ( u 11 , u 21 , u 31 ) T = ( 0 , 1 , 0 ) T .
Moreover, direct computation shows that
F β ( m 1 , β ) = ( 0 0 0 ) U 1 T F β ( m 1 , β ) = 0 .
D F β ( m 1 , β ) . V 1 = ( S ^ S ^ 0 ) U 1 T [ D F β ( m 1 , β ) . V 1 ] = S ^ .
Using Equation (43) gives that
D 2 F ( m 1 , β ) ( V 1 , V 1 ) = ( 0 2 β ( α + β α ) 0 )          U 1 T [ D 2 F ( m 1 , β ) ( V 1 , V 1 ) ] = 2 β ( α + β α ) .
Hence, according to Sotomayor’s theorem, the transcritical bifurcation takes place, and the proof is complete. □
Theorem 12. 
A transcritical bifurcation takes place near the P.F.E.P when the parameter  h 3  moves through the point  h 3 = e 1 a 1 S ¯ b + S ¯ + e 2 a 2 I ¯ d 2 ( h 3 )  provided that the following condition is met
2 b e 1 a 1   ( b + S ¯ ) 2 ρ 1 + 2 e 2 a 2 ρ 2 0 .
Otherwise, a pitchfork bifurcation takes place.
Proof. 
As h 3 = h 3 the Jacobian matrix at P.F.E.P. becomes
J 2 = J ( m 2 , h 3 ) = [ α S ¯ ( α + β ) S ¯ ( r ( 1 θ ) μ + a 1 b + S ¯ ) S ¯ β I ¯ 0 a 2 I ¯ 0 0 0 ] = [ b i j ] .
Hence, the eigenvalues of J 2 are given in Equation (19), with λ 23 = 0 . Accordingly, the eigenvectors corresponding to λ 23 = 0 of J 2 and J 2 T can be determined as follows, respectively.
V 2 = ( v 12 , v 22 , v 32 ) T = ( ρ 1 , ρ 2 , 1 ) T ;   U 2 = ( u 12 , u 22 , u 32 ) T = ( 0 , 0 , 1 ) T ,
where ρ 1 = b 23 b 21 > 0 and ρ 2 = b 11 b 23 b 13 b 21 b 12 b 21 < 0 .
Moreover, direct computation shows that
F h 3 ( X ,   h 3 ) = ( 0 0 Y ) F h 3 ( m 2 , h 3 ) = ( 0 0 0 ) U 2 T F h 3 ( m 2 , h 3 ) = 0 . D F h 3 ( X ,   h 3 ) . V 2 = ( 0 0 1 ) D F h 3 ( m 2 , h 3 ) . V 2 = ( 0 0 1 ) U 2 T [ D F h 3 ( m 2 , h 3 ) . V 2 ] = 1 .
Using Equation (43) gives that
D 2 F ( m 2 , h 3 ) . ( V 2 , V 2 ) = [ k i 1 ( m 2 , h 3 ) ] 3 × 1 ,
where
k 11 ( m 2 , h 3 ) = 2 α ρ 1 2 2 ( r ( 1 θ ) μ + a 1 b   ( b + S ¯ ) 2 ) ρ 1 2 ( α + β ) ρ 1 ρ 2 + 2 ( r ( 1 θ ) S ¯ μ 2 ) , k 21 ( m 2 , h 3 ) = 2 β ρ 1 ρ 2 2 a 1 ρ 2 , k 31 ( m 2 , h 3 ) = 2 b e 1 a 1   ( b + S ¯ ) 2 ρ 1 + 2 e 2 a 2 ρ 2 .
Thus, it is obtained that
U 2 T [ D 2 F ( m 2 , h 3 ) . ( V 2 , V 2 ) ] = 2 b e 1 a 1   ( b + S ¯ ) 2 ρ 1 + 2 e 2 a 2 ρ 2 .
Clearly, under Condition (45), it is found that U 2 T [ D 2 F ( m 2 , h 3 ) . ( V 2 , V 2 ) ] 0 . Hence, according to Sotomayor’s theorem, transcritical bifurcation takes place. However, if Condition (45) is violated, then Equation (44) yields that
D 3 F ( m 2 , h 3 ) . ( V 2 , V 2 , V 2 ) = ( 6 r ( 1 θ ) ρ 1 μ 2 6 r ( 1 θ ) S ¯ μ 3 + 6 b a 1 ρ 1 2 ( b + S ¯ ) 3 0 6 b a 1 e 1 ρ 1 2 ( b + S ¯ ) 3 ) .
Consequently, U 2 T [ D 3 F ( m 2 , h 3 ) . ( V 2 , V 2 , V 2 ) ] = 6 b a 1 e 1 ρ 1 2 ( b + S ¯ ) 3 , Hence, a pitchfork bifurcation takes place, and the proof is complete. □
Theorem 13. 
Assume that Condition (24) holds. Then, a transcritical bifurcation takes place near the D.F.E.P when the parameter  h 2  moves through the point  h 2 = β S ˜ + a 2 Y ˜ d 1 ( h 2 )  provided that the following condition is met
2 β ρ 3 2 a 1 ρ 4 0 .
Otherwise, a pitchfork bifurcation does not occur.
Proof. 
As h 3 = h 3 the Jacobian matrix at D.F.E.P. becomes
J 3 = J ( m 3 , h 2 ) = [ α S ˜ + a 1 Y ˜ ( b + S ˜ ) 2 ( α + β ) S ˜ ( r μ ( 1 θ ) ( μ + Y ˜ ) 2 + a 1 b + S ˜ ) S ˜ 0 0 0 b e 1 a 1 ( b + S ˜ ) 2 Y ˜ e 2 a 2 Y ˜ 0 ] = [ c i j ] .
Hence, the eigenvalues of J 3 are given in Equation (22) with λ 32 = 0 . Accordingly, due to Condition (24), the eigenvectors corresponding to λ 32 = 0 of J 3 and J 3 T can be determined as follows, respectively.
V 3 = ( v 13 , v 23 , v 33 ) T = ( ρ 3 , 1 , ρ 4 ) T ;   U 3 = ( u 13 , u 23 , u 33 ) T = ( 0 , 1 , 0 ) T ,
where ρ 3 = c 32 c 31 < 0 and ρ 4 = c 11 c 32 c 12 c 31 c 13 c 31 .
Moreover, direct computation shows that
F h 2 ( X ,   h 2 ) = ( 0 I 0 ) F h 2 ( m 3 , h 2 ) = ( 0 0 0 ) U 3 T F h 2 ( m 3 , h 2 ) = 0 . D F h 2 ( X ,   h 2 ) . V 3 = ( 0 1 0 ) D F h 2 ( m 3 , h 2 ) . V 3 = ( 0 1 0 ) U 3 T [ D F h 2 ( m 3 , h 2 ) . V 3 ] = 1 .
Using Equation (43) gives that:
D 2 F h 2 ( m 3 , h 2 ) . ( V 3 , V 3 ) = [ k i 1 ( m 3 , h 2 ) ] 3 × 1 ,
where,
k 11 ( m 3 , h 2 * ) = 2 ( α a 1 b   Y ˜ ( b + S ˜ ) 3 ) ρ 3 2 2 ( r μ ( 1 θ ) ( μ + Y ˜ ) 2 + a 1 b   ( b + S ˜ ) 2 ) ρ 3 ρ 4 2 ( α + β ) ρ 3 + 2 ( r μ ( 1 θ ) S ˜ ( μ + Y ˜ ) 3 ) ρ 4 2 k 21 ( m 3 , h 2 ) = 2 β ρ 3 2 a 1 ρ 4 , k 31 ( m 3 , h 2 ) = 2 b e 1 a 1   Y ˜ ( b + S ˜ ) 3 ρ 3 2 + 2 b e 1 a 1   ( b + S ˜ ) 2 ρ 3 ρ 4 + 2 e 2 a 2 ρ 4 .
Thus, it is obtained that
U 3 T [ D 2 F h 2 ( m 3 , h 2 ) . ( V 3 , V 3 ) ] = 2 β ρ 3 2 a 1 ρ 4 .
Clearly, under Condition (46), it is found that U 3 T [ D 2 F ( m 3 , h 2 ) . ( V 3 , V 3 ) ] 0 . Hence, according to Sotomayor’s theorem, transcritical bifurcation takes place. However, if Condition (46) is violated, then Equation (44) yields that:
D 3 F ( m 3 , h 2 ) . ( V 3 , V 3 , V 3 ) = ( 6 r ( 1 θ ) μ ρ 3 ρ 4 2 ( μ + Y ˜ ) 3 6 r ( 1 θ ) μ S ˜ ρ 4 3 ( μ + Y ˜ ) 4 6 b a 1 Y ˜ ρ 3 3 ( b + S ˜ ) 4 + 6 b a 1 ρ 3 2 ρ 4 ( b + S ˜ ) 3 0 6 b a 1 e 1 Y ˜ ρ 3 3 ( b + S ˜ ) 4 6 b a 1 e 1 ρ 3 2 ρ 4 ( b + S ˜ ) 3 ) .
Consequently, U 3 T [ D 3 F ( m 3 , h 2 ) . ( V 3 , V 3 , V 3 ) ] = 0 , Hence a pitchfork bifurcation does not occur, and the proof is complete. □
Theorem 14. 
A saddle-node bifurcation takes place near the P.E.P. when the parameter  e 2  moves through the point  e 2 = a 12 a 23 a 31 a 2 Y ( a 11 a 23 a 21 a 13 ) ( e 2 ) provided that the following conditions are met
a 1 Y ( b + S ) 2 < α .
ρ 7 k 11 ( m 4 , e 2 ) + ρ 8 k 21 ( m 4 , e 2 ) + k 31 ( m 4 , e 2 ) 0 .
Proof. 
As e 2 = e 2 the Jacobian matrix at P.E.P. becomes
J 4 = J ( m 4 , e 2 ) = [ α + a 1 Y ( b + S ) 2 ( α + β ) S ( r μ ( 1 θ ) ( μ + Y ) 2 + a 1 b + S ) S β I 0 a 2 I b e 1 a 1 Y ( b + S ) 2 e 2 a 2 Y 0 ] = [ a i j ] .
Direct computation shows that the determinant of J 4 that A 3 gives in the characteristic Equation (26) at e 2 = e 2 equals zero. Therefore, the characteristic Equation (26) has zero eigenvalue, say λ 3 = 0 , along with the other two negative real parts eigenvalues due to Condition (47), which are given by λ 1 , 2 = A 1 2 ± A 1 2 4 A 2 2 .
Accordingly, the eigenvectors corresponding to λ 3 = 0 of J 4 and J 4 T can be determined as follows, respectively.
V 4 = ( v 14 , v 24 , v 34 ) T = ( ρ 5 , ρ 6 , 1 ) T ;   U 4 = ( u 14 , u 24 , u 34 ) T = ( ρ 7 , ρ 8 , 1 ) T ,
where ρ 5 = a 23 a 21 > 0 , ρ 6 = a 11 a 23 a 13 a 21 a 12 a 21 < 0 , ρ 7 = a 32 a 12 > 0 and ρ 8 = a 11 a 32 a 12 a 31 a 12 a 21 .
Moreover, direct computation shows that
F e 2 ( X ,   e 2 ) = ( 0 0 a 2 I Y ) F e 2 ( m 4 , e 2 ) = ( 0 0 a 2 I Y ) U 4 T F e 2 ( m 4 , e 2 ) = a 2 I Y .
Using Equation (43) gives that
D 2 F ( m 4 , e 2 ) ( V 4 , V 4 ) = [ k i 1 ( m 4 , e 2 ) ] 3 × 1
where
k 11 ( m 4 , e 2 ) = 2 ( α a 1 b   Y ( b + S ) 3 ) ρ 5 2 2 ( r μ ( 1 θ ) ( μ + Y ) 2 + a 1 b   ( b + S ) 2 ) ρ 5 2 ( α + β ) ρ 5 ρ 6 + 2 ( r μ ( 1 θ ) S ( μ + Y ) 3 ) , k 21 ( m 4 , e 2 ) = 2 β ρ 5 ρ 6 2 a 1 ρ 6 , k 31 ( m 4 , e 2 ) = 2 b e 1 a 1   Y ( b + S ) 3 ρ 5 2 + 2 b e 1 a 1   ( b + S ) 2 ρ 5 + 2 e 2 a 2 ρ 6 .
Therefore, it is obtained that
U 4 T [ D 2 F ( m 4 , e 2 ) ( V 4 , V 4 ) ] = ρ 7 k 11 ( m 4 , e 2 ) + ρ 8 k 21 ( m 4 , e 2 ) + k 31 ( m 4 , e 2 ) .
Clearly, under Condition (47), it is found that U 4 T [ D 2 F ( m 4 , e 2 ) ( V 4 , V 4 ) ] 0 . Hence, according to Sotomayor’s theorem, saddle-node bifurcation takes place. This completes the proof. □

6. Numerical Simulation

The dynamics of the system (1) are numerically examined in the following. Confirming the analytical results and identifying the control set of parameters on the dynamics of the system are the goals. With the numerical solution, the following hypothetical parameter values are employed.
r = 1 ,   θ = 0.5 ,   μ = 0.4 ,   α = 0.2 , β = 0.3 ,   a 1 = 0.5 ,   b = 0.2 ,   a 2 = 0.6 , d 1 = 0.1 ,   d 2 = 0.1 ,   e 1 = 0.4 ,   e 2 = 0.4 ,   h 1 = 0.2 ,   h 2 = 0.2 ,   h 3 = 0.2 .  
It is observed that for the above set of data, the system (1) approaches asymptotically to the P.E.P. that is determined by m 4 = ( 1.45 , 0.51 , 0.22 ) as presented in Figure 1 below.
It is obvious from Figure 1 that system (1) has a unique asymptotically stable P.E.P. Now, the influence of varying the value of r on the dynamic of the system (1) is studied, and the obtained results are presented at selected values in Figure 2. It is obtained that, for r belongs to ( 0 , 0.2 ) , [ 0.2 , 0.4 ] , ( 0.4 , 0.68 ) , [ 0.68 , 2.06 ) , [ 2.06 , 6.61 ] , and r 6.62 , the system’s (1) solution approaches m 0 , m 1 , m 2 , m 4 , bistable between m 4 and 3D periodic attractor, and 3d periodic attractor, respectively.
According to Figure 2, there are five bifurcation points belonging to the range of r , including the transcritical point. For the α to belong to ( 0 , 0.13 ] , ( 0.13 , 0.42 ) , [ 0.42 , 0.8 ) , and α 0.8 , the solution of system (1) converges 3D periodic attractor, m 4 , m 2 , and m 1 , see Figure 3 for the selected values of α .
According to Figure 3, there are three bifurcation points in the range of α . Similar behaviors are obtained when varying the parameter h 2 as those obtained with α .
Now, the influence of varying the parameters θ , μ and b in the case of the system (1) having periodic dynamics is investigated in Figure 4, Figure 5 and Figure 6. It is obtained that for the data set (49) with α = 0.1 and θ belonging to [ 0 , 0.9 ) , and [ 0.9 , 1 ] the solution converges 3D periodic attractor, and m 4 respectively. While for the data set (49) with α = 0.1 and μ belonging to [ 0 , 11 ) , and μ 11 , the solution converges 3D periodic attractor, and m 4 respectively, see Figure 4 for the selected values of θ , and μ . On the other hand, for the data set (49) with α = 0.1 and varying b as shown in Figure 5 and Figure 6.
According to Figure 4, Figure 5 and Figure 6, the parameters θ , μ and b control the stability of the system (1). It is observed that the system converges to m 1 , m 2 , m 4 , and m 2 when the parameter β belongs to ( 0 , 0.07 ] , ( 0.07 , 0.1 ) , [ 0.1 , 0.82 ] , and β > 0.82 , respectively; see Figure 7 for selected values of β . However, the system converges, m 2 , m 4 , m 3 , a 3D periodic attractor, and a 2D periodic attractor in the S Y -plane when the parameter a 1 belongs to ( 0 , 0.03 ] , ( 0.03 , 0.83 ) , [ 0.83 , 0.86 ] , [ 0.88 , 1.95 ) , and a 1 1.95 , respectively, as shown in Figure 8 at the selected values. On the other hand, for a 2 in the ranges ( 0 , 0.28 ] and a 2 > 0.28 it is noted that system (1) converges m 2 and m 4 respectively as presented in Figure 9 at typical values of a 2 . However, for d 1 in the ranges ( 0 , 0.64 ) and d 1 0.64 , it is noted that system (1) converges m 4 and m 2 respectively, as shown in Figure 10 at the typical value of d 1 . Finally, it is obtained that varying the parameters d 1 and d 2 have a similar influence on the dynamics of the system (1).
Now, it is found that the system converges to m 4 , m 3 , a 2D periodic attractor in the S Y -plane, and a 3D periodic attractor when the parameter e 1 belongs to [ 0 , 0.66 ) , [ 0.66 , 0.7 ] , ( 0.7 , 0.74 ] , and ( 0.74 , 1 ] respectively, see Figure 11 for selected values of e 1 . However, the system converges, m 2 , m 4 , and bistable between the 3D periodic attractor and m 4 when the parameter e 2 belongs to [ 0 , 0.18 ] , ( 0.18 , 0.89 ) , and [ 0.89 , 1 ] , respectively, as shown in Figure 12 using selected values of e 2 . On the other hand, for h 1 in the ranges ( 0 , 0.53 ) , [ 0.53 , 0.8 ) , ( 0.8 , 1 ) and h 1 1 , it is noted that system (1) converges m 4 , m 2 , m 1 , and m 0 , respectively, as presented in Figure 13 at typical values of h 1 . However, for h 3 in the ranges ( 0 , 0.07 ) , [ 0.07 , 0.08 ] , ( 0.08 , 0.35 ] and h 3 > 0.35 it is noted that system (1) converges chaotic attractor, m 3 , m 4 , and m 2 , respectively, as shown in Figure 14 at a typical value of h 3 . Finally, the bifurcation diagram is given in Figure 15 shows clearly the existence of chaos in the range h 3 ( 0 , 0.07 ) .

7. Conclusions

In this study, the mathematical formulation and investigation of the effect of fear on a harvested prey–predator system in the presence of the disease among the prey population are carried out. Due to the difference in prey speed between healthy and diseased prey, two different types of Holling’s functional responses are used. The proposed model’s solution characteristics are investigated. The equilibrium points that are biologically viable have all been identified. Each equilibrium point’s local and global stability requirements are identified. The system’s persistence requirements are computed. Near the equilibrium points, all conceivable local bifurcation types are investigated. It is observed that system (1) is seen to undergo a transcritical bifurcation near the A.E.P. However, near the P.F.E.P. and D.F.E.P., the system has a transcritical bifurcation if a certain condition is met; otherwise, the system has a pitchfork bifurcation. Finally, saddle-node bifurcation occurs near the P.E.P. if certain criteria are met.
Finally, using a hypothetical set of parameter values, numerical research is conducted to confirm the analytical evaluation results and specify the impact of changing the parameter values. Regarding the impact of changing the parameter values on the system’s dynamic, the following conclusions are reached numerically (1). System (1) exhibits a variety of behaviors, including stable point, periodic, bistable, and chaotic behavior, because it is extremely sensitive to changes in its parameter values. From the perspective of how the parameters affect the system’s dynamic, the parameters can be divided into three groups: stabilizing, destabilizing, and extinction. The system (1) is stabilized at the positive equilibrium point in the stabilizing group by raising the parameters (the degree of fear and the minimal cost of fear, the predator’s half-saturation constant, and the predator’s maximum attack rate on the infected prey). Increasing the parameters, however, in the destabilizing group creates a destabilizing system (1), and the solution approaches asymptotically to 3D periodic attractors (the susceptible prey’s intrinsic growth rate and the conversion rates of the susceptible and infected biomass to predator biomass). Finally, every other parameter in the system (1) belongs to the extinction group because, when increased, it results in the extinction of one or more species.

Author Contributions

The authors H.A.I. and R.K.N. contributed to the conception, design model, acquisition of 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.

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to thank the editor and referees for their valuable comments and suggestions, which improved the quality of this manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Hsu, S.-B.; Hwang, T.-W.; Kuang, Y. Global analysis of the Michaelis–Menten type ratio-dependent predator–prey system. J. Math. Biol. 2001, 42, 489–506. [Google Scholar] [CrossRef] [PubMed]
  2. Naji, R.K.; Balasim, A.T. On the dynamical behavior of three species food web model. Chaos Solitons Fractals 2007, 34, 1636–1648. [Google Scholar] [CrossRef]
  3. Liu, X.; Lou, Y. Global dynamics of a predator–prey model. J. Math. Anal. Appl. 2010, 371, 323–340. [Google Scholar] [CrossRef] [Green Version]
  4. Naji, R.K.; Upadhyay, R.K.; Rai, V. Dynamical consequences of predator interference in a tri-trophic model food chain. Nonlinear Anal. Real World Appl. 2010, 11, 809–818. [Google Scholar] [CrossRef]
  5. Gonzalez-Olivares, E. Multiple limit cycles in a Gause type predator-prey model with Holling type III functional response and Allee effect on prey. Bull. Math. Biol. 2011, 73, 1378–1397. [Google Scholar] [CrossRef]
  6. Xiaolei, Y.; Fuqin, S. Global dynamics of a predator-prey model incorporating a constant prey refuge. Electron. J. Differ. Equ. 2013, 2013, 1–8. [Google Scholar]
  7. Llibre, J.; Xiao, D. Global Dynamics of a Lotka-Volterra Model with Two Predators Competing for One Prey. SIAM J. Appl. Math. 2014, 74, 434–453. [Google Scholar] [CrossRef] [Green Version]
  8. Xu, C.; Yuan, S.; Zhang, T. Global dynamics of a predator–prey model with defense mechanism for prey. Appl. Math. Lett. 2016, 62, 42–48. [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. Volterra, V. Fluctuations in the abundance of a species considered mathematically. Nature 1926, 118, 558–560. [Google Scholar] [CrossRef] [Green Version]
  11. Lotka, A.J. Elements of physical biology. Nature 1926, 21, 341–343. [Google Scholar] [CrossRef]
  12. Holling, S.C. Some characteristics of simple types of predation and parasitism. Can. Entomol. 1959, 91, 385–398. [Google Scholar] [CrossRef]
  13. Anderson, R.M.; May, R. The invasion, persistence, and spread of infectious diseases within animal and plant communities. Philos. Trans. R. Soc. B Biol. Sci. 1986, 314, 533–570. [Google Scholar] [CrossRef]
  14. Pal, A.K.; Samanta, G.P. Stability analysis of an eco-epidemiological model incorporating a prey refuge. Nonlinear Anal. Model. Control 2010, 15, 473–491. [Google Scholar] [CrossRef]
  15. Naji, R.K.; Mustafa, A.N. The Dynamics of an Eco-Epidemiological Model with Nonlinear Incidence Rate. J. Appl. Math. 2012, 2012, 852631. [Google Scholar] [CrossRef]
  16. Sahoo, B.; Poria, S. Diseased prey-predator model with general Holling type interactions. Appl. Math. Comput. 2014, 226, 83–100. [Google Scholar] [CrossRef]
  17. Jana, S.; Kar, T.K. Modeling and analysis of a prey-predator system with disease in the prey. Chaos Solitons Fractals 2013, 47, 42–53. [Google Scholar] [CrossRef]
  18. Bera, S.P.; Maiti, A.; Samanta, G.P. A Prey-predator Model with Infection in Both Prey and Predator. Filomat 2015, 29, 1753–1767. [Google Scholar] [CrossRef]
  19. 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, 17. [Google Scholar] [CrossRef]
  20. Hugo, A.; Simanjilo, E. Analysis of an Eco-Epidemiological Model under Optimal Control Measures for Infected Prey. Appl. Appl. Math. Int. J. 2019, 14, 117–138. [Google Scholar]
  21. Ibrahim, H.A.; Naji, R.K. A Prey-Predator Model with Michael-Mentence Type of Predator Harvesting and Infectious Disease in Prey. Iraqi J. Sci. 2020, 61, 1146–1163. [Google Scholar] [CrossRef]
  22. Hussein, W.; Satar, H.A. The Dynamics of a Prey-Predator Model with Infectious Disease in Prey: Role of Media Coverage. Iraqi J. Sci. 2021, 62, 4930–4952. [Google Scholar] [CrossRef]
  23. Bezabih, A.F.; Edessa, G.K.; Rao, K.P. Ecoepidemiological Model and Analysis of Prey-Predator System. J. Appl. Math. 2021, 2021, 6679686. [Google Scholar] [CrossRef]
  24. Savadogo, A.; Sangaré, B.; Ouedraogo, H. A mathematical analysis of prey-predator population dynamics in the presence of an SIS infectious disease. Res. Math. 2022, 9, 2020399. [Google Scholar] [CrossRef]
  25. Ganguli, C.; Kar, T.K. Optimal harvesting of a prey-predator model with variable carrying capacity. Int. J. Biomath. 2017, 10, 1750069. [Google Scholar] [CrossRef]
  26. Raymond, C.; Hugo, A.; Kung’aro, M. Modeling Dynamics of Prey-Predator Fishery Model with Harvesting: A Bioeconomic Model. J. Appl. Math. 2019, 2019, 2601648. [Google Scholar] [CrossRef]
  27. Satar, H.A.; Naji, R.K. Stability and Bifurcation in a Prey–Predator–Scavenger System with Michaelis–Menten Type of Harvesting Function. Differ. Equ. Dyn. Syst. 2022, 30, 933–956. [Google Scholar] [CrossRef]
  28. Liu, C.; Zhang, Q.; Zhang, Y.; Duan, X. Bifurcation and control in a differential-algebraic harvested prey-predator model with stage structure for predator. Int. J. Bifurc. Chaos 2008, 18, 3159–3168. [Google Scholar] [CrossRef]
  29. Taylor, R. Predation; Chapman & Hall: New York, NY, USA, 1984. [Google Scholar]
  30. Lima, S.; Dill, L.M. Behavioral decisions made under the risk of predation: A review and prospectus. Can. J. Zool. 1990, 68, 619–640. [Google Scholar] [CrossRef]
  31. Wang, X.; Zanette, L.; Zou, X. Modelling the fear effect in predator-prey interactions. J. Math. Biol. 2016, 75, 1179–1204. [Google Scholar] [CrossRef]
  32. Wang, X.; Zou, X. Modeling the fear effect in predator-prey interactions with adaptive avoidance of predators. Bull. Math. Biol. 2017, 79, 13–25. [Google Scholar] [CrossRef] [PubMed]
  33. Das, A.; Samanta, G.P. Modeling the fear effect on a stochastic prey-predator system with additional food for the predator. J. Phys. A Math. Theor. 2018, 51, 465601. [Google Scholar] [CrossRef]
  34. Wang, J.; Cai, Y.; Fu, S.; Wang, W. The effect of the fear factor on the dynamics of a predator-prey model incorporating the prey refuge. Chaos 2019, 29, 83109. [Google Scholar] [CrossRef] [PubMed]
  35. Zhanga, H.; Cai, Y.; Fu, S.; Wan, W. Impact of the fear effect in a prey-predator model incorporating a prey refuge. Appl. Math. Comput. 2019, 356, 328–337. [Google Scholar] [CrossRef]
  36. 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]
  37. 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]
  38. Abbas, Z.S.; Naji, R.K. Modeling and Analysis of the Influence of Fear on a Harvested Food Web System. Mathematics 2022, 10, 3300. [Google Scholar] [CrossRef]
  39. 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]
  40. Hale, J.K. Theory of Functional Differential Equations; Springer: Heidelberg, Germany, 1977. [Google Scholar]
  41. Perko, L. Differential Equations and Dynamical Systems, 3rd ed.; Springer: New York, NY, USA, 2001. [Google Scholar]
Figure 1. The trajectories of system (1) using data set (49) and starting from different initial points approach asymptotically m 4 = ( 1.45 , 0.51 , 0.22 ) : (a) 3D Phase plot. (b) Time series of solutions.
Figure 1. The trajectories of system (1) using data set (49) and starting from different initial points approach asymptotically m 4 = ( 1.45 , 0.51 , 0.22 ) : (a) 3D Phase plot. (b) Time series of solutions.
Mathematics 11 02909 g001
Figure 2. The trajectories of system (1) using data set (49) with different values of r and starting from various initial points. (a) Three-dimensional Phase plot converges m 0 . (b) Three-dimensional Phase plot converges m 1 = ( 0.75 , 0 , 0 ) . (c) Three-dimensional Phase plot converges m 2 = ( 1 , 0.4 , 0 ) . (d) Three-dimensional Phase plot converges m 4 = ( 3.3 , 0.46 , 1.15 ) . (e) Three-dimensional Phase plot for bistable between 3D periodic and m 4 = ( 4.38 , 0.45 , 1.69 ) . (f) Three-dimensional Phase plot converges 3D periodic attractor.
Figure 2. The trajectories of system (1) using data set (49) with different values of r and starting from various initial points. (a) Three-dimensional Phase plot converges m 0 . (b) Three-dimensional Phase plot converges m 1 = ( 0.75 , 0 , 0 ) . (c) Three-dimensional Phase plot converges m 2 = ( 1 , 0.4 , 0 ) . (d) Three-dimensional Phase plot converges m 4 = ( 3.3 , 0.46 , 1.15 ) . (e) Three-dimensional Phase plot for bistable between 3D periodic and m 4 = ( 4.38 , 0.45 , 1.69 ) . (f) Three-dimensional Phase plot converges 3D periodic attractor.
Mathematics 11 02909 g002aMathematics 11 02909 g002b
Figure 3. The trajectories of system (1) using data set (49) with different values of α and starting from various initial points. (a) Three-dimensional Phase plot converges 3D periodic attractor. (b) Three-dimensional Phase plot converges m 4 = ( 1.17 , 0.53 , 0.09 ) . (c) Three-dimensional Phase plot converges m 2 = ( 1 , 0.37 , 0 ) . (d) Three-dimensional Phase plot converges m 1 = ( 0.94 , 0 , 0 ) .
Figure 3. The trajectories of system (1) using data set (49) with different values of α and starting from various initial points. (a) Three-dimensional Phase plot converges 3D periodic attractor. (b) Three-dimensional Phase plot converges m 4 = ( 1.17 , 0.53 , 0.09 ) . (c) Three-dimensional Phase plot converges m 2 = ( 1 , 0.37 , 0 ) . (d) Three-dimensional Phase plot converges m 1 = ( 0.94 , 0 , 0 ) .
Mathematics 11 02909 g003
Figure 4. The trajectories of system (1) using data set (49) with α = 0.1 and selected values of θ , and μ . (a) Three-dimensional Phase plot converges 3D periodic attractor. (b) Three-dimensional Phase plot converges m 4 = ( 3.67 , 0.45 , 1.33 ) . (c) Three-dimensional Phase plot converges 3D periodic attractor. (d) Three-dimensional Phase plot converges m 4 = ( 3.83 , 0.45 , 1.42 ) .
Figure 4. The trajectories of system (1) using data set (49) with α = 0.1 and selected values of θ , and μ . (a) Three-dimensional Phase plot converges 3D periodic attractor. (b) Three-dimensional Phase plot converges m 4 = ( 3.67 , 0.45 , 1.33 ) . (c) Three-dimensional Phase plot converges 3D periodic attractor. (d) Three-dimensional Phase plot converges m 4 = ( 3.83 , 0.45 , 1.42 ) .
Mathematics 11 02909 g004
Figure 5. The trajectories of system (1) using data set (49) with α = 0.1 and selected values of b . (a) Three-dimensional Phase plot converges chaotic attractor. (b) Three-dimensional Phase plot converges m 4 = ( 1.85 , 0.71 , 0.42 ) .
Figure 5. The trajectories of system (1) using data set (49) with α = 0.1 and selected values of b . (a) Three-dimensional Phase plot converges chaotic attractor. (b) Three-dimensional Phase plot converges m 4 = ( 1.85 , 0.71 , 0.42 ) .
Mathematics 11 02909 g005
Figure 6. The bifurcation diagram as a function of b ( 0 , 0.6 ) with the data set (49) and α = 0.1 . (a) The trajectory of S as a function of b . (b) The trajectory of I as a function of b . (c) The trajectory of Y as a function of b .
Figure 6. The bifurcation diagram as a function of b ( 0 , 0.6 ) with the data set (49) and α = 0.1 . (a) The trajectory of S as a function of b . (b) The trajectory of I as a function of b . (c) The trajectory of Y as a function of b .
Mathematics 11 02909 g006
Figure 7. The trajectories of system (1) using data set (49) with different values of β and starting from various initial points. (a) Three-dimensional Phase plot converges m 1 = ( 4 , 0 , 0 ) . (b) Three-dimensional Phase plot converges m 2 = ( 3.74 , 0.17 , 0 ) . (c) Three-dimensional Phase plot converges m 4 = ( 0.8 , 0.58 , 0.16 ) . (d) Three-dimensional Phase plot converges m 2 = ( 0.33 , 0.66 , 0 ) .
Figure 7. The trajectories of system (1) using data set (49) with different values of β and starting from various initial points. (a) Three-dimensional Phase plot converges m 1 = ( 4 , 0 , 0 ) . (b) Three-dimensional Phase plot converges m 2 = ( 3.74 , 0.17 , 0 ) . (c) Three-dimensional Phase plot converges m 4 = ( 0.8 , 0.58 , 0.16 ) . (d) Three-dimensional Phase plot converges m 2 = ( 0.33 , 0.66 , 0 ) .
Mathematics 11 02909 g007
Figure 8. The trajectories of system (1) using data set (49) with different values of a 1 and starting from various initial points. (a) Three-dimensional Phase plot converges m 2 = ( 1 , 1.19 , 0 ) . (b) Three-dimensional Phase plot converges m 4 = ( 1.56 , 0.36 , 0.28 ) . (c) Three-dimensional Phase plot converges m 3 = ( 1.5 , 0 , 0.46 ) . (d) Three-dimensional Phase plot converges to a 3D periodic attractor. (e) Three-dimensional Phase plot converges to a 2D periodic attractor in the S Y -plane.
Figure 8. The trajectories of system (1) using data set (49) with different values of a 1 and starting from various initial points. (a) Three-dimensional Phase plot converges m 2 = ( 1 , 1.19 , 0 ) . (b) Three-dimensional Phase plot converges m 4 = ( 1.56 , 0.36 , 0.28 ) . (c) Three-dimensional Phase plot converges m 3 = ( 1.5 , 0 , 0.46 ) . (d) Three-dimensional Phase plot converges to a 3D periodic attractor. (e) Three-dimensional Phase plot converges to a 2D periodic attractor in the S Y -plane.
Mathematics 11 02909 g008aMathematics 11 02909 g008b
Figure 9. The trajectories of system (1) using data set (49) with different values of a 2 and starting from various initial points. (a) Three-dimensional Phase plot converges m 2 = ( 1 , 1.2 , 0 ) . (b) Three-dimensional Phase plot converges m 4 = ( 1.91 , 0.29 , 0.27 ) .
Figure 9. The trajectories of system (1) using data set (49) with different values of a 2 and starting from various initial points. (a) Three-dimensional Phase plot converges m 2 = ( 1 , 1.2 , 0 ) . (b) Three-dimensional Phase plot converges m 4 = ( 1.91 , 0.29 , 0.27 ) .
Mathematics 11 02909 g009
Figure 10. The trajectories of system (1) using data set (49) with different values of d 1 and starting from various initial points. (a) Three-dimensional Phase plot converges m 4 = ( 2.19 , 0.48 , 0.09 ) . (b) Three-dimensional Phase plot converges m 2 = ( 3 , 0.4 , 0 ) .
Figure 10. The trajectories of system (1) using data set (49) with different values of d 1 and starting from various initial points. (a) Three-dimensional Phase plot converges m 4 = ( 2.19 , 0.48 , 0.09 ) . (b) Three-dimensional Phase plot converges m 2 = ( 3 , 0.4 , 0 ) .
Mathematics 11 02909 g010
Figure 11. The trajectories of system (1) using data set (49) with different values of e 1 and starting from various initial points. (a) Three-dimensional Phase plot converges m 4 = ( 1.64 , 0.32 , 0.32 ) . (b) Three-dimensional Phase plot converges m 3 = ( 1.5 , 0 , 0.64 ) . (c) Three-dimensional Phase plot converges a 2D periodic attractor in the S Y -plane. (d) Three-dimensional Phase plot converges to a 3D periodic attractor.
Figure 11. The trajectories of system (1) using data set (49) with different values of e 1 and starting from various initial points. (a) Three-dimensional Phase plot converges m 4 = ( 1.64 , 0.32 , 0.32 ) . (b) Three-dimensional Phase plot converges m 3 = ( 1.5 , 0 , 0.64 ) . (c) Three-dimensional Phase plot converges a 2D periodic attractor in the S Y -plane. (d) Three-dimensional Phase plot converges to a 3D periodic attractor.
Mathematics 11 02909 g011aMathematics 11 02909 g011b
Figure 12. The trajectories of system (1) using data set (49) with different values of e 2 and starting from various initial points. (a) Time series of the solutions that converge m 2 = ( 1 , 1.19 , 0 ) . (b) Time series of the solutions that converge m 4 = ( 1.74 , 0.22 , 0.37 ) . (c) Time series of the solutions that represent bistable between a 3D periodic attractor in the S Y -plane and m 4 = ( 1.74 , 0.22 , 0.37 ) . (d) Time series of the solutions that represent bistable between a 3D periodic attractor in the S Y -plane and m 4 = ( 1.75 , 0.2 , 0.39 ) .
Figure 12. The trajectories of system (1) using data set (49) with different values of e 2 and starting from various initial points. (a) Time series of the solutions that converge m 2 = ( 1 , 1.19 , 0 ) . (b) Time series of the solutions that converge m 4 = ( 1.74 , 0.22 , 0.37 ) . (c) Time series of the solutions that represent bistable between a 3D periodic attractor in the S Y -plane and m 4 = ( 1.74 , 0.22 , 0.37 ) . (d) Time series of the solutions that represent bistable between a 3D periodic attractor in the S Y -plane and m 4 = ( 1.75 , 0.2 , 0.39 ) .
Mathematics 11 02909 g012
Figure 13. The trajectories of system (1) using data set (49) with different values of h 1 and starting from various initial points. (a) Time series of the solutions that converge m 4 = ( 1.14 , 0.54 , 0.07 ) . (b) Time series of the solutions that converge m 2 = ( 1 , 0.4 , 0 ) . (c) Time series of the solutions that converge m 1 = ( 0.5 , 0 , 0 ) . (d) Time series of the solutions that converge m 0 = ( 0 , 0 , 0 ) .
Figure 13. The trajectories of system (1) using data set (49) with different values of h 1 and starting from various initial points. (a) Time series of the solutions that converge m 4 = ( 1.14 , 0.54 , 0.07 ) . (b) Time series of the solutions that converge m 2 = ( 1 , 0.4 , 0 ) . (c) Time series of the solutions that converge m 1 = ( 0.5 , 0 , 0 ) . (d) Time series of the solutions that converge m 0 = ( 0 , 0 , 0 ) .
Mathematics 11 02909 g013aMathematics 11 02909 g013b
Figure 14. The trajectories of system (1) using data set (49) with different values of h 3 and starting from various initial points. (a) Three-dimensional Phase plot converges chaotic attractor. (b) Three-dimensional Phase plot converges m 3 = ( 1.4 , 0 , 0.66 ) . (c) Three-dimensional Phase plot converges m 4 = ( 1.66 , 0.29 , 0.33 ) (d) Three-dimensional Phase plot converges m 2 = ( 1 , 1.2 , 0 ) .
Figure 14. The trajectories of system (1) using data set (49) with different values of h 3 and starting from various initial points. (a) Three-dimensional Phase plot converges chaotic attractor. (b) Three-dimensional Phase plot converges m 3 = ( 1.4 , 0 , 0.66 ) . (c) Three-dimensional Phase plot converges m 4 = ( 1.66 , 0.29 , 0.33 ) (d) Three-dimensional Phase plot converges m 2 = ( 1 , 1.2 , 0 ) .
Mathematics 11 02909 g014aMathematics 11 02909 g014b
Figure 15. The bifurcation diagram as a function of h 3 ( 0 , 0.1 ) with the data set (49). (a) The trajectory of S as a function of h 3 . (b) The trajectory of I as a function of h 3 . (c) The trajectory of Y as a function of h 3 .
Figure 15. The bifurcation diagram as a function of h 3 ( 0 , 0.1 ) with the data set (49). (a) The trajectory of S as a function of h 3 . (b) The trajectory of I as a function of h 3 . (c) The trajectory of Y as a function of h 3 .
Mathematics 11 02909 g015
Table 1. The model terms and parameters description.
Table 1. The model terms and parameters description.
Parameter or TermDescription
r The susceptible prey’s intrinsic growth rate.
θ + μ ( 1 θ ) μ + Y [ 0 , 1 ] The fear function, where μ > 0 is the level of fear and θ [ 0 , 1 ] is the minimum cost of fear [36].
α The intraspecific competition rate among the prey individuals.
β The infection rate.
a 1 The maximum attack rate of the predator to the susceptible prey.
b The half saturation constant of the predator.
a 2 The maximum attack rate of the predator to the infected prey.
d i ; i = 1 , 2 The death rates of the infected prey and a predator, respectively.
e i [ 0 , 1 ] ; i = 1 , 2 The conversion rates of the susceptible and infected biomass to predator biomass.
h i ; i = 1 , 2 , 3 The harvested rates for the susceptible prey, infected prey, and predator.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Ibrahim, H.A.; Naji, R.K. The Impact of Fear on a Harvested Prey–Predator System with Disease in a Prey. Mathematics 2023, 11, 2909. https://doi.org/10.3390/math11132909

AMA Style

Ibrahim HA, Naji RK. The Impact of Fear on a Harvested Prey–Predator System with Disease in a Prey. Mathematics. 2023; 11(13):2909. https://doi.org/10.3390/math11132909

Chicago/Turabian Style

Ibrahim, Hiba Abdullah, and Raid Kamel Naji. 2023. "The Impact of Fear on a Harvested Prey–Predator System with Disease in a Prey" Mathematics 11, no. 13: 2909. https://doi.org/10.3390/math11132909

APA Style

Ibrahim, H. A., & Naji, R. K. (2023). The Impact of Fear on a Harvested Prey–Predator System with Disease in a Prey. Mathematics, 11(13), 2909. https://doi.org/10.3390/math11132909

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