361 

© 2025 The Author(s). Published by College of Education for Pure Science (Ibn Al-Haitham), 

University of Baghdad. This is an open-access article distributed under the terms of the Creative Commons 

Attribution 4.0 International License 

 
 

   

The Time Delay Effect and  Harvesting on  the  Predator-Prey: 

Analysis And Simulation 

 

Nidhal Faisal Ali 1,4* , Hassan F. Al-Husseiny 2  , and Yassine Sabbar 3  
1Department of  Mathematics, College of Science, University of Baghdad, Baghdad, Iraq. 
2Department of  Mathematics, College of Science, University of Baghdad, Baghdad, Iraq. 

3MAIS Laboratory, MAMCS Group, FST Errachidia, Moulay Ismail University of Meknes, Morocco. 
4Department of Electrical Engineering Techniques, College of  Electrical Engineering Technical,                      

University Middle Technical, Baghdad, Iraq. 

*Corresponding Author. 

 

Received:14 September 2023               Accepted:12 December 2023              Published: 20 April 2025  

doi.org/10.30526/38.2.3717 

 

Abstract    

     Mathematical modeling based on time-delay differential equations is an important tool to 

understand the effects of delays in biological systems and to analyze how they influence the 

dynamics of the asymptotic behavior of these systems. The prey-predator model described in 

this paper includes diseases in the prey species, harvests in each population, and time lags in 

predation and gestation of the predator. The solutions of the model are positive and bounded 

for all times within a realistic region. The existence of all fixed points has been proven. When 

a time lag is present, the essential conditions for the local stability of the positive equilibrium 

and the occurrence of Hopf bifurcations can be determined by analyzing the associated 

characteristic equation. The characteristics of the Hopf bifurcation are determined by applying 

normal form theory and the central manifold theorem. Finally, we use numerical simulations 

to validate our analytical results. A Hopf bifurcation in the system occurs when the delay 

exceeds a certain threshold. 

Keywords: harvest, predator–prey, time delay, stability, Hopf bifurcation. 

 

1. Introduction 

     The dynamics of interacting populations are often studied using mathematical models. 

Mathematical models have become essential tools to understand how diseases spread and are 

controlled. These models, referred to as "epidemiological models," are employed to study 

disease transmission and management in human or animal populations. On the other hand, the 

term "ecological model" pertains to mathematical models depicting the dynamic interactions 

among species in ecological systems. In 1925 and 1926, (1, 2) independently developed 

mathematical models within ecology, elucidating interactions among biological species. Eco-

epidemiological models encompass mathematical representations of dynamic behaviors within 

ecological systems, including disease dynamics. Anderson and May (3) studied the dynamics of 

an eco-epidemiology model that included interactions between infected prey populations and 

https://creativecommons.org/licenses/by/4.0/
https://creativecommons.org/licenses/by/4.0/
https://orcid.org/0000-0002-2136-7384
mailto:nidhal.f1980@yahoo.com
https://orcid.org/0009-0007-7883-5875
mailto:hassan.fadhil.r@sc.uobaghdad.edu%20.iq
https://orcid.org/0000-0002-1127-4395
mailto:y.sabbar@umi.ac.ma


IHJPAS. 2025,38(2) 

362 
 

 

predators in 1986. The researchers had studied the dynamics of the Eco-epidemiological 

models independently for many years, for example (4-6). It is well known that prey-predator 

models, which may involve a variety of natural factors, can be used to describe the most 

significant relationships between the individuals of species in the environment; modeling 

predator-prey interactions is becoming the most crucial topic for ecologists and applied 

mathematics research. Predator-prey models may often be divided into three types based on the 

infectious diseases that affect the population. The first type is that models only include infected 

prey (7-9). In 2012, a prey-predator model of the Lotka-Volterra type without harvesting was 

considered by Johri et al.(10). They assumed that there is a disease only in prey and that the 

susceptible prey's conversion rate is the same as that of the infected prey.  The eco-

epidemiological model with two prey populations where one prey species has an infectious 

disease is suggested and studied (11). The second type only includes the diseased in the predator 

(12-14) the prey-predator model with disease in the predator and the functional response studied 

and proposed (14). The third type of prey-predator model is a disease in both populations (15 -

17). Kant and Kumar (16)  investigated a prey-predator system where the prey migrates, and 

both populations encounter disease infection. The prey-predator model and assumed there is a 

disease in both populations have been studied (17). Harvesting strongly influences a model's 

behavior. It refers to reducing a population through hunting or individual capture. It may have 

a detrimental impact on harvested population density. In general, there are two kinds: linear 

and nonlinear. However, from a biological and economic viewpoint, nonlinear harvesting 

functions are better suited for use in reality. Research studies have used linear harvesting 

functions, constant-yield prey, or predator harvesting. In 2001, a mathematical model that 

included immature and mature supposes was studied by Song and Chen (18). The prey-predator 

system considers harvesting for the prey and predator species studied (19).The modeling of 

prey-predator involving the harvested with incorporating a prey refuge suggested and 

formulated (20). However, it is essential to consider the effect of past life history when 

analyzing the system's stability. Additionally, because time delays occur in so many biological 

situations (maturation, gestation, capture, or other factors.) in prey-predator systems and 

ignoring them means ignoring reality, delay differential equations are widely used in the 

literature on ecological interactions between predator and prey. On the other hand, several 

studies have been proposed to show the effect of the time delay (21-27). The influence of a 

delayed incubation period on disease transmission using a nonlinear incidence rate in prey-

predator systems was examined (21). Al-Jubouri and Naji (23) proposed a mathematical model 

incorporating a time delay in the disease transmission process. The effect of time delay on the 

dynamics of the prey-predator model with prey harvesting was studied  (25). In 2011, Naji and 

Ibrahim (28) formulated the following mathematical model. 
𝑑𝑆

𝑑𝑡 
= 𝑟𝑆 (1 −

𝑆+𝐼

𝐾
) − 𝐶𝑆𝐼 − 𝐸1𝑆

𝑑𝐼

𝑑𝑡
= 𝐶𝑆𝐼 −

𝛼𝐼𝑍

𝛾+𝐼
− 𝜆𝐼 − 𝐸2𝐼       

𝑑𝑍

𝑑𝑡
= −𝜃𝑍 +

ℎ𝐼𝑍

𝛾+𝐼
− 𝐸3𝑧              

                                                                                              (1) 

 where 𝑆(𝑡)  , 𝐼(𝑡) 𝑎𝑛𝑑 𝑍(𝑡) represent the number of susceptible prey, infected prey and 

predator respectively;  𝑟 intrinsic growth rate; 𝐾  carrying capacity of the prey in absence the 

predator and harvesting; 𝐶 infection rate; 𝛼 maximum attack rate; 𝜆 death rate of   𝐼(𝑡) ; 𝛾 half 

saturation level coefficient; 𝜃 death rate of 𝑍(𝑡) ; ℎ growth rate of the 𝑍(𝑡) due to predation of 

the   𝐼(𝑡);  𝐸1 , 𝐸2 𝑎𝑛𝑑 𝐸3  are the harvesting efforts for 𝑆(𝑡)  , 𝐼(𝑡) 𝑎𝑛𝑑 𝑍(𝑡) respectively. The 

time delay impact is concentrated. The organization of this paper is as follows: First, The given 



IHJPAS. 2025,38(2) 

363 
 

 

system (1) is modified. Second, it illustrates the positivity and bound of solutions of the system 

(2). Third, we essentially investigate the stability and existence of Hopf bifurcation. Fourth, 

investigates the properties of the Hopf bifurcation. Fifth, the major theoretical results and a 

discussion are shown by numerical simulation. Finally, the conclusion.   

Now, the improved system (1) can be expressed as follows. 
𝑑𝑆

𝑑𝑡 
= 𝑟𝑆 (1 −

𝑆+𝐼

𝐾
) − 𝐶𝑆𝐼 − 𝐸1𝑆                     

𝑑𝐼

𝑑𝑡
= 𝐶𝑆𝐼 −

𝛼𝐼(𝑡−𝜏)𝑍(𝑡−𝜏)

𝛾+𝐼(𝑡−𝜏)
− 𝜆𝐼 − 𝐸2𝐼             

𝑑𝑍

𝑑𝑡
= −𝜃𝑍 +

ℎ𝐼(𝑡−𝜏)𝑍(𝑡−𝜏)

𝛾+𝐼(𝑡−𝜏)
− 𝐸3𝑧                    

                                                                                   (2) 

Here, 𝜏 represents time delay. The system described above has the same biological interpretation 

for its parameters as those identified in the system (1). 

 

2. Materials and Methods 

2.1. Positive and Boundedness 

Before embarking on the study, it is essential to verify the biological integrity of the proposed 

model. Accordingly, within this section, we introduce the following theorem, which investigates 

the system's positivity and boundedness.  

Theorem 1. The solutions of system (2) are positive and bounded. 

Proof: First, it is proven that for 𝑡 ≥ 0   all solutions to system (2) are positive. 

𝑑𝑆

𝑑𝑡
≥ −𝑆 (

𝑟(𝑆 + 𝐼)

𝐾
+ 𝑐𝐼 + 𝐸1) 

Consequently, it is calculated to obtain it. 

𝑆(𝑡) ≥ 𝑆(0)𝑒𝑥𝑝 − {∫ (
𝑟(𝑆(ℴ)+𝐼(ℴ))

𝐾
+ 𝑐𝐼(ℴ) + 𝐸1) 𝑑(ℴ)

𝑡

0
}  

Because 𝑆(0) > 0. We get   𝑆(𝑡) > 0 for any 𝑆(0) > 0. 

The proof of  𝐼(𝑡) > 0 and 𝑍(𝑡) > 0 for all  𝑡 ≥ 0  can be done in the same way. 

Next, the following is the proof that the solutions of system (2) are bounded for all  𝑡 ≥ 0  . 

Define  𝒫(𝑡) = 𝑆(𝑡) + 𝐼(𝑡) + 𝑍(𝑡) 

As a result, the following is obtained: 
𝑑𝒫

𝑑𝑡
+𝓃𝒫 ≤ 𝑟𝑆, 

where 𝓃 = 𝑚𝑖𝑛 {𝐸1, 𝜆 + 𝐸2 , 𝜃 + 𝐸3}. 

Since 
𝑑𝑆

𝑑𝑡 
≤ 𝑟𝑆 (1 −

𝑆+𝐼

𝐾
) − 𝐸1𝑆  , we obtain according the comparison theorem (29)  

lim
𝑡→∞

𝑠𝑢𝑝 𝑆(𝑡) ≤
𝐾(𝑟 − 𝐸1)

𝑟
 

Thus, for 𝑡 ≥ 0 , we have 𝑆(𝑡) ≤
𝐾(𝑟−𝐸1)

𝑟
 . 

which implies that: 

 
𝑑𝒫

𝑑𝑡
+𝓃𝒫 ≤ 𝐾(𝑟 − 𝐸1), 

Then, after applying the Gronwall lemma to the above inequality, we have for 𝑡 → ∞ 

 𝒫(𝑡) ≤
𝐾(𝑟−𝐸1)

𝓃
 

This implies that   the solutions are bounded. 

2.2. Local  Stability and Hopf bifurcation  

 In this subsection, we will determine the local stability of each equilibrium point within system 

(2).  It is well-known that an equilibrium point's location and number unchanged with time 

delay.. As a result, system (2) has four equilibrium points. For further details, see (28) 



IHJPAS. 2025,38(2) 

364 
 

 

 The first equilibrium point, namely the vanishing equilibrium point      

           denoted  𝐹0 = (0,0,0) always exists. 

 The second equilibrium point, namely the axial equilibrium point      

           denoted  𝐹1 = (  
𝐾(𝑟−𝐸1)

𝑟
, 0,0) always exists. 

 The third equilibrium point, namely the planer equilibrium point      

            denoted  𝐹2 = (𝑆̅  , 𝐼,̅ 0) , where        

    𝑆̅  = 𝐸2 + 𝜆 𝐶⁄                                                                                                                       (3) 

    𝐼 ̅  = 𝐶𝐾(𝑟 − 𝐸1) − 𝑟(𝜆 + 𝐸2) (𝑟 + 𝐶𝐾)𝐶⁄                                                                           (4) 

It is clear that 𝐹2 is exists if the following condition satisfied  

         𝑟(𝜆 + 𝐸2) < 𝑐𝐾(𝑟 − 𝐸1)                                                                                                     (5) 

 The fourth equilibrium point, namely the positive equilibrium    

            point denoted  𝐹3 = (𝑆∗ , 𝐼∗, 𝑍∗) ,  

where   

          𝑆∗ =
ℎ𝐾(𝑟−𝐸1)−(𝜃+𝐸3)[𝐾(𝑟−𝐸1)−𝛾(𝑟+𝑐𝐾)]

𝑟(ℎ−(𝜃+𝐸3)
                                                                                     (6)                                                                       

𝐼∗ =
𝛾(𝐸3+𝜃)

ℎ−(𝜃+𝐸3)
; ℎ − (𝜃 + 𝐸3) ≠ 0                                                                                             (7) 

 𝑍∗ =
(𝛾+𝐼∗)(𝐶𝑆∗−(𝜆+𝐸3))

𝛼
                                                                                                             (8)          

 It is clear that 𝐹3 is exists if the following conditions are hold   

𝛾 <
𝐾(𝑟−𝐸1)(ℎ−(𝜃+𝐸3))

(𝑟+𝐶𝐾)(𝜃+𝐸3)
                                                                                                                (9) 

𝜃 + 𝐸3 < ℎ                                                                                                                              (10) 

  
𝜆+𝐸2

𝐶
< 𝑆∗                                                                                                                                (11) 

Now ,the linearization method  is used to determine the stability of the above equilibrium 

points. 

The general Jacobian matrix (𝐽𝑚)  for system (2) at any equilibrium point   𝐹 = (𝑆, 𝐼, 𝑍) is  

𝐽𝐹 =  [  𝑎𝑖𝑗] , 𝑖𝑗 = 0,1,2,3                                                                                                       (12) 

where  

𝑎11 = 𝑟 − (
𝑟

𝐾
(2𝑆 + 𝐼) + 𝐶𝐼 + 𝐸1)  ; 𝑎12 = −(

𝑟

𝐾
+ 𝐶) 𝑆   ; 𝑎13 = 0 ; 

𝑎21 = 𝐶𝐼 ; 𝑎22 =  𝐶𝑆 − ( 
𝛼𝛾𝑍𝑒−𝜇𝜏

(𝛾+𝐼)2
− (𝜆 + 𝐸2) ; 𝑎23 =

−𝛼𝐼𝑒−𝜇𝜏

𝛾+𝐼
; 

𝑎31 = 0; 𝑎32 =
ℎ𝛾𝑍𝑒−𝜇𝜏

(𝛾+𝐼)2
; 𝑎33 =

ℎ𝐼𝑒−𝜇𝜏

𝛾+𝐼
− (𝜃 + 𝐸3). 

Then , the characteristic equation corresponding to the matrix above can be expressed as:  

𝑞1(𝜇) + 𝑞2(𝜇)𝑒
−𝜇𝜏 = 0                                                                                                                 (13) 

where  𝑞1(𝜇) 𝑎𝑛𝑑 𝑞2(𝜇) are the polynomial of 𝜇 

The (𝐽𝑚)   for system (2) at 𝐹0is as follows:   

𝐽𝐹0 = [

𝑟 − 𝐸1 0 0
0 −(𝜆 + 𝐸2) 0
0 0 −(𝜃 + 𝐸3)

]                                                                             (14) 

The eigenvalues of  𝐽𝐹0  are 

𝜇01 = 𝑟 − 𝐸1,  𝜇02 = −(𝜆 + 𝐸2) < 0 and 𝜇03 = −(𝜃 + 𝐸3) < 0. 

The necessary condition of the coexistence of all species is  𝑟 − 𝐸1 > 0 see in (24)  

Obtained that  𝜇01 > 0, Thus  𝐹0 is unstable saddle point for all 𝜏 ≥ 0  

The (𝐽𝑚)  for system (2) at 𝐹1  is as follows: 

 



IHJPAS. 2025,38(2) 

365 
 

 

𝐽𝐹1 =

[
 
 
 −(𝑟 − 𝐸1) −

(𝑟−𝐸1)(𝑟+𝐶𝐾)

𝑟
0

0
𝐶𝐾(𝑟−𝐸1)

𝑟
− (𝜆 + 𝐸2) 0

0 0 −(𝜃 + 𝐸3)]
 
 
 

                                                         (15) 

 

The eigenvalues of  𝐽𝐹1  are  

𝜇11 = −(𝑟 − 𝐸1) < 0,  𝜇12 =
𝐶𝐾(𝑟−𝐸1)

𝑟
− (𝜆 + 𝐸2) and  𝜇13 = −(𝜃 + 𝐸3) < 0. 

It is clear that 𝐹1 is asymptotically stable for all 𝜏 ≥ 0  Provided that the following condition is 

met. 

𝐶𝐾(𝑟 − 𝐸1) < 𝑟(𝜆 + 𝐸2)                                                                                                       (16) 

The (𝐽𝑚)  for system (2) at  𝐹2is as follows: 

 𝐽𝐹2 =

[
 
 
 
 −

𝑟

𝐶𝐾
(𝜆 + 𝐸2)      −

(𝜆+𝐸2)(𝑟+𝐶𝐾)

𝐶𝐾
          0

𝜂

𝑟+𝐶𝐾
           0  − 

𝛼𝜂𝑒−𝜇𝜏

(𝛾+𝐼)̅(𝑟+𝐶𝐾)

0 0
ℎ 𝐼�̅�−𝜇𝜏

𝛾+ 𝐼̅
− (𝜃 + 𝐸3);

        

]
 
 
 
 

= 𝑏𝑖𝑗;  𝑖, 𝑗 =

1, . . ,3                                                                                                                                      (17) 

where 𝜂 = 𝐶𝐾(𝑟 − 𝐸1) − 𝑟(𝜆 + 𝐸2) 

Clearly, the roots of the following equation represent two eigenvalues of  𝐽𝐹2 

𝜇2 − (𝑏11 + 𝑏22)𝜇 + 𝑏11𝑏22 − 𝑏12𝑏21 = 0                                                                          (18) 

The above equation have negative real part for all 𝜏 ≥ 0  under the existence condition (5) of 

𝐹2 . 

While other eigenvalue  of  𝐽𝐹2 is given by the root of  

𝜇31 + (𝜃 + 𝐸3) −
ℎ 𝐼�̅�−𝜇𝜏

𝛾+ 𝐼̅
= 0                                                                                                                    (19) 

Thus, 

1- if  𝜏 = 0  equation  (19)   has eigenvalue  𝜇31 =
ℎ 𝐼̅

𝛾+ 𝐼̅
− (𝜃 + 𝐸3)   which is negative under 

the condition  
ℎ 𝐼̅

𝛾+ 𝐼̅
< (𝜃 + 𝐸3)                                                                                                                                           (20)     

  Hence the 𝐹2 is locally asymptotically stable under the conditions (5) and (20) when   𝜏 = 0   

2- for 𝜏 > 0,  if equation (19) has the roots which a pair of purely imaginary must 

intersect the imaginary axis, now let  𝜇 = 𝑖�̅�  (�̅� > 0)  be the root of equation (19) 

By substituting 𝜇 = 𝑖�̅�  in equation (19), we obtain  

𝜃 + 𝐸3 =
ℎ 𝐼 ̅

𝛾+ 𝐼 ̅
cos �̅� 𝜏

− �̅� =
ℎ 𝐼 ̅

𝛾+ 𝐼 ̅
sin �̅� 𝜏         

                                                                                          (21) 

Squaring and adding both sides of equation (21) ,we obtain the following result:   

�̅� = ±√(
ℎ 𝐼̅

𝛾+ 𝐼
)̅
2

− (𝜃 + 𝐸3)2                                                                                                (22) 

Note that, from the condition (20), we have �̅�(𝜏) when 𝜏 > 0  it cannot be real, which is in                  

opposition to the assumption. As a result, the root of the characteristic equation (19) cannot be 

purely imaginary, and is asymptotically stable for all. For all 𝜏 ≥ 0, the equilibrium point 𝐹2 

demonstrates asymptotic stability. 

The (𝐽𝑚)  for system (2) at  𝐹3 is as follows: 



IHJPAS. 2025,38(2) 

366 
 

 

𝐽𝐹3  = [

𝑟− 𝑅1 −𝑅2𝑆∗ 0
𝐶𝐼∗ 𝑅5 −  𝛼𝑅3𝑒

−𝜇𝜏 − 𝛼𝑅4𝑒−𝜇𝜏

0  ℎ𝑅3𝑒
−𝜇𝜏  ℎ𝑅4𝑒

−𝜇𝜏 − 𝑅6

 ] = 𝑐𝑖𝑗;  𝑖, 𝑗 = 1, . .3                           (23)                                                                                           

where  

𝑅1 =
𝑟

𝐾
(2𝑆∗ + 𝐼∗) + 𝐶𝐼∗ + 𝐸1 > 0  ; 𝑅2 =

𝑟

𝐾
+ 𝐶 > 0   ; 𝑅3 =

𝛾𝑍∗

(𝛾+𝐼∗)2
> 0; 

 𝑅4 =
𝐼∗

𝛾+𝐼∗
   ; 𝑅5 = 𝐶𝑆∗ − (𝜆 + 𝐸3) > 0 ; 𝑅6 = (𝜃 + 𝐸3)> 0. 

The characteristic equation of  𝐽𝐹3  is   

𝜇3 +𝒜1𝜇
2 +𝒜2𝜇+𝒜3 + (𝒜4𝜇

2 +𝒜5𝜇+𝒜6)𝑒
−𝜇𝜏 = 0                            (24) 

Where  

𝒜1 = 𝑅1 + 𝑅6 − (𝑟 + 𝑅5)   

𝒜2 = (𝑟− 𝑅1)(𝑅5 − 𝑅6) + 𝑅2𝐶 𝑆
∗𝐼∗ − 𝑅5𝑅6 

𝒜3 = (𝑟− 𝑅1)𝑅5𝑅6+ 𝑅2𝑅6𝐶 𝑆
∗𝐼∗ 

𝒜4 =  𝛼𝑅3 − ℎ𝑅2 

𝒜5 = ℎ𝑅5𝑅6 +  𝛼 𝑅3𝑅6  + (𝑟− 𝑅1)( ℎ𝑅4− 𝛼𝑅3) 

𝒜6 = − ℎ𝑅4(𝑅2𝐶 𝑆
∗𝐼∗ + 𝑅5(𝑟− 𝑅1))− 𝛼𝑅3𝑅6𝑒

−𝜇𝜏(𝑟− 𝑅1) 

Thus,  

1. If 𝜏 = 0, then equation (24) becomes: 

𝜇3 +  𝜇2(𝒜1 +𝒜4) + 𝜇(𝒜2 +𝒜5) + (𝒜3 +𝒜6) = 0                                    

(25) 

According to the Hurwitz criterion, equation (25) have three negative roots if the following 

conditions are satisfied. 

 𝑟 <  𝑅1                                                                                                                                 (26) 

 𝑅5 < 𝛼 𝑅3                                                                                                                              (27)             

 ℎ𝑅4 < 𝑅6                                                                                                                               (28) 

 𝑅8 <  𝑅7                                                                                                                                 (29) 

Here  

 𝑅7 = ( 𝑐11 +  𝑐22 +  𝑐33)[− 𝑐11( 𝑐22 +  𝑐33)−  𝑐22 𝑐33 +  𝑐23 𝑐32 +  𝑐12 𝑐21] 
 𝑅8 = −[ 𝑐11( 𝑐22 𝑐33−  𝑐23 𝑐32) −  𝑐33 𝑐12 𝑐21] 
and  

𝑐𝑖𝑗;  𝑖, 𝑗 = 1, . .3 define in equation (23) 

Hence the 𝐹3 is locally asymptotically stable under the conditions (26-29) 

2- When 𝜏 > 0, Assume that the root of equation (24) is purely imaginary, namely 𝜇 = ±𝑖𝜔  

  (  𝜔 > 0) if in addition to condition (26, 27) and the following condition hold 

𝒜6 > 𝒜3                                                                                                              (30)  

Let  𝜇 = 𝑖𝜔  is the root of equation (24) and by separating equation (24) to the real and 

imaginary part, yields 

 
( 𝒜4𝜔

2 −𝒜6 ) 𝑠𝑖𝑛𝜔𝜏 +𝒜5 𝜔 𝑐𝑜𝑠 𝜔𝜏 = 𝜔
3 −𝒜2 𝜔   

𝒜5 𝜔 𝑠𝑖𝑛𝜔𝜏  +( 𝒜6 −𝒜4𝜔
2 ) 𝑐𝑜𝑠 𝜔𝜏 =  𝒜1 𝜔

2 −𝒜3 

                                                    (31) 

The result of squaring and summing the above equations yields:  

  𝜔6 + 𝒞1𝜔
4 + 𝒞2𝜔

2 + 𝒞3 = 0                                                                                 (32) 

Where 

𝒞1 = 𝒜1
2   − 𝒜4

2 − 2𝒜2; 

𝒞2 = 𝒜2
2   − 𝒜5

2 − 2𝒜1𝒜3 + 2 𝒜4 𝒜6 ; 



IHJPAS. 2025,38(2) 

367 
 

 

𝒞3 = 𝒜3
2   − 𝒜6

2 

Put  ℋ = 𝜔2 , then equation (32) becomes 

ℋ3 + 𝒞1ℋ
2 + 𝒞2ℋ + 𝒞3 = 0                                                                                  (33)  

 Under the conditions (27,28) and condition (30) we have   𝒞3 < 0 . The equation (33) have a 

positive root which is unique say 𝜔0   by using Descartes' rule of sign. Hence 𝜔0 is also the 

positive root of equation (32). Hence, there exists at least a pair of imaginary roots, denoted as 

±𝑖𝜔0 , that satisfies equation (24). From equation (31) after substituting 𝜔0 , we obtain : 

 cos  𝜔0𝜏 =
𝒷1

𝒷2 
 

Here 

𝒷1 = (  𝒜5
−  𝒜1 𝒜4 )𝜔0

4 + ( 𝒜1 𝒜 6 +  𝒜3 𝒜 4 −  𝒜2 𝒜5 )𝜔0
2  −  𝒜3 𝒜6  

 𝒷2 =   𝒜 4
2𝜔0

4 + (  𝒜 5
2 −  2𝒜4 𝒜6 )𝜔0

2  +  𝒜 6
2
 

Then, 𝜏𝑚  corresponding to 𝜔0  as below 

𝜏𝑚 =
1

 𝜔0 
(cos−1(

ℓ1

ℓ2
) + 2𝜋𝑚) ;𝑚 = 0,1,2, …                                                                        (34) 

Define  𝜏0 = 𝑚𝑖𝑛
𝑚≥0

𝜏𝑚 

Then, we obtain the following theorem 

Theorem 2. The system (2)  asymptotic stability at  𝐹3within the τ ∈ [0, τ0) and a Hopf 

bifurcation occurs at 𝜏 = τ0    when specific conditions are met. 

3(𝜔0
2)2 + 2𝒞1𝜔0

2 + 𝒞2 ≠ 0                                                                                                        (35) 

Proof.  

for 𝜏 ∈  [0, τ0)    𝐹3  is asymptotically stable as shown in conditions (26) - (29). 

However, when 𝜏 = τ0, we can demonstrate the presence of a Hopf bifurcation by establishing 

that 𝐹3  is conditionally stable, specifically, by confirming that equation (24) exhibits purely 

imaginary roots ±𝑖 𝜔0    at 𝜏 = τ0, This condition can be expressed [
𝑑(𝑅𝑒𝜆(𝜏))

𝑑𝜏
]
𝜏=𝜏0

≠ 0  

If we assume that equation (24) has the eigenvalue that is 𝜇(𝜏) = ð(𝜏) + 𝑖𝜔(𝜏) such that 

 ð(τ0) = 0 and  𝜔(τ0) = 𝜔0 > 0. τ0 define in equation (34).  

When we differentiate equation (24) with respect to 𝜏, and  apply the chain rule. This results in 

the following expression: 

[3𝜇2 + 2 𝒜1𝜇 +  𝒜2 + (2 𝒜4𝜇 +  𝒜5)𝑒
−𝜇𝜏 − 𝜏( 𝒜4𝜇

2 +  𝒜5𝜇 +  𝒜6)𝑒
−𝜇𝜏]

𝑑𝜇

𝑑𝜏
 

= 𝜇( 𝒜4𝜇
2 +  𝒜5𝜇 +  𝒜6)𝑒

−𝜇𝜏                                                                               (36) 

From equation (24), we have  

[
𝑑𝜇

𝑑𝜏
]
−1

=
(3𝜇2+2 𝒜1𝜇+ 𝒜2)𝑒

𝜇𝜏

𝜇( 𝒜4𝜇2+ 𝒜5𝜇+ 𝒜6)
+

2 𝒜4𝜇 + 𝒜5

𝜇( 𝒜4𝜇2+ 𝒜5𝜇+ 𝒜6)
−

𝜏

𝜇
                                                            (37)     

Since for 𝜏 = 𝜏0,  and  𝜇 = 𝑖𝜔0  we get  

[
𝑑𝜇

𝑑𝜏
]
−1

=
(( 𝒜2 − 3𝜔0

2) + 2𝑖 𝒜1𝜔0) (cos𝜔0𝜏 + 𝑖 sin𝜔0𝜏) 

− 𝒜5𝜔0
2 + 𝑖𝜔0( 𝒜6 −  𝒜4𝜔0

2)
+

 𝒜5 + 2𝑖 𝒜4𝜔0

− 𝒜5𝜔0
2 + 𝑖𝜔0( 𝒜6 −  𝒜4𝜔0

2)

−
𝜏0
𝑖𝜔0 

 

Now since 

𝑠𝑖𝑔𝑛 [
𝑑(𝑅𝑒𝜇)

𝑑𝜏
]
𝜏=𝜏0

= 𝑠𝑖𝑔𝑛 [𝑅𝑒( 
𝑑𝜇

𝑑𝜏
 )−1]

 𝜇=𝑖𝜔0
  .                                                                      (38)                                                                 

It is clear that : 



IHJPAS. 2025,38(2) 

368 
 

 

[
𝑑𝜇

𝑑𝜏
]
−1

=
(( 𝒜2 − 3𝜔0

2) + 2𝑖 𝒜1𝜔0) (cos𝜔0𝜏 + 𝑖 sin𝜔0𝜏) 

− 𝒜5𝜔0
2 + 𝑖𝜔0( 𝒜6 −  𝒜4𝜔0

2)
+

 𝒜5 + 2𝑖 𝒜4𝜔0

− 𝒜5𝜔0
2 + 𝑖𝜔0( 𝒜6 −  𝒜4𝜔0

2)
−

𝜏0
𝑖𝜔0 

 

Hence, we have 

𝑅𝑒 [
𝑑𝜇

𝑑𝜏
]
𝜏=𝜏0

−1

= 𝑅𝑒
(( 𝒜2−3𝜔0

2)+2𝑖 𝒜1𝜔0) (cos𝜔0𝜏+𝑖 sin𝜔0𝜏) 

− 𝒜5𝜔0
2+𝑖𝜔0( 𝒜6− 𝒜4𝜔0

2)
+ 𝑅𝑒

 𝒜5+2𝑖 𝒜4𝜔0

− 𝒜5𝜔0
2+𝑖𝜔0( 𝒜6− 𝒜4𝜔0

2)
  

𝑅𝑒 [
𝑑𝜇

𝑑𝜏
]
𝜏=𝜏0

−1

=
𝜔0
2[3𝜔0

4 + (2𝒜1
2 − 4𝒜2 − 2𝒜4

2)𝜔0
2 + (𝒜2

2 − 2𝒜1𝒜3 −𝒜5
2 + 2𝒜4𝒜6)]

(𝒜5𝜔0
2)2 + (𝒜6 −𝒜4𝜔0

2)2

=
𝑓′(𝜔0

2)

(𝒜5𝜔0
2)2 + (𝒜6 −𝒜4𝜔0

2)2
 

 

Here  𝑓′(𝜔0
2) = 3(𝜔0

2)2 + 2𝒞1𝜔0
2 + 𝒞2 ≠ 0 due to condition (35). So, we have  

𝑆𝑖𝑔𝑛 {
𝑑

𝑑𝜏
(𝑅𝑒𝜇)|𝜏=𝜏0} = 𝑆𝑖𝑔𝑛 {𝑅𝑒 (

𝑑𝜇

𝑑𝜏
)
𝜏=𝜏0

−1

} = 𝑆𝑖𝑔𝑛{𝑓′(𝜔0
2)}  

Assuming that 
𝑑

𝑑𝜏
(𝑅𝑒𝜇)|𝜏=𝜏0 < 0,  it is implies that the roots of the characteristic has   roots 

with positive real parts at 𝜏 = 𝜏0. This contradicts the local stability of the positive equilibrium 

point .Therefore  , we can deduce that [
𝑑(𝑅𝑒𝜇)

𝑑𝜏
]
𝜏=τ0

> 0 under condition (35).Consequently ,  

the transversality condition is satisfied , leading  to  a Hopf bifurcation happens at 𝜏 = 𝜏0,  

and  𝜇 = 𝑖𝜔0. 
2.3. The Direction and  Stability of  the Hopf  Bifurcation. 

In this section, we investigate the orientation of the Hopf bifurcation near 𝐹3  at 𝜏 = 𝜏0and 

establish the prerequisites for the stability of the resulting periodic solution in the system (2). 

We accomplish this by applying Hassard's center manifold theorem and normal form theory 

)30(. 

Theorem 3. 

     )i)  If   ℳ2 > 0 , then the Hopf bifurcation is supercritical and the bifurcating periodic 

solutions   

           exist for 𝜏 > 𝜏0  , and   If  ℳ2 < 0 ,then the Hopf bifurcation is subcritical  and the      

             bifurcating  periodic  solutions exist for  𝜏 < 𝜏0   .     

     (ii)   If  𝒰2 < 0, then the bifurcating periodic solution are stable, and if  𝒰2 > 0, then the  

              bifurcating  periodic solution are unstable 

     (iii)  If 𝒯2 > 0, the period of the bifurcating cyclic solutions increases, and if   𝒯2 < 0, the   

                period decreases.  

        where ℳ2, 𝒰2and 𝒯2 are given 

𝐶1(0) =
𝑖

2𝜔0𝜏0
(𝒢11  𝒢20 − 2|𝒢11|

2 −
|𝒢02|

2

3
) +

𝒢21

2
,

ℳ2 = −
𝑅𝑒{𝐶1(0)}

𝑅𝑒{
𝑑𝜇

𝑑𝜏
(𝜏0)}

,                                                          

𝒰2 = 2𝑅𝑒{𝐶1(0)},                                                       

𝒯2 =
−𝐼𝑚{𝐶1(0)}+ℳ2  𝐼𝑚{

𝑑𝜇

𝑑𝜏
(𝜏0)}

𝜔0𝜏0
.                                    }

  
 

  
 

                                                        (39)                                                          

and 𝒢11, 𝒢20, 𝒢02 and 𝒢21 are given in the proof 

Proof. Let 𝔘1(𝑡) = 𝑆(𝑡) − 𝑆
∗, 𝔘2(𝑡) = 𝐼(𝑡) − 𝐼∗, 𝔘3(𝑡) = 𝑍(𝑡) − 𝑍

∗, and 𝜏 = 𝜏0 + 𝒮, here  

𝜏0 is define by equation  (34 ) and 𝒮 ∈ 𝑅  . It is possible to convert system (2) into a 

functional differential equation in 𝐶 = 𝐶([−1,0],𝑅3)  then 

𝔘′(𝑡) = 𝐿𝒮(𝔘𝑡) + ℱ(𝒮, 𝔘𝑡),                                                                                                  (40) 



IHJPAS. 2025,38(2) 

369 
 

 

 Where 

𝔘(𝑡) = (𝔘1(𝑡), 𝔘2(𝑡), 𝔘3(𝑡))
𝑇 ∈ 𝐶 = 𝐶([−1,0], 𝑅3) and 𝐿𝒮: 𝐶 → 𝑅3, ℱ:𝑅 × 𝐶 → 𝑅3   

are given by:  

𝐿𝒮(Γ) = (𝒮 + 𝜏0)[𝔇1Γ(0) + 𝔇2Γ(−1)]                                                                                (41) 

 The nonlinear is  

 ℱ(𝒮, Γ) = (𝒮 + 𝜏0) (
ℋ1

ℋ2

ℋ3

) 

where 

𝔇1 = [

ℱ10
(1)

ℱ01
(1)

 0

ℱ1000
(2)

ℱ0100
(2)

0

0 0 ℱ100
(3)

] = [  

  𝑟 − 𝑅1  −𝑅2𝑆∗                     0
  𝐶𝐼∗      𝑅5                    0     
    0        0                 −𝑅6

],             

         

𝔇2 = [

0 0 0

ℱ0010
(2)

0 ℱ0001
(2)

ℱ010
(3)

0 ℱ001
(3)

] = [

0 0 0
−𝛼𝑅3 0 − 𝛼𝑅4
 ℎ𝑅3 0 ℎ𝑅4  

], 

 

with 𝑅1, 𝑅2, 𝑅3  and 𝑅4 are define in the 𝐽𝐹3, while 

ℋ1 = ∑
1

𝑖!𝑗!𝑖+𝑗≥2 ℱ𝑖𝑗
(1)Γ1

𝑖(0)Γ2
𝑗(0), 

ℋ2 = ∑
1

𝑖! 𝑗!𝑚! 𝑛!
𝑖+𝑗+𝑚+𝑛≥2

ℱ𝑖𝑗𝑚𝑛
(2) Γ1

𝑖(0)Γ2
𝑗(0)Γ̃1

𝑚(−1)Γ̃2
𝑛(−1), 

ℋ3 = ∑
1

𝑘!𝑚! 𝑛!
𝑘+𝑚+𝑛≥2

ℱ𝑘𝑚𝑛
(3) Γ3

𝑘(0)Γ̃1
𝑚(−1)Γ̃2

𝑛(−1), 

where, Γ(𝜐0) = (Γ1(𝜐0), Γ2(𝜐0), Γ3(𝜐0)) ∈ 𝐶,−1 ≤ 𝜐0 ≤ 0,  and  

ℱ𝑖𝑗
(1)𝛤1

𝑖(0)𝛤2
𝑗(0), =

𝜕𝑖+𝑗ℱ(1)

𝜕𝛤1
𝑖𝛤2
𝑗 |

(𝛤1,𝛤2)=(0,0)

, 

ℱijmn
(2) Γ1

i(0)Γ2
j(0)Γ̃1

m(−1)Γ̃2
n(−1) =

∂i+j+m+nℱ(2)

∂Γ1
i Γ2

j
Γ̃1
mΓ̃2

n
|
(Γ1,Γ2,Γ̃1,Γ̃2)=(0,0,−1,−1)

,         

ℱ𝑘𝑚𝑛
(3) Γ3

𝑘(0)Γ̃1
𝑚(−1)Γ̃2

𝑛(−1) =
𝜕𝑘+𝑚+𝑛ℱ(3)

𝜕Γ3
𝑘Γ̃1

𝑚Γ̃2
𝑛 |

(Γ3,Γ̃1,Γ̃2)=(0,−1,−1)
. 

Based on the Riesz representation theorem, a 3×3 matrix function ℳ0 (𝜐0, 𝒮) exists for −1 ≤

𝜐0 ≤ 0  such that. 

              𝐿𝒮(Γ) = ∫   𝑑
0

−1
ℳ0(𝜐0,𝒮)Γ(𝜐0)  𝑓𝑜𝑟  Γ𝜖𝐶.                                                                             (42)                              

In actuality, it can be chosen. 

ℳ0(𝜐0,𝒮) = (𝜏0 + 𝒮) (𝔇1𝜎(𝜐0)−𝔇2𝜎(𝜐0+ 1)) ,                                                                (43) 

here, 𝜎 is  called the Dirac delta function and  

𝜎(𝜐0) = {
1 𝜐0 = 0
0 𝜐0 ≠ 0

    .     

For  Γ ∈ 𝐶([−1,0], 𝑅3), define        

 

𝒜(𝒮)Γ(𝜐0) = {

𝑑Γ(𝜐0)

𝑑𝜐0
 ,                                  −1 ≤ 𝜐0 < 0  ,

∫  𝑑
0

−1
𝜂
0
(ℌ0,𝜎0)𝜑0(ℌ0),        𝜐0 = 0   ,

                                                (44) 

and  



IHJPAS. 2025,38(2) 

370 
 

 

ℛ(𝒮)Γ(𝜐0) = {
0,               −1 ≤ 𝜐0 < 0,      

ℱ(𝒮 , Γ), 𝜐0 = 0.
                                                                         (45) 

 Hence, the system (39) is equivalent  

𝔘′(𝑡) = 𝒜(𝒮)𝔘𝑡 + ℛ(𝒮)𝔘𝑡  .                                                                                                (46) 

Where, 𝔘𝑡 = 𝔘(t + 𝜐0),−1 ≤ 𝜐0 ≤ 0 .                             

For 𝛹0 ∈ 𝐶1([−1,0], 𝑅3), the adjoint operator 𝒜∗of 𝒜(0) is  

𝒜∗ 𝛹0(ℌ0) = {
−
𝑑𝛹0(ℌ0)

𝑑ℌ0
 ,                                 0 < ℌ0 ≤ 1 ,     

∫ 𝑑
0

−1
Γ𝑇(ℰ0, 0)𝛹0(−ℰ0), ℌ0 = 0 .

                                            (47) 

                                                            

For Γ ∈ 𝐶([−1,0], 𝑅3), and   𝛹0 ∈ (𝐶1[−1,0], (𝑅3)∗) . we define the bilinear inner product 

  〈𝛹0(ℌ0), Γ(𝜐0)〉 = 𝛹0(0)Γ(0) − ∫ ∫ 𝛹0
𝑇
(𝜍0 − 𝜐0) 𝑑ℳ0(𝜐0)Γ(𝜍0)𝑑𝜍0 ,

𝜐0

𝜍0=0

0

𝜐0=−1
                   (48)                                         

Given ℳ0(𝜐0) =  ℳ0(𝜐0, 0., it follows that 𝒜 = 𝒜 (0) and 𝒜∗ are adjoint operators. Referring 

to the previous theorem 2, we can deduce that ±𝑖 𝜔0     are eigenvalues of 𝒜 (0) and 𝒜∗, 

respectively. By performing a straightforward calculation, it becomes evident that. 

                         𝑝(𝜐0) = (1, 𝑝1, 𝑝2)
𝑇 𝑒𝑖𝑤0𝜏0𝜐0 

                         𝑝∗(ℌ0) = 𝐷0(1, 𝑝1
∗, 𝑝2

∗)𝑇  𝑒−𝑖𝑤0𝜏0ℌ0 

Here  

                  𝑝1 =
𝑖𝑤0−ℱ10

(1)

ℱ01
(1) ,     𝑝2 = −

ℱ1000
(2)

+ℱ0010
(2)

𝑒−𝑖𝑤0𝜏0+(ℱ0100
(2)

  −𝑖𝑤0)𝑝1

𝑓0001
(2)

  𝑒𝑖𝑤0𝜏0
, 

                  𝑝1
∗ = −

ℱ10
(1)

ℱ0100
(2)

++𝑖𝑤0
, 𝑝1

∗ = −
ℱ10
(2)
+𝑖𝑤0+(ℱ1000

(2) +ℱ0010
(2)

  𝑒−𝑖𝜏0𝑤0)𝑝1
∗

ℱ010
(2)

  𝑒−𝑖𝑤0𝜏0
. 

From equation (48), we can get  

 
⟨𝑝∗(ℌ0), 𝑝(𝜐0)⟩ = 𝐷0̅̅ ̅ [1 + �̅�1

∗𝑝1 + �̅�2
∗𝑝2 + �̅�1

∗𝜏0𝑒
−𝑖𝑤0𝜏0 (ℱ00010

(2)
+ℱ00001

(2)
𝑝2)

+�̅�2
∗𝜏0𝑒

−𝑖𝑤0𝜏0 (ℱ010
(3)

+ℱ001
(3)
𝑝2)] .

                  (49) 

Let, 𝐷0 = [1 + 𝑝1
∗ �̅�1 + 𝑝2

∗ �̅�2 + 𝜏0  𝑒
−𝑖𝑤0𝜏0 [𝑝1

∗ (ℱ00010
(2)

+ℱ00001
(2)

�̅�2)+ 𝑝2
∗ (ℱ010

(3)
+

ℱ001
(3)
�̅�2)]]

−1

, where, 𝐷0̅̅ ̅  represent the conjugate complex number of 𝐷0, such that ⟨𝑝∗, 𝑝⟩ = 1 

and 

 ⟨𝑝∗, �̅�⟩ = 0. 

Next, using Hassard et al.'s algorthems )26(, we can obtain the Hopf bifurcation's properties:  

𝒢
20
= 2𝜏0𝐷0̅̅ ̅(ℏ1 + ℏ5�̅�1

∗ + ℏ9�̅�2
∗)

𝒢
11
= 𝜏0𝐷0̅̅ ̅(ℏ2 + ℏ6�̅�1

∗ + ℏ10�̅�2
∗)

𝒢
02
= 2𝜏0𝐷0̅̅ ̅(ℏ3 + ℏ7�̅�1

∗ + ℏ11�̅�2
∗)

𝒢
21
= 2𝜏0𝐷0̅̅ ̅(ℏ4 + ℏ8�̅�1

∗ + ℏ12�̅�2
∗)}
 
 

 
 

                                                                            (50)                                                      

Where 

ℏ1 = ℱ11
(1)
  𝑃1 +ℱ20

(1)
  , 

ℏ2 = ℱ11
(1)
  (𝑃1 + �̅�1  ) + 2ℱ20

(1)
  , 

ℏ3 = ℱ11
(1)
  �̅�2 +ℱ20

(1)
  , 

ℏ4 = ℱ11
(1)
( 𝑃1 𝑤11

(1)
(0) +

1

2
  �̅�1  𝑤20

(1)
(0) +

1

2
  𝑤20

(2)
(0) + 𝑤11

(2)
(0))

+ℱ02
(1)
( 𝑤20

(1)
(0) + 2 𝑤11

(1)
(0))

, 

ℏ5 = ℱ1100
(2)

  𝑃1 +ℱ0020
(2)

  𝑃1
(2)
𝑒−2𝑖𝑤0𝜏0 +ℱ0011

(2)
  𝑃1𝑃2  𝑒

−2𝑖𝑤0𝜏0, 



IHJPAS. 2025,38(2) 

371 
 

 

ℏ6 = ℱ1100
(2)

(𝑃1 + �̅�1 ) + 2ℱ0020
(2)

  𝑃1�̅�1 + 2 𝑃1𝑃2ℱ0011
(2)

    , 

ℏ7 = ℱ1100
(2)

  �̅�1 + 𝑓0020 
(2)

�̅�1
2𝑒2𝑖𝑤0𝜏0   +   �̅�1�̅�2𝑒

2𝑖𝑤0𝜏0, 

ℏ8 = ℱ1100
(2)

( 
1

2
�̅�1  𝑤20

(1)
(0) + 𝑃1𝑤11

(1)
(0) +

1

2
𝑤20
(2)
(0) + 𝑤11

(2)
(0))                                         

  +ℱ0020
(2)

(�̅�1  𝑤20
(2)
(−1)𝑒𝑖𝑤0𝜏0 + 2𝑃1  𝑤11

(2)
(−1)𝑒−𝑖𝑤0𝜏0)             

                            +ℱ0011
(2)

(
1

2
(𝑃̅̅ 1̅   + �̅�2)𝑤20

(2)
(−1)  𝑒𝑖𝑤0𝜏0 + (𝑃1 + 𝑃2)  𝑤11

(3)
(−1)  𝑒−𝑖𝑤0𝜏0)

 

ℏ9 = ℱ020
(3)
  𝑃1

2  𝑒−2𝑖𝑤0𝜏0 +ℱ011
(3)
  𝑃1𝑃2  𝑒

−2𝑖𝑤0𝜏0, 

ℏ10 = 2ℱ020
(3)
  𝑃1�̅�1 + 2ℱ011

(3)
  𝑃1𝑃2, 

ℏ11 = ℱ020
(3)
  �̅�1

2  𝑒2𝑖𝑤0𝜏0 +ℱ011
(3)
 �̅�1�̅�2  𝑒

2𝑖𝑤0𝜏0, 

ℏ12 = ℱ020
(3)
(  �̅�1  𝑤20

(1)
(−1)  𝑒𝑖𝑤0𝜏0 + 2𝑃1𝑤11

(1)
(−1)𝑒−𝑖𝑤0𝜏0)                                     

+ℱ011
(3)

(
1

2
(�̅�1   + �̅�2)𝑤20

(2)
(−1)  𝑒𝑖𝑤0𝜏0 + (𝑃1 + 𝑃2)  𝑤11

(3)
(−1)  𝑒−𝑖𝑤0𝜏0)

 

with 

 𝑤20(𝜐0) =
𝑖𝒢
20

𝑤0𝜏0
𝑃(0)𝑒𝑖𝑤0𝜏0𝜐0 +

𝑖�̅�
02

3𝑤0𝜏0
�̅�(0)𝑒−𝑖𝑤0𝜏0𝜐0 + ℒ1𝑒

2𝑖𝑤0𝜏𝜐0.                                     (51) 

𝑤11(𝜐0) = −
𝑖𝒢11
𝑤0𝜏0

  𝑃(0)  𝑒𝑖𝑤0𝜏0𝜐0 +
𝑖�̅�11
𝑤0𝜏0

�̅�(0)  𝑒−𝑖𝑤0𝜏0𝜐0 + ℒ2.                                            (52) 

             

Her ℒ1 = (ℒ1
(1)
, ℒ1

(2)
, ℒ1

(3)
)
𝑇

 and ℒ2 = (ℒ2
(1)
, ℒ2

(2)
, ℒ2

(3)
)
𝑇

 can be calculated by the following  

equations: 

𝒿1
∗ℒ1 = 2𝜏0𝒿1.                                                                                                        (53)                                                                                              

𝒿2
∗  ℒ2 = −𝜏0𝒿2 .                                                                                                                      (54) 

Where, 

𝒿1
∗ = (2𝑖𝑤0𝜏0𝐼 − ∫

0

−1

𝑑ℳ0(𝜐0)  𝑒
2𝑖𝑤0𝜏0𝜐0), 

𝒿2
∗ = (∫

0

−1

𝑑ℳ0(𝜐0)), 

𝒿1 = (ℏ1  ℏ5  ℏ9)
𝑇,  

𝒿2 = (ℏ2  ℏ6  ℏ10)
𝑇. 

Accordingly, it is determined that: 

 ℒ1 =

 2 𝒿1

(

 
 
 

2𝑖𝑤0 − ℱ10
(1)
    −ℱ01

(1)
  0

−ℱ10000
(2)

− ℱ0010
(2)

𝑒2𝑖𝑤0𝜏0𝜐0    2𝑖𝑤0 − ℱ0100
(2)

    ℱ0001
(2)

𝑒2𝑖𝑤0𝜏0𝜐0

−ℱ010
(3)
  𝑒2𝑖𝑤0𝜏0𝜐0     0    2𝑖𝑤0 − ℱ100

(3)
− ℱ001

(3)
𝑒2𝑖𝑤0𝜏0𝜐0)

 
 
 

−1

. 

  

ℒ1 = − 𝒿2

(

 
 
 

−ℱ10
(1)
    −ℱ01

(1)
     0

−ℱ1000
(2)

−ℱ0010
(2)

    −ℱ0100
(2)

     −ℱ0001
(2)

  

−ℱ010
(3)
    0    −ℱ100

(3)
−ℱ001

(3)
)

 
 
 

−1

      . 



IHJPAS. 2025,38(2) 

372 
 

 

Thus, equations (51)– (54) can be used to calculate 𝑤20(𝜐0) and 𝑤11(𝜐0) . Following that, the 

proof can be completed by determining the expressions given in equation (39) based on those 

in equation (50).  

 

 3. Numerical Simulation  

     In this section, we present essential findings through numerical representation, utilizing a set 

of biologically reasonable hypothetical values as provided below. The goal is to validate the 

theoretically generated outcomes and gain insights into the parameters' impact on the system 

dynamics (2). 

   𝑟 = 0.009 , 𝑘 = 40 , 𝑐 = 0.0002 , 𝐸1 = 0.001, 𝛼 = 0.001, 𝛾 = 1.2 ,   𝜆 = 0.02               
     𝐸2 = 0.005; 𝜃 = 0.001, ℎ = 0.03, 𝐸3 = 0.02, 𝜏 = 15.0, 𝑛 = 5000                          

  (55) 

It is noted that for all  𝜏 ≥ 0  and the data provided the solution of system (2) has globally 

asymptotically stable   𝐹1  as shown in Figure 1. 

     

Figure 1. The system's (2) trajectory based on the data provided in equation (55). 

 

(A) The system's  (2) time series   gradually converge to 𝐹1. (B) 3D phase diagram representing 

the globally asymptotic stability of point 𝐹1. 

It is observed that the same data from equation (55), with 𝑐 = 0.02   and ℎ = 0.001  system (2) 

exhibits global asymptotic stability for𝐹2  as depicted in Figure 2. 

 
 

Figure 2. The system's (2) trajectory based on the data provided in equation (55 ) with 𝑐 = 0.02   and ℎ = 0.001  

.(A) The system's (2) time series   gradually converge to 𝐹2. (B) 3D phase diagram representing the globally 

asymptotic stability of point 𝐹2. 

Here, we discuss the impact of time delays on the behavior of the system (2) near the 𝐹3 point 

It is noticed that the conditions of theorem 2 are met for the parameters in the data of equation 

(55) with = 0.02 . it is observed that when  𝜏 = 15  the  𝐹3  point is asymptotic stable as shown 

in Fig. (3) while for  𝜏 = τ0 = 45  a Hopf bifurcation occurs at   𝐹3 as shown in Figure 4. On 

other  hand  for   τ =  60 > τ0 =  45 increasing period as 𝜏 increases as shown in Figure 5. 

 



IHJPAS. 2025,38(2) 

373 
 

 

 
 

Figure 3.The system's (2) trajectory based on the data provided in equation (55 ) with 𝑐 = 0.02  and 𝜏 = 15   

(A)The system's (2) time series   gradually converge to 𝐹3. (B) 3D phase diagram representing the globally 

asymptotic stability of point. 

 
 

Figure 4.The system's (2) trajectory based on the data provided in equation(55 ) with 𝑐 = 0.02 𝑎𝑛𝑑 𝜏 = τ0 = 45  

.(A) a periodic solution's existence near 𝐹3    near 𝐹3 (B) 3D periodic solution. 

  

Figure 5. The system's (2) trajectory based on the data provided in equation(55) with 𝑐 = 0.02 and τ =  60 >

τ0 =  45 .(A) a periodic solution's existence near 𝐹3  (B) 3D periodic solution. 

 

4. Conclusions  

     In this study, a delayed predator-prey model with harvests is proposed. Our aim is to 

understand how the stability of the model is affected by the latent time of predation and the 

gestation period of the predators. All properties of the solution, such as positivity and 

boundedness, were investigated. It was found that the system (2) can have four equilibrium 

points. The stability analysis has shown that the discrete time delay has no effect on the stability 



IHJPAS. 2025,38(2) 

374 
 

 

of the axial equilibrium point and the planer equilibrium point, so it is still locally asymptotically 

stable for all 𝜏≥0. It has been demonstrated that the coexistence equilibrium point is 

characterized by asymptotic stability when the delay does not approach a critical value of ,τ-0..  

However, a Hopf bifurcation occurs when  =,τ-0. , making it an unstable point, and the solution 

approaches asymptotically to periodic dynamics for τ>,τ-0.. Additionally, the center manifold 

theorem was applied to investigate the stability and direction of bifurcating periodic dynamics. 

Finally, the Matlab program is employed for a numerical investigation of the system's global 

dynamics. For the given data in equation (55), ,-1. is globally asymptotically stable. When the 

same data from equation (55) is used with  =0.02   and ℎ=0.001  , it is observed that ,𝐹-2 . is 

globally asymptotically stable. For the data in equation (55) with =0.02 ,  it is noted that ,𝐹-3. 

is asymptotically stable when 𝜏=15. However, when ,=τ-0.=45, a Hopf bifurcation occurs. 

 

Acknowledgements  

     We would like to express our gratitude other referees for their valuable comments and 

suggestions that led to a truly significant  improvement of the paper. 

 

Conflict of Interest 

     The authors declare that there are no competing interests regarding the publication of this 

paper. 

 

Funding 

     This work is  not supported  by  any the Foundation. 

 

References 

1. Lotka AJ. Elements of Physical Biology. Baltimore: Williams & Wilkins; 1925.   

2. Volterra V. Variazioni e fluttuazioni del numero d'individui in specie animali conviventi. Rome: 

Società Anonima Tipografica "Leonardo da Vinci"; 1927. 

3. Anderson RM, May RM. The invasion, persistence, and spread of infectious diseases within animal 

and plant communities. Philos Trans R Soc Lond B Biol Sci. 1986;314(1167):533-

570. https://doi.org/10.1098/rstb.1986.0072 

4. Al-Momen SM, Naji RK. The Dynamics of Modified Leslie-Gower Predator-Prey Model Under 

the Influence of Nonlinear Harvesting and Fear Effect. Iraqi J Sci. 2022;63(1):259–

282. https://doi.org/10.24996/ijs.2022.63.1.27 

5. Dehingia K, Mohsen AA, Alharbi SA, Alsemiry RD, Rezapour S. Dynamical Behavior of a 

Fractional Order Model for Within-Host SARS-CoV-2. Mathematics. 

2022;10(13):2344. https://doi.org/10.3390/math10132344 

6. Hussien RM, Naji RK. The Dynamics of a Delayed Ecological Model with Predator Refuge and 

Cannibalism. Commun Math Biol Neurosci. 2023;2023:7988. https://doi.org/10.28919/cmbn/7988 

7. Ali NF. Modeling and Stability of Prey-Predator System Involving Infectious Disease in Each 

Population with Harvesting of the Prey. Al-Nahrain J Sci. 2015;18(4):144-

152. https://doi.org/10.22401/JNUS.18.4.20 

8. Hassan K, Mustafa A, Hama M. An Eco-Epidemiological Model Incorporating Harvesting Factors. 

Symmetry. 2021;13(11):2179. https://doi.org/10.3390/sym13112179 

9. Ali NF, Aaid AA. On the Dynamics of Prey-Predator Model Involving Treatment and Infections 

Disease in Prey Population. Iraqi J Sci. 2015;56(3C):2654-

2673. https://doi.org/10.24996/ijs.2015.56.3C.27 

10. Johri A, Trivedi N, Sisodiya A, Sing B, Jain S. Study of a prey-predator model with diseased prey. 

Int J Contemp Math Sci. 2012;7(10):489-498. 

https://doi.org/10.1098/rstb.1986.0072
https://doi.org/10.24996/ijs.2022.63.1.27
https://doi.org/10.3390/math10132344
https://doi.org/10.28919/cmbn/7988
https://doi.org/10.22401/JNUS.18.4.20
https://doi.org/10.3390/sym13112179
https://doi.org/10.24996/ijs.2015.56.3C.27


IHJPAS. 2025,38(2) 

375 
 

 

11. Sharma S, Samanta GP. Analysis of a Two Prey One Predator System with Disease in the First Prey 

Population. Int J Dyn Control. 2015;3(3):210–224. https://doi.org/10.1007/s40435-014-0107-4 

12. Ahmed LS, AL-Husseiny HF. Dynamical Behavior of an Eco-Epidemiological Model Involving 

Disease in Predator and Stage Structure in Prey. Iraqi J Sci. 2019;60(8):1766–

1782. https://doi.org/10.24996/ijs.2019.60.8.14 

13. Hugo A, Simanjilo E. Analysis of an eco-epidemiological model under optimal control measures 

for infected prey. Appl Appl Math. 2019;14(1):117-

138. https://digitalcommons.pvamu.edu/aam/vol14/iss1/8 

14. YP, Ma MJ, Zuo P, Liang X. Analysis of an eco-epidemiological model with disease in the predator. 

Appl Mech Mater. 2014;536:861-864. https://doi.org/10.4028/www.scientific.net/AMM.536-

537.861 

15. Abdul Satar H, Naji RK. Stability and bifurcation of a prey-predator-scavenger model in the 

existence of toxicant and harvesting. Int J Math Math Sci. 

2019;2019:1540015. https://doi.org/10.1155/2019/1540015 

16. Kant S, Kumar V. Stability analysis of predator–prey system with migrating prey and disease 

infection in both species. Appl Math Model. 2017;42:509-

539. https://doi.org/10.1016/j.apm.2017.07.017 

17. Naji RK, Ali NF. Modeling and Stability of Lotka-Volterra Prey-Predator System Involving 

Infectious Disease in Each Population. Iraqi J Sci. 2014;55(2):491-505. 

18. Song X, Chen L. Optimal harvesting and stability for a two-species competitive system with stage 

structure. Math Biosci. 2001;170(2):173-186. https://doi.org/10.1016/S0025-5564(01)00060-0 

19. Das A, Pal M. Theoretical analysis of an imprecise prey-predator model with harvesting and optimal 

control. J Optim. 2019;2019:1-12. 

20. Kar TK. Modelling and analysis of a harvested prey–predator system incorporating a prey refuge. J 

Comput Appl Math. 2006;185(1):19-33. 

21. Hale JH. Ordinary Differential Equations. New York: Wiley-Interscience; 1969. 

22. Hassard BD, Kazarinoff ND, Wan Y-H. Theory and Applications of Hopf Bifurcation. Cambridge: 

Cambridge University Press; 1981. 

23. Al-Jubouri KQ, Naji RK. Delay in eco-epidemiological prey-predator model with predation fear 

and hunting cooperation. Commun Math Biol Neurosci. 

2023;2023:89. https://doi.org/10.28919/cmbn/8081 

24. Pal AK, Bhattacharyya A, Pal S. Study of delay induced eco-epidemiological model incorporating 

a prey refuge. Filomat. 2022;36(2):557-578. https://doi.org/10.2298/FIL2202557P 

25. Tankam I, Tchinda Mouofo P, Mendy A, Lam M, Tewa JJ, Bowong S. Local Bifurcations and 

Optimal Theory in a Delayed Predator–Prey Model with Threshold Prey Harvesting. Int J Bifurcat 

Chaos. 2015;25(7):1540015. https://doi.org/10.1142/S021812741540015X 

26. Samanta S, Tiwari PK, Alzahrani AK, Alshomrani AS. Chaos in a nonautonomous eco-

epidemiological model with delay. Appl Math Model. 2020;79:865-

880. https://doi.org/10.1016/j.apm.2019.10.040 

27. Zhang X, Liu Z. Hopf bifurcation analysis in a predator-prey model with predator-age structure and 

predator-prey reaction time delay. Appl Math Model. 2021;91:530-

548. https://doi.org/10.1016/j.apm.2020.10.040 

28. Kamel Naji R, Abdullah Ibrahim H. Chaos in a harvested prey-predator model with infectious 

disease in the prey. J Al-Qadisiyah Comput Sci Math. 2011;3(2):1-21. 

29. Hale JH. Ordinary Differential Equations. New York: Wiley-Interscience; 1969. 

30. Hassard BD, Kazarinoff ND, Wan Y-H. Theory and Applications of Hopf Bifurcation. Cambridge: 

Cambridge University Press; 1981. 

 

 

 

https://doi.org/10.1007/s40435-014-0107-4
https://doi.org/10.24996/ijs.2019.60.8.14
https://digitalcommons.pvamu.edu/aam/vol14/iss1/8
https://doi.org/10.4028/www.scientific.net/AMM.536-537.861
https://doi.org/10.4028/www.scientific.net/AMM.536-537.861
https://doi.org/10.1155/2019/1540015
https://doi.org/10.1016/j.apm.2017.07.017
https://doi.org/10.1016/S0025-5564(01)00060-0
https://doi.org/10.28919/cmbn/8081
https://doi.org/10.2298/FIL2202557P
https://doi.org/10.1142/S021812741540015X
https://doi.org/10.1016/j.apm.2019.10.040
https://doi.org/10.1016/j.apm.2020.10.040