445 

© 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 

 

 
Influence of Infection Delay on the Covid-19 Pandemic with Vaccination 

Control: Modeling and Simulation 

 

Rami Raad Saadi 1*  and Hassan F. Al-Husseiny 2   
1Department of Networks and cyber security, College of Engineering, AL-Iraqia University, Baghdad, Iraq. 

1,2Department of Mathematics, College of Science, University of Baghdad, Baghdad, Iraq. 

*Corresponding Author.  

 

 Received:1 July 2023                       Accepted:4 September 2023                   Published:20 January 2025  

 doi.org/10.30526/38.1.3638 

 
 Abstract   

     In this work, we formulate a mathematical model of the killer COVID-19 pandemic with time 

delay and some governmental measures that include vaccine subsidies to understand the dynamic 

behavior of COVID-19. For the dynamic study, a new model, 𝑆𝑉1𝑉2𝐼𝑅 was used purposed in which 

infectious individuals were divided into five sub compartments. Our aim is to construct a more 

reliable and realistic model for a complete mathematical and computational analysis and design of 

different control strategies for the proposed deterministic model. We first obtain the basic 

reproduction number for the model is computed using a next-generation technique to predict the 

future dynamics of the pandemic. The local stability of the model was also investigated at each 

equilibrium point. The findings show that the time delay can produce a Hopf bifurcation for a 

𝑆𝑉1𝑉2𝐼𝑅 model. The obtained numerical results are discussed and predict through graphs. 

Keywords: Time delayed (T.D.), Equilibrium points (E.Ps.), COVID-19 Pandemic, Stability, 

Hopf Bifurcation (H.B.), Numerical simulation (N.S.). 

 

1. Introduction 

       Mathematical models of infectious disease transmission are increasingly being used to guide 

public health policy also they are used characterize the complex interactions, and enable 

information from diverse sources to give a clear understanding of the behavior of these diseases. 

Infectious disease epidemiology models are inherently multidisciplinary because the transmission 

of infection within a population is affected not just by the biologicl characteristics of the infectious 

agent and its host but also by the patterns of contact between hosts and the environment. There are 

many examples of the spread and control of epidemics, such as of examples Ferguson et al. (1) 

where studied the transmission intensity and impact of control policies on the foot and mouth 

epidemic in Great Britain. Donnelly et al. (2) suggested the epidemiological and genetic analysis 

of severe acute respiratory syndrome. Cauchemez et al. (3) studied the Middle East respiratory 

https://creativecommons.org/licenses/by/4.0/
https://creativecommons.org/licenses/by/4.0/
https://doi.org/10.30526/38.1.3501
https://orcid.org/0000-0001-7892-237X
mailto:rami.raad1103a@sc.uobaghdad.edu.iq
https://orcid.org/0009-0007-7883-5875
mailto:hassan.fadhil.r@sc.uobaghdad.edu.iq%20.


IHJPAS. 2025, 38 (1) 
 

446 
 

syndrome coronavirus: quantification of the extent of the epidemic, surveillance biases, and 

transmissibility.  

Mohsen et al. (4) studied the global stability of COVID-19 model involving the quarantine strategy 

and media coverage effects. AL-Husseiny et al. (5) have discussed the effect of individuals 

asymptomatic (Carrier) on the dynamical behavior of a COVID-19 virus. Hattaf et al. (6) they 

suggested modeling the dynamics of COVID-19 with carrier effect and environmental 

contamination. Abdulkadhim and Al-Husseiny (7) studied the global stability and bifurcation of a 

COVID-19 virus modeling with possible loss of the immunity.  

The dynamics of populations are significantly influenced by time delays. The dynamics of state 

variables in many real-world processes, notably in many biological phenomena, depend on the 

phenomenon's history, or on the state variables' former values, in addition to the phenomenon's 

current state. Time delays may have an impact on the dynamics of infectious diseases, as shown 

in Zuo et al. (8) formulated the relationship between media coverage and an epidemic's recruitment 

and spread. Aekabut et al. (9)investigated a delayed SEIR epidemic model in which the diseased 

and latent phases are contagious. Rasha M. Yaseen et al. (10) studied the Stability and Hopf 

bifurcation of an epidemiological model with the effect of delaying the awareness programms and 

vaccination: analysis and simulation. Zhe Yin et al. (11) investigated how an age-structured SEIRS 

model was affected by time delays. Mohsen et al. (12) investigated the dynamics of a curfew 

technique in a model of the coronavirus pandemic epidemic. Zizhen et al. (13) suggested SVIRS 

epidemic model includes a number of delays along with incidence and treatment rates for Holling 

type II. Dehingia et al. (14) investigated the dynamic behavior of a SARS-CoV-2 within-host 

fractional order model. Shurowq et al. (15) study of the COVID-19 epidemic and the bifurcation 

analysis of a mathematical model for vaccination. Naji and Hussien (16) proposed and analyzed 

the epidemic model type of SEIR with nonlinear incidence and treatment rates and also used time 

delays owing to the incubation period. Ahmed et al. (17) discussed a mathematical model for the 

dynamics of COVID-19 pandemic involving infected immigrants. Naji and Mohsen [18] studied 

the stability analysise with the bifurcation of an SVIRE epidemic model involving immigrants. 

Ahmed and AL- Husseiny (19) studied the dynamical behavior of an eco-epidemiological model 

involving disease in predators and stage structure in prey. Mohsen and Hattaf (20) studied the 

dynamics of a generalized fractional epidemic model of COVID-19 with carrier effect. 

This study presents and evaluates a mathematical model that represents the dynamics of the 

COVID-19 pandemic's delayed infection and includes two stages of vaccination. This paper is 

organized as follows: Section 2 illustrates the innovative coronavirus mathematical modeling and 

the two steps of immunization with delayed infection. Local stability and H.B. are discussed in 

Section 3. In Section 4, a numerical simulation is utilized to analyze the effects of altering every 

system parameter. 

 

 

 

 

 



IHJPAS. 2025, 38 (1) 
 

447 
 

2. Mathematical model: see in(21) 

     The system of the time delayed is: 

 

      
𝑑𝑆

𝑑𝑡
= 𝛬 −

𝛼𝑆

1+𝑛𝑉1
− 𝛽1𝑆(𝑡 − 𝜏)𝐼(𝑡 − 𝜏) − 𝜇𝑆                                                     

𝑑𝑉1

𝑑𝑡
=

𝛼𝑆

1+𝑛𝑉1
− 𝛾𝑉1 − 𝛽2𝑉1𝐼 − 𝜇𝑉1                                                                  

𝑑𝑉2

𝑑𝑡
= 𝛾𝑉1 − 𝛽3𝑉2𝐼 − 𝜇𝑉2                                                                                

𝑑𝐼

𝑑𝑡
= 𝛽1𝑆(𝑡 − 𝜏)𝐼(𝑡 − 𝜏) + 𝛽2𝑉1𝐼 + 𝛽3𝑉2𝐼 − (𝜇 + 𝜇1)𝐼 − 𝜃𝐼                

𝑑𝑅

𝑑𝑡
= 𝜃𝐼 − 𝜇𝑅                                                                                                    

 

                                                                                                                                                       (1) 

Now, see(21) we get; 

(𝑡) =
𝜃𝐼

𝜇
,                                                                                                                                            (2) 

As a result, the system below will be studied rather than system (1); 

 
𝑑𝑆

𝑑𝑡
= Λ −

𝛼𝑆

1+𝑛𝑉1
− 𝛽1𝑆(𝑡 − 𝜏)𝐼(𝑡 − 𝜏) − 𝜇𝑆                                                         

𝑑𝑉1

𝑑𝑡
=

𝛼𝑆

1+𝑛𝑉1
− 𝛾𝑉1 − 𝛽2𝑉1𝐼 − 𝜇𝑉1                                                                           

𝑑𝑉2

𝑑𝑡
= 𝛾𝑉1 − 𝛽3𝑉2𝐼 − 𝜇𝑉2                                                                                         

 
𝑑𝐼

𝑑𝑡
= 𝛽1𝑆(𝑡 − 𝜏)𝐼(𝑡 − 𝜏) + 𝛽2𝑉1𝐼 + 𝛽3𝑉2𝐼 − (𝜇 + 𝜇1)𝐼 − 𝜃𝐼                        

                               (3) 

 

                       

3. Local stability analysis (L.S.A.) and hopf bifurcation (H.B.) 

     In this section, the L.S. and H.B. of system (3) are studied. The position and quantity of 

equilibrium points are known and don't alter with time delays. 

 Accordingly, from (21)  system (3) have six E.Ps., say  

 

 

𝐸0 = (𝑆,̅ 0,0,0) 𝑤ℎ𝑒𝑛 𝛼 = 0 ,   𝐸1 = (𝑆̿, �̿�1, 0,0) 𝑤ℎ𝑒𝑛 𝛾 = 0,

𝐸2 = (�̂�, 0,0, 𝐼) 𝑤ℎ𝑒𝑛 𝛼 = 0 , 𝐸3 = (�̆�, �̆�1, �̆�2, 0),                     

𝐸4 = (�̃̃�, �̃̃�1, 0, 𝐼)  𝑤ℎ𝑒𝑛 𝛾 = 0 𝑎𝑛𝑑  𝐸5 = (𝑆∗, 𝑉1
∗, 𝑉2

∗, 𝐼∗).       
 
 

 

It is well known that, the Jacobian matrix(J.M) of system (3) at any E.Ps. 𝐸 = (𝑆, 𝑉1, 𝑉2, 𝐼) is 

𝐽(𝐸) = (𝑎𝑖𝑗)4×4    ; 𝑖, 𝑗 = 1,2,3,4. 

where 

𝑎11 = −[
𝛼

1+𝑛𝑉1
+ 𝛽1𝐼𝑒

−𝜆𝜏 + 𝜇]   , 𝑎12 =
𝑛𝛼𝑆

(1+𝑛𝑉1)2
 , 𝑎14 = −𝛽1𝑆𝑒

−𝜆𝜏 , 𝑎21 =
𝛼

1+𝑛𝑉1
 ,                     

𝑎22 = −[
𝑛𝛼𝑆

(1+𝑛𝑉1)2
+ 𝛾 + 𝛽2𝐼 + 𝜇]    , 𝑎24 = −𝛽2𝑉1 , 𝑎32 = 𝛾     , 𝑎33 = −[𝛽3𝐼 + 𝜇],                      

𝑎34 = −𝛽3𝑉2 , 𝑎41 = 𝛽1𝐼𝑒
−𝜆𝜏 , 𝑎42 = 𝛽2𝐼 ,   𝑎43 = 𝛽3𝐼, 𝑎13 = 𝑎23 = 𝑎31 = 0,                          

𝑎44 = 𝛽1𝑆𝑒
−𝜆𝜏 + 𝛽2𝑉1 + 𝛽3𝑉2 − (𝜇 + 𝜇1 + 𝜃)                                                                                   

    

     

                                                                                                                                                       (4)                                                

….    

while its associated characteristic equation(C.E.) takes the form 



IHJPAS. 2025, 38 (1) 
 

448 
 

 𝑃(𝜆) + 𝑄(𝜆)𝑒−𝜆𝜏 = 0                                                                                                       (5) 

here 𝑃(𝜆) and 𝑄(𝜆) are polynomials of 𝜆. Accordingly the L.S. properties of system (3) at all 

feasible E.Ps. are determined by the roots of the above equation for all 𝜏 ≥ 0. 

 

• For the (FEP) when  𝛼 = 0, Equation (4) reduces to  

𝐽(𝐸0) =

[
 
 
 
−𝜇 0 0 −𝛽1𝑆̅𝑒

−𝜆𝜏

0 −(𝛾 + 𝜇) 0 0
0 𝛾 −𝜇 0

0 0     0     𝛽1𝑆̅𝑒
−𝜆𝜏 − (𝜇 + 𝜇1 + 𝜃)]

 
 
 

                                                   (6) 

The (C.E.) of 𝐽(𝐸0) is: 

(𝛽1𝑆̅𝑒
−𝜆𝜏 − (𝜇 + 𝜇1 + 𝜃) − 𝜆)(−𝜇 − 𝜆)(−(𝛾 + 𝜇) − 𝜆)(−𝜇 − 𝜆) = 0                                   (6a)     

Equation (6a) represents the eigenvalues of 𝐽(𝐸0) and has 4-roots: 

𝜆1 = 𝛽1𝑆̅𝑒
−𝜆𝜏 − (𝜇 + 𝜇1 + 𝜃)

𝜆2 = −𝜇                                      

𝜆3 = −(𝛾 + 𝜇)

𝜆4 = −𝜇            
                          

}
 

 
                                                                                                 (6b) 

Now, for 𝜏 = 0 we get the eigenvalues will be negative and FEP is locally asymptotically stable 

(L.A.S.) if  

𝛽1𝑆̅ < 𝜇 + 𝜇1 + 𝜃                                                                                                                        (7) 

Now, for  𝜏 > 0  the Equation (6b), has two wholly imaginary roots, namely= ±𝑖𝜔 ( 𝜔 > 0). 
By substituting 𝜆 = ±𝑖𝜔  in Equation (6b) we get: 

 𝛽1𝑆̅(𝑐𝑜𝑠𝜔𝜏 − 𝑖𝑠𝑖𝑛𝜔𝜏) = 𝜇 + 𝜇1 + 𝜃 + 𝑖𝜔 

So, separating the real and imagined parts yields 

𝛽1𝑆̅𝑐𝑜𝑠𝜔𝜏 = 𝜇 + 𝜇1 + 𝜃

𝛽1𝑆̅𝑠𝑖𝑛𝜔𝜏 = −𝜔              
}                                                                                                             (8) 

Squaring each equation and then adding them, we get that 

 𝜔 = ∓√𝛽1
2(𝑆̅)2 − (𝜇 + 𝜇1 + 𝜃)2 

Note that, under the condition (7),𝜔(𝜏) with 𝜏 > 0 cannot be real, which contradicts with the 

assumption. Therefore, the C.E. (6a) can’t have purely imaginary root, and FEP is L.A.S. for all 

𝜏 ≥ 0 if the condition (7) hold. 

• For the (SEP) when  𝛾 = 0, Equation (4) reduces to  

𝐽(𝐸1) =

[
 
 
 
 
 
 
 −(

𝛼

1 + 𝑛�̿�1
+ 𝜇)

𝑛𝛼𝑆̿

(1 + 𝑛�̿�1)
2 0 −𝛽1𝑆̿𝑒

−𝜆𝜏

𝛼

1 + 𝑛�̿�1
−(

𝑛𝛼𝑆̿

(1 + 𝑛�̿�1)
2 + 𝜇) 0 −𝛽2�̿�1

0 0 −𝜇 0

0 0     0      𝛽1𝑆̿𝑒
−𝜆𝜏 + 𝛽2�̿�1 − (𝜇 + 𝜇1 + 𝜃)]

 
 
 
 
 
 
 

 

                                                                

 (9) 

The C.E. of 𝐽(𝐸1) is 

[𝜆2 + 𝐴1𝜆 + 𝐴2][−𝜇 − 𝜆][𝛽1𝑆̿𝑒
−𝜆𝜏 + 𝛽2�̿�1 − (𝜇 + 𝜇1 + 𝜃) − 𝜆] = 0                                    (10a) 

Where  



IHJPAS. 2025, 38 (1) 
 

449 
 

 𝐴1 =
𝛼

1+𝑛�̿�1
+

𝑛𝛼�̿�

(1+𝑛𝑉1)
2 + 2𝜇 

 𝐴2 = 𝜇 (
𝑛𝛼�̿�

(1+𝑛�̿�1)
2 +

𝛼

1+𝑛�̿�1
+ 𝜇) 

 The Equation (10a) represents the eigenvalues of 𝐽(𝐸1) and has 4-roots: 

𝜆1,2 = −
𝐴1

2
∓
1

2
√𝐴1

2 − 4𝐴2                   

𝜆3 = −𝜇                                                     

𝜆4 = 𝛽1𝑆̿𝑒
−𝜆𝜏 + 𝛽2�̿�1 − (𝜇 + 𝜇1 + 𝜃)

}                                                                                   (10b) 

Now, for 𝜏 = 0 we get all the above eigenvalues will be negative and the SEP is L.A.S. if   

𝛽1𝑆̿ + 𝛽2�̿�1 < 𝜇 + 𝜇1 + 𝜃                                                                                                             (11) 

Now, for  𝜏 > 0  the Equation (10b), has two wholly imaginary roots, namely= ±𝑖𝜔 ( 𝜔 > 0). 
By substituting 𝜆 = ±𝑖𝜔  in Equation (10b)  we get : 

 𝛽1𝑆̿(𝑐𝑜𝑠𝜔𝜏 − 𝑖𝑠𝑖𝑛𝜔𝜏) = 𝜇 + 𝜇1 + 𝜃 − 𝛽2�̿�1 + 𝑖𝜔 

So, separating the real and imagined parts yields 

𝛽1𝑆̿𝑐𝑜𝑠𝜔𝜏 = 𝜇 + 𝜇1 + 𝜃 − 𝛽2�̿�1

𝛽1𝑆̿𝑠𝑖𝑛𝜔𝜏 = −𝜔                            
}                                                                                              (12) 

Squaring each equation and then adding them, we get that 

 𝜔 = ∓√𝛽1
2(𝑆̅̅)

2
− (𝜇 + 𝜇1 + 𝜃 − 𝛽2�̿�1)

2
 

Note that, under the condition (11),𝜔(𝜏) with 𝜏 > 0 cannot be real, which contradicts with the 

assumption. Therefore, the C.E. (10a) can’t have purely imaginary root, and SEP is L.A.S. for all 

𝜏 ≥ 0 if the condition (11) hold.  

 

For the (TEP) when  𝛼 = 0, Equation (4) reduces to  

𝐽(𝐸2) =

[
 
 
 
 
−(𝑅1𝑒

−𝜆𝜏 + 𝑅2) 0 0 −𝑅3𝑒
−𝜆𝜏

0 −(𝛾 + 𝛽2𝐼 + 𝜇) 0 0

0 𝛾 −(𝛽3𝐼 + 𝜇) 0

𝑅1𝑒
−𝜆𝜏 𝛽2𝐼     𝛽3𝐼      𝑅3𝑒

−𝜆𝜏 + 𝑅4]
 
 
 
 

                                                                                                                   

                                                                                                                                                                                             

                                                                                                                                                     (13) 

Where, 𝑅1 = 𝛽1𝐼   , 𝑅2 = 𝜇  , 𝑅3 = 𝛽1�̂� , 𝑅4 = −(𝜇 + 𝜇1 + 𝜃).  
The C.E. of 𝐽(𝐸2) is given by 

[𝜆2 + 𝐵1𝜆 + 𝐵2 + (𝐵3𝜆 + 𝐵4)𝑒
−𝜆𝜏][−(𝛾 + 𝛽2𝐼 + 𝜇) − 𝜆]

                                                                                  [−(𝛽3𝐼 + 𝜇) − 𝜆] = 0
                                         (14)                                  

With  

 
𝐵1 = 𝑅2 − 𝑅4 ,   𝐵2 = −𝑅2𝑅4,                                 
𝐵3 = 𝑅1 − 𝑅3  ,   𝐵4 = −𝑅2𝑅3 − 𝑅1𝑅4.               

 

Now, when = 0 , Equation (14) is; 

[𝜆2 + (𝐵1 + 𝐵3)𝜆 + 𝐵2 + 𝐵4][−(𝛾 + 𝛽2𝐼 + 𝜇) − 𝜆][−(𝛽3𝐼 + 𝜇) − 𝜆] = 0 

                                                                                                                                                     (15) 

So, either 

[−(𝛾 + 𝛽2𝐼 + 𝜇) − 𝜆][−(𝛽3𝐼 + 𝜇) − 𝜆] = 0                                                                           (16a) 

Or,  



IHJPAS. 2025, 38 (1) 
 

450 
 

[𝜆2 + (𝐵1 + 𝐵3)𝜆 + 𝐵2 + 𝐵4] = 0                                                                                           (16b) 

From Equation (16a) we obtain that  

 
𝜆2 = −(𝛾 + 𝛽2𝐼 + 𝜇) < 0

𝜆3 = −(𝛽3𝐼 + 𝜇) < 0         
 

Which is always negative eigenvalue. 

Now, it is easy to verify that 𝐵1 + 𝐵3 > 0 and 𝐵2 + 𝐵4 > 0 under the following sufficient 

conditions  

𝛽1�̂� < 𝛽1𝐼 + 2𝜇 + 𝜇1 + 𝜃,                                                                                                          (17a) 

𝜇𝛽1�̂� < (𝛽1𝐼 + 𝜇)(𝜇 + 𝜇1 + 𝜃).                                                                                                (17b) 

We see that all roots of Equation (16b), which represent the eigenvalues of (13), have negative 

real parts. Consequently, under the conditions (17a-17b), TEP is L.A.S. for system (3) when 𝜏 =

0. 

Now, for 𝜏 > 0, then  

Either, Equation (16a) Which is always negative eigenvalue. 

Or,  

[𝜆2 + 𝐵1𝜆 + 𝐵2 + (𝐵3𝜆 + 𝐵4)𝑒
−𝜆𝜏] = 0                                                                                    (18) 

From (18) has two wholly imaginary roots, namely= ±𝑖𝜔 ( 𝜔 > 0). 

By substituting 𝜆 = ±𝑖𝜔  in Equation (18) and separating the real and imagined parts yields

 
𝐵4𝑐𝑜𝑠𝜔𝜏 + 𝐵3𝜔𝑠𝑖𝑛𝜔𝜏 = 𝜔

2 − 𝐵2
𝐵3𝜔𝑐𝑜𝑠𝜔𝜏 − 𝐵4𝑠𝑖𝑛𝜔𝜏 = −𝐵1𝜔    

                                                                                             (19) 

Through squaring and adding both equations in Equation (19), we get  

𝜔4 + ℎ1𝜔
2 + ℎ2 = 0                                                                                                                  (20) 

Further, by letting ℏ = 𝜔2 in Equation (20) 

𝑔(ℏ) = ℏ2 + 𝑔1ℏ + 𝑔2 = 0                                                                                                      (21a) 

Where 

 
𝑔1 = 𝑅2

2 + 𝑅4
2 − (𝑅1 − 𝑅3)

2,                              

𝑔2 = 𝑅2
2𝑅4

2−(𝑅2𝑅3 + 𝑅1𝑅4)
2.                             

  

Straightforward computation shows that due to the following condition 

𝑅2𝑅4 < 𝑅2𝑅3 + 𝑅1𝑅4                                                                                                               (21b) 

we obtain 𝑔2 < 0. So, there exists a unique positive root in accordance with Descartes' rule of 

signs ℏ0 = 𝜔0
2 providing Equation (21a). Which is, Equation (20) has a positive𝜔0. As a result, 

Equation (18) has at least two roots ±𝜔0 that are purely imaginary and which correspond to the 

time delay 𝜏. Additionally, by substituting 𝜔0 in (19) and solving the system for 𝜏, yields : 

𝜏𝑘 =
1

𝜔0
𝑐𝑜𝑠−1

(𝐵4−𝐵1𝐵3)𝜔0
2−𝐵2𝐵4

𝐵4
2+𝐵3

2𝜔0
2 +

2𝑘𝜋

𝜔0
                                                                                        (22) 

Where, k = 0,1,2,…. Hence we get the corresponding 𝜏𝑘 > 0 for which system (3) has two wholly 

imaginary roots ±𝜔0. 

Let 𝜆(𝜏) = 𝜇(𝜏) + 𝑖𝜔(𝜏) be a root of Equation (18) near 𝜏 = 𝜏𝑘 with  𝜇(𝜏𝑘) = 0 and 𝜔(𝜏𝑘) =
𝜔0. Then comes the next theorem: 

 

 

Theorem1. The C.E. (18), roots satisfy the following transversality requirement; 

[
𝑑(𝑅𝑒𝜆(𝜏))

𝑑𝜏
]
𝜏=𝜏𝑘

> 0                                                                                                                      (23a) 

Provided that  



IHJPAS. 2025, 38 (1) 
 

451 
 

𝜔0
2 > 𝐵2                                                                                                                                     (23b) 

 

Proof. By substituting 𝜆(𝜏) in (18) and differentiating the resulting equation in 𝜏,we can get  

[2𝜆 + 𝐵1 + 𝐵3𝑒
−𝜆𝜏 − 𝜏(𝐵3𝜆 + 𝐵4)𝑒

−𝜆𝜏]
𝑑𝜆

𝑑𝜆
= 𝜆(𝐵3𝜆 + 𝐵4)𝑒

−𝜆𝜏                                              (24) 

Thus,  

(
𝑑𝜆

𝑑𝜏
)
−1

=
2𝜆+𝐵1

𝜆(𝐵3𝜆+𝐵4)𝑒−𝜆𝜏
+

𝐵3

𝜆(𝐵3𝜆+𝐵4)
−

𝜏

𝜆
                                                                                    (25) 

Since for = 𝜏0 , and 𝜆 = 𝑖𝜔0, we have got  

[
𝑑𝜆

𝑑𝜏
]
𝜏=𝜏0

−1

=
𝐵1+2𝑖𝜔0

𝐵1𝜔0
2+𝑖𝜔0(𝜔0

2−𝐵2)
+

𝐵3

𝐵3𝜔0
2+𝑖𝐵4𝜔0

−
𝑖𝜏0

𝜔0
                                                                          (26) 

 

Now, since  

𝑠𝑖𝑔𝑛 [
𝑑(𝑅𝑒𝜆)

𝑑𝜏
]
𝜏=𝜏0

= 𝑠𝑖𝑔𝑛 [𝑅𝑒 (
𝑑𝜆

𝑑𝜏
)
−1

]
𝜆=𝑖𝜔0

                                                                               (27) 

It is clear that: 

 

𝑅𝑒 [
𝐵1+2𝑖𝜔0

𝐵1𝜔0
2+𝑖𝜔0(𝜔0

2−𝐵2)
] =

𝐵1
2+2(𝜔0

2−𝐵2)

𝐵1
2𝜔0

2+(𝜔0
2−𝐵2)

2  ,                                     

𝑅𝑒 [
𝐵3

𝐵3𝜔0
2+𝑖𝐵4𝜔0

] =
𝐵3
2

𝐵3
2+𝐵4

2 ,                                                                    

𝑅𝑒 [
𝑖𝜏0

𝜔0
] = 𝑧𝑒𝑟𝑜 .                                                                            

 

Hence, we have  

 𝑅𝑒 [
𝑑𝜆

𝑑𝜏
]
𝜏=𝜏0

−1

=
𝐵1
2+2(𝜔0

2−𝐵2)

φ0
+

𝐵3
2

φ1
 

Where  

 
φ0 = 𝐵1

2𝜔0
2 + (𝜔0

2 − 𝐵2)
2 > 0 ,

φ1 = 𝐵3
2 + 𝐵4

2 > 0 .                     
 

 We obtain [
𝑑(𝑅𝑒𝜆(𝜏))

𝑑𝜏
]
𝜏=τ0

> 0    under (23b). 

This outcome demonstrates how the roots of C.E. (18), as 𝜏 passes through τ0, traverse the 

imaginary axis from left to right. As a result, the system (3) experiences an H.B. at 𝜏 = τ0 and 

loses its stability. 

 

For the (FOEP), Equation (5) reduces to  

𝐽(𝐸3) =

[
 
 
 
 
 
 
 − (

𝛼

1 + 𝑛�̆�1
+ 𝜇)

𝑛𝛼�̆�

(1 + 𝑛�̆�1)
2 0 −𝛽1�̆�𝑒

−𝜆𝜏

𝛼

1 + 𝑛�̆�1
−(

𝑛𝛼�̆�

(1 + 𝑛�̆�1)
2 + 𝛾 + 𝜇) 0 −𝛽2�̆�1

0 𝛾 −𝜇 −𝛽3�̆�2
0 0     0      𝛽1�̆�𝑒

−𝜆𝜏 + 𝛽2�̆�1 + 𝛽3�̆�2 − (𝜇 + 𝜇1 + 𝜃)]
 
 
 
 
 
 
 

 

                                                               

(28) 

The C.E. of 𝐽(𝐸3) is given by 

[𝜆2 + 𝐶1𝜆 + 𝐶2][−𝜇 − 𝜆][ 𝛽1�̆�𝑒
−𝜆𝜏 + 𝛽2�̆�1 + 𝛽3�̆�2 − (𝜇 + 𝜇1 + 𝜃) − 𝜆] = 0     

                                                                                                                                                   (29a) 



IHJPAS. 2025, 38 (1) 
 

452 
 

Where  

 𝐶1 =
𝛼

1+𝑛𝑉1
(1 +

𝑛�̆�

1+𝑛𝑉1
) + 𝛾 + 2𝜇 

 𝐶2 =
𝛼

1+𝑛𝑉1
(
𝑛𝜇�̆�

1+𝑛𝑉1
+ 𝛾 + 𝜇) + 𝜇(𝛾 + 𝜇) 

The eq. (29a) represents the eigenvalues of 𝐽(𝐸3) and has 4-roots:

 

𝜆1,2 = −
𝐶1

2
∓
1

2
√𝐶1

2 − 4𝐶2                                  

𝜆3 = −𝜇                                                                   

𝜆4 = 𝛽1�̆�𝑒
−𝜆𝜏 + 𝛽2�̆�1 + 𝛽3�̆�2 − (𝜇 + 𝜇1 + 𝜃)

}                                                          (29b) 

Now, for 𝜏 = 0 we get all the above eigenvalues will be negative and the FOEP is L.A.S. if the 

following condition  

𝛽1�̆� + 𝛽2�̆�1 + 𝛽3�̆�2 < 𝜇 + 𝜇1 + 𝜃                                                                                              (30) 

Now, for  𝜏 > 0  suppose that (29b), has two wholly imaginary roots, namely𝜆 = ±𝑖𝜔 ( 𝜔 > 0). 

By substituting 𝜆 = ±𝑖𝜔  in Equation (29b) we get: 

 𝛽1�̆�(𝑐𝑜𝑠𝜔𝜏 − 𝑖𝑠𝑖𝑛𝜔𝜏) = 𝜇 + 𝜇1 + 𝜃 − 𝛽2�̆�1 − 𝛽3�̆�2 + 𝑖𝜔  

So, separating the real and imagined parts yields 

𝛽1�̆�𝑐𝑜𝑠𝜔𝜏 = 𝜇 + 𝜇1 + 𝜃 − 𝛽2�̆�1 − 𝛽3�̆�2
𝛽1�̆�𝑠𝑖𝑛𝜔𝜏 = −𝜔                                          

}                                                                                  (31) 

Squaring each equation and then adding them, we get that 

 𝜔 = ∓√𝛽1
2(�̆�)

2
− (𝜇 + 𝜇1 + 𝜃 − 𝛽2�̆�1 − 𝛽3�̆�2)

2
 

Note that, under the condition (30),𝜔(𝜏) with 𝜏 > 0 cannot be real, which contradicts with the 

assumption. Therefore, the C.E. (29a) can’t have purely imaginary root, and 𝐸1 is L.A.S. for all   

𝜏 ≥ 0 if the condition (30) hold. 

 

For the (FIEP) when  𝛾 = 0, Equation (4) reduces to                    

         𝐽(𝐸4) = (𝑑𝑖𝑗)4×4   ;  𝑖, 𝑗 = 1,2,3,4 

Here 

𝑑11 = −[𝑅1 + 𝑅2𝑒
−𝜆𝜏] , 𝑑12 =

𝑛𝛼�̃̃�

(1 + 𝑛�̃̃�1)
2     , 𝑑14 = −𝑅3𝑒

−𝜆𝜏, 𝑑21 =
𝛼

1 + 𝑛�̃̃�1
 ,                                

  𝑑22 = −(
𝑛𝛼�̃̃�

(1 + 𝑛�̃̃�1)
2 + 𝛽2𝐼 + 𝜇 ) , 𝑑24 = −𝛽2�̃̃�1  , 𝑑33 = −(𝛽3𝐼 + 𝜇), 𝑑41 = 𝑅2𝑒

−𝜆𝜏,               

  
𝑑42 = 𝛽2𝐼 , 𝑑43 = 𝛽3𝐼  ,   𝑑44 = 𝑅3𝑒

−𝜆𝜏 + 𝑅4, 𝑑13 = 𝑑23 = 𝑑31 = 𝑑32 = 𝑑34 = 0 .                                                                                  

                                                                                            

 

                                                                                                                                                     (32) 

Where, 𝑅1 = 
𝛼

1+𝑛�̃̃�1
+ 𝜇  , 𝑅2 = 𝛽1𝐼  , 𝑅3 = 𝛽1�̃̃� , 𝑅4 = 𝛽2�̃̃�1 − (𝜇 + 𝜇1 + 𝜃).  

The C.E of 𝐽(𝐸4) is  

[−(𝛽3𝐼 + 𝜇) − 𝜆] [𝜆
3 + 𝐷1𝜆

2 + 𝐷2𝜆 + 𝐷3 + (𝐷4𝜆
2 + 𝐷5𝜆 + 𝐷6)𝑒

−𝜆𝜏] = 0                                        

                                                                                                                                                     (33) 

  With  



IHJPAS. 2025, 38 (1) 
 

453 
 

𝐷1 = 𝑅1 − 𝑑22 − 𝑅4 , 𝐷2 = −(𝑑22𝑅1 + 𝑅1𝑅4 + 𝑑12𝑑21 + 𝑑23𝑑32) + 𝑑22𝑅4 ,                                                                                                                          
𝐷3 = 𝑅1𝑅4𝑑22 + 𝑅4𝑑12𝑑21 − 𝑅1𝑑23𝑑32 ,  𝐷4 = 𝑅2 − 𝑅3  ,                                                                                                                                                            

    

𝐷5 = −(𝑅2𝑑22 + 𝑅1𝑅3 + 𝑅2𝑅4) + 𝑑22𝑅3 ,                                                                                                                      

 𝐷6 = 𝑅1𝑅3𝑑22 + 𝑅3𝑑12𝑑21 + 𝑅2𝑅4𝑑22 + 𝑅3𝑑21𝑑32 − (𝑅2𝑑12𝑑23 + 𝑅2𝑑12𝑑23 + 𝑅1𝑑23𝑑32 + 𝑅2𝑑232𝑑32).
                                                                        

 

Now, when = 0 , Equation (33) becomes  

[−(𝛽3𝐼 + 𝜇) − 𝜆] [𝜆
3 + (𝐷1 + 𝐷4)𝜆

2 + (𝐷2 +𝐷5)𝜆 + 𝐷3 + 𝐷6] = 0 

                                                                                                                                                   (34) 

So, either 

[−(𝛽3𝐼 + 𝜇) − 𝜆] = 0                                                                                                              (35a) 

or,  

 [𝜆3 + (𝐷1 + 𝐷4)𝜆
2 + (𝐷2 + 𝐷5)𝜆 + 𝐷3 + 𝐷6] = 0                                                      (35b) 

From Equation (35a) then    

 𝜆3 = −𝛽3𝐼 − 𝜇 < 0 

Which is always negative eigenvalue. 

Now, it is easy to verify that 𝐷1 + 𝐷4 > 0 and 𝐷3 + 𝐷6 > 0 under the following sufficient 

conditions  

𝛽1�̃̃� + 𝛽2�̃̃�1 < 𝜇 + 𝜇1 + 𝜃                                                                                                         (36a) 

𝑛𝛼2�̃̃�

(1+𝑛�̃̃�1)
3 < (

𝛼

1+𝑛�̃̃�1
+ 𝛽1𝐼 + 𝜇)(

𝑛𝛼�̃̃�

(1+𝑛�̃̃�1)
2 + 𝛽2𝐼 + 𝜇)                                                               (36b) 

We have  

 (𝐷1 + 𝐷4)(𝐷2 + 𝐷5) − (𝐷3 +𝐷6) > 0. 

Under the following sufficient conditions 

𝑛𝛽1𝛼�̃̃�

(1+𝑛�̃̃�1)
2 < (

𝑛𝛼�̃̃�

(1+𝑛�̃̃�1)
2 + 𝛽2𝐼 + 𝜇)𝛽2                                                                                        (36c) 

𝛽2�̃̃�1 (𝛽1�̃̃� + 𝛽2�̃̃�1) +
𝛽1𝛼�̃̃�

1+𝑛�̃̃�1
< 𝛽2�̃̃�1(𝜇 + 𝜇1 + 𝜃)                                                                   (36d) 

So, according to Routh-Hurwitz criterion, we see that all roots of Equation (35b), which represent 

the eigenvalues of (32), have negative real parts. Consequently, under the conditions (36a)-(36d), 

𝐸4 is L.A.S for system (3) when 𝜏 = 0. 

 

Now, for 𝜏 > 0, then  

Either, Equation (35a) which is always negative eigenvalue. 

or,  

[𝜆3 + 𝐷1𝜆
2 + 𝐷2𝜆 + 𝐷3 + (𝐷4𝜆

2 + 𝐷5𝜆 + 𝐷6)𝑒
−𝜆𝜏] = 0                                                         (37) 

From (37) has two wholly imaginary roots, namely= ±𝑖𝜔 ( 𝜔 > 0). 

By substituting 𝜆 = ±𝑖𝜔  in Equation (37) and separating the real and imagined parts yields 

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

2 −𝐷3
𝐷5𝜔𝑐𝑜𝑠𝜔𝜏 + (𝐷4𝜔

2 − 𝐷6)𝑠𝑖𝑛𝜔𝜏 = 𝜔
3 − 𝐷2𝜔

                                                                      (38) 

Through squaring and adding both equations in Equation (38), we get  

𝜔6 + ℎ1𝜔
4 + ℎ2𝜔

2 + ℎ3 = 0                                                                                                   (39) 

Further , by letting 𝜐 = 𝜔2 in Equation (39) 

ℎ(𝜐) = 𝜐3 + ℎ1𝜐
2 + ℎ2𝜐 + ℎ3 = 0                                                                                         (40a) 

Where 



IHJPAS. 2025, 38 (1) 
 

454 
 

 

ℎ1 = 𝐷1
2 − 2𝐷2 − 𝐷4

2,                                                  

ℎ2 = 2𝐷4𝐷6 + 𝐷2
2 − 𝐷5

2 − 2𝐷1𝐷3,                             

ℎ3 = 𝐷3
2 − 𝐷6

2.                                                                 

  

Straightforward computation shows that due to the following condition 

𝐷3 < 𝐷6                                                                                                                                     (40b) 

we obtain ℎ3 < 0. So, there exists a unique positive root in accordance with Descartes' rule of 

signs 𝜐0 = 𝜔0
2 providing Equation (40a). Which is, Equation (39) has a positive 𝜔0. As a result, 

Equation (37) has at least two roots ±𝜔0 that are purely imaginary and which correspond to the 

time delay 𝜏. Additionally, by substituting 𝜔0 in (38) and solving the system for 𝜏, yields : 

𝜏𝑘 =
1

𝜔0
𝑐𝑜𝑠−1

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

2−𝐷3𝐷6

𝐷4
2𝜔0

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

2+𝐷6
2 +

2𝑘𝜋

𝜔0
                                                     (41) 

where, k = 0,1,2,…. Hence we get the corresponding 𝜏𝑘 > 0 for which system (3) has two wholly 

imaginary roots ±𝜔0. 

Let 𝜆(𝜏) = 𝜇(𝜏) + 𝑖𝜔(𝜏) be a root of Equation (37) near 𝜏 = 𝜏𝑘 with  𝜇(𝜏𝑘) = 0 and 𝜔(𝜏𝑘) =

𝜔0. Then comes the next theorem: 

 

Theorem2. The C.E. (37), roots satisfy the following transversality requirement; 

[
𝑑(𝑅𝑒𝜆(𝜏))

𝑑𝜏
]
𝜏=𝜏𝑘

> 0                                                                                                                     (42a) 

Provided that  

2𝐷2 > 𝐷1
2                                                                                                                                  (42b) 

𝐷2
2 > 2𝐷1𝐷3                                                                                                                              (42c) 

2𝐷4𝐷6 > 2𝐷4
2𝜔0

2 + 𝐷5
2                                                                                                             (42d) 

 

Proof. By substituting 𝜆(𝜏) in Equation (37) and differentiating the resulting equation in 𝜏,we can 

get  

[3𝜆2 + 2𝐷1𝜆 + 𝐷2 + (2𝐷4𝜆 + 𝐷5)𝑒
−𝜆𝜏 − 𝜏(𝐷4𝜆

2 + 𝐷5𝜆 + 𝐷6)𝑒
−𝜆𝜏]

𝑑𝜆

𝑑𝜆
= 𝜆(𝐷4𝜆

2 + 𝐷5𝜆 + 𝐷6)𝑒
−𝜆𝜏 

                                                                                                                                        (43) 

Thus,  

(
𝑑𝜆

𝑑𝜏
)
−1

=
3𝜆2+2𝐷1𝜆+𝐷2

𝜆(𝐷4𝜆2+𝐷5𝜆+𝐷6)𝑒−𝜆𝜏
+

2𝐷4𝜆+𝐷5

𝜆(𝐷4𝜆2+𝐷5𝜆+𝐷6)
−

𝜏

𝜆
                                                                    (44) 

Since for = 𝜏0 , and 𝜆 = 𝑖𝜔0, we have got  

 

[
𝑑𝜆

𝑑𝜏
]
𝜏=𝜏0

−1

=
(𝐷2 − 3𝜔0

2) + 2𝑖𝐷1𝜔0

−𝜔0
2(𝜔0

2 − 𝐷2) + 𝑖𝜔0(𝐷1𝜔0
2 − 𝐷3)

+
𝐷5 + 2𝑖𝐷4𝜔0

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

2)
−
𝑖𝜏0
𝜔0

 

Now, since  

𝑠𝑖𝑔𝑛 [
𝑑(𝑅𝑒𝜆)

𝑑𝜏
]
𝜏=𝜏0

= 𝑠𝑖𝑔𝑛 [𝑅𝑒 (
𝑑𝜆

𝑑𝜏
)
−1

]
𝜆=𝑖𝜔0

                                                                              (45) 

It is clear that: 



IHJPAS. 2025, 38 (1) 
 

455 
 

 

𝑅𝑒 [
(𝐷2−3𝜔0

2)+2𝑖𝐷1𝜔0

−𝜔0
2(𝜔0

2−𝐷2)+𝑖𝜔0(𝐷1𝜔0
2−𝐷3)

] =
2𝐷1(𝐷1𝜔0

2−𝐷3)−(𝜔0
2−𝐷2)(𝐷2−3𝜔0

2)

𝜔0
2(𝜔0

2−𝐷2)
2
+(𝐷1𝜔0

2−𝐷3)
2

𝑅𝑒 [
𝐷5+2𝑖𝐷4𝜔0

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

2)
] =

2𝐷4(𝐷6−𝐷4𝜔0
2)−𝐷5

2

𝐷5
2𝜔0

2+(𝐷6−𝐷4𝜔0
2)
2                                  

𝑅𝑒 [
𝑖𝜏0

𝜔0
] = 𝑧𝑒𝑟𝑜                                                                             

 

Hence, we have  

 𝑅𝑒 [
𝑑𝜆

𝑑𝜏
]
𝜏=𝜏0

−1

=
2𝐷1(𝐷1𝜔0

2−𝐷3)−(𝜔0
2−𝐷2)(𝐷2−3𝜔0

2)

Ψ0
+
2𝐷4(𝐷6−𝐷4𝜔0

2)−𝐷5
2

Ψ1
 

Where  

 
Ψ0 = 𝜔0

2(𝜔0
2 − 𝐷2)

2 + (𝐷1𝜔0
2 − 𝐷3)

2 > 0 ,

Ψ1 = 𝐷5
2𝜔0

2 + (𝐷6 − 𝐷4𝜔0
2)2 > 0 .                

 

We obtain [
𝑑(𝑅𝑒𝜆(𝜏))

𝑑𝜏
]
𝜏=τ0

> 0    (42b- 42d). 

This outcome demonstrates how the roots of C.E. (37), as 𝜏 passes through τ0, traverse the 

imaginary axis from left to right. As a result, the system (3) experiences an H.B. at 𝜏 = τ0 and 

loses its stability. 

 

For the (SIEP), Equation (4) reduces to                

             𝐽(𝐸5) = (𝑐𝑖𝑗)4×4   ;  𝑖, 𝑗 = 1,2,3,4 

Here 

𝑐11 = −[𝑅1 + 𝑅2𝑒
−𝜆𝜏]  , 𝑐12 =

𝑛𝛼𝑆∗

(1 + 𝑛𝑉1
∗)2
    , 𝑐14 = −𝑅3𝑒

−𝜆𝜏,   𝑐21 =
𝛼

1 + 𝑛𝑉1
∗ ,                                                                               

 𝑐22 = −(
𝑛𝛼𝑆∗

(1 + 𝑛𝑉1
∗)2

+ 𝛾 + 𝛽2𝐼
∗ + 𝜇) , 𝑐24 = −𝛽2𝑉1

∗ , 𝑐32 = 𝛾,   𝑐33 = −(𝛽3𝐼
∗ + 𝜇) ,                                                                  

  
 𝑐34 = −𝛽3𝑉2

∗   , 𝑐41 = 𝑅2𝑒
−𝜆𝜏, 𝑐42 = 𝛽2𝐼

∗  , 𝑐43 = 𝛽3𝐼
∗ ,   𝑐44 = 𝑅3𝑒

−𝜆𝜏 + 𝑅4,                                                                                                                            
𝑐13 = 𝑐31 = 𝑐23 = 0.                                                                                                                                                                                           

  

                                                                                                                                                     (46) 

Where, 

𝑅1 =
𝛼

1 + 𝑛𝑉1
∗ + 𝜇 ,  𝑅2 = 𝛽1𝐼

∗  𝑎𝑛𝑑 𝑅3 = 𝛽1𝑆
∗ , 𝑅4 = 𝛽2𝑉1

∗ + 𝛽3𝑉2
∗ − (𝜇 + 𝜇1 + 𝜃)  

The C.E. of 𝐽(𝐸5) is  

[𝜆4 + 𝐶1𝜆
3 + 𝐶2𝜆

2 + 𝐶3𝜆 + 𝐶4 + (𝐶5𝜆
3 + 𝐶6𝜆

2 + 𝐶7𝜆 + 𝐶8)𝑒
−𝜆𝜏] = 0                       (47)                                      



IHJPAS. 2025, 38 (1) 
 

456 
 

With 
𝐶1 = 𝑅1 − (𝑅4 + 𝑐22 + 𝑐33), 𝐶2 = (𝑐22 + 𝑐33)(𝑅4 − 𝑅1) +   𝑐22𝑐33 −                                      

                                                                                                (𝑅1𝑅4 + 𝑐12𝑐21 + 𝑐24𝑐42 + 𝑐34𝑐43),
                                                                                                          

𝐶3 = (𝑐22𝑐33 + 𝑅4𝑐22 + 𝑅4𝑐33 − (𝑐24𝑐42 + 𝑐34𝑐43))𝑅1 + (𝑐12𝑐21 − 𝑐22𝑐33)𝑅4 +                                                        

                                                                                (𝑐12𝑐21 + 𝑐24𝑐42)𝑐33 + (𝑐22𝑐34 − 𝑐24𝑐32)𝑐43,    

𝐶4 = (𝑐22𝑐34𝑐43 + 𝑐24𝑐33𝑐42 − 𝑅4𝑐22𝑐33 − 𝑐24𝑐32𝑐43)𝑅1 + 𝑐12𝑐21(𝑐34𝑐43 − 𝑐33𝑅4),

 𝐶5 = 𝑅2 − 𝑅3 , 𝐶6 = (𝑐22 + 𝑐33)(𝑅3 − 𝑅2) − (𝑅1𝑅3 + 𝑅2𝑅4),                                         
                                                    

                                                                  

    

  𝐶7 = (𝑐33(𝑐22 + 𝑅4) − 𝑐24(𝑐42 + 𝑐12) − 𝑐34𝑐43)𝑅2 +                                                                                                 

                                                                 
(𝑅3𝑐22 + 𝑅4𝑐22 + 𝑅3𝑐33)𝑅1 + (𝑐12𝑐21 + 𝑐21𝑐42 − 𝑐22𝑐33)𝑅3,

                    
  

                                                                                  

𝐶8 = (𝑐21𝑐32𝑐43 − (𝑅1𝑐22 + 𝑐12𝑐21 + 𝑐21𝑐42)𝑐33)𝑅3 +                    

(𝑐22𝑐34𝑐43 + 𝑐24𝑐33𝑐42 + 𝑐12𝑐24𝑐33 − (𝑅4𝑐22𝑐33 + 𝑐24𝑐32𝑐43))𝑅2
                                                                                                                                                        

 

Now, when = 0 , Equation (47) is:  

[𝜆4 + (𝐶1 + 𝐶5)𝜆
3 + (𝐶2 + 𝐶6)𝜆

2 + (𝐶3 + 𝐶7)𝜆 + 𝐶4 + 𝐶8] = 0                                             (48) 

all the eigenvalues of Equation (48) will be present in the left half plane and the SIEP is L.A.S. of 

system (3) under the following condition: 

2(𝛽1𝑆
∗ + 𝛽2𝑉1

∗ + 𝛽3𝑉2
∗) < 𝜇 + 𝜇1 + 𝜃                                                                                      (49) 

Now, for 𝜏 > 0, then from Equation (47) has two wholly imaginary roots, namely 𝜆 = ±𝑖𝜔 ( 𝜔 >

0). 

By substituting 𝜆 = ±𝑖𝜔  in Equation (47) and separating the real and imaginary parts, which 

gives 

(𝐶8 − 𝐶6𝜔
2)𝑐𝑜𝑠𝜔𝜏 + (𝐶7𝜔 − 𝐶5𝜔

3)𝑠𝑖𝑛𝜔𝜏 = 𝐶2𝜔
2 − 𝜔4 − 𝐶4,

(𝐶7𝜔 − 𝐶5𝜔
3)𝑐𝑜𝑠𝜔𝜏 + (𝐶6𝜔

2 − 𝐶8)𝑠𝑖𝑛𝜔𝜏 = 𝐶1𝜔
3 − 𝐶3𝜔.       

                                              (50) 

Through squaring and adding both equations in Equation (50), we get  

𝜔8 + ℎ1𝜔
6 + ℎ2𝜔

4 + ℎ3𝜔
2 + ℎ4 = 0                                                                                        (51) 

Further, by letting 𝜐 = 𝜔2 in Equation (51) 

ℎ(𝜐) = 𝜐4+ℎ1𝜐
3 + ℎ2𝜐

2 + ℎ3𝜐 + ℎ4 = 0,                                                                                (52a) 

where 

 

ℎ1 = 𝐶1
2 − 2𝐶2 − 𝐶5

2,                                                                                      

ℎ2 = 2𝐶5𝐶7 + 𝐶2
2 + 2𝐶4 − 2𝐶1𝐶3 − 𝐶6

2,                                                    

ℎ3 = 2𝐶6𝐶8 + 𝐶3
2 − 𝐶7

2 − 2𝐶1𝐶4,

ℎ4 = 𝐶4
2 − 𝐶8

2.                                  
                                                                 

  

Straightforward computation shows that due to the following condition 

𝐶4 < 𝐶8                                                                                                                                     (52b) 

 

we obtain ℎ4 < 0. So, there exists a unique positive root in accordance with Descartes' rule of 

signs 𝜐0 = 𝜔0
2 providing Equation (52a). Which is, Equation (51) has a positive 𝜔0. As a result, 

Equation (47) has at least two roots ±𝜔0 that are purely imaginary and which correspond to the 

time delay 𝜏. Additionally, by substituting 𝜔0 in (50) and solving the system for 𝜏, yields : 

𝜏𝑘 =          
1

𝜔0
𝑐𝑜𝑠−1

(𝐶6−𝐶1𝐶5)𝜔0
6+(𝐶1𝐶7+𝐶5+𝐶3−𝐶6𝐶2−𝐶8)𝜔0

4+(𝐶6𝐶4+𝐶2𝐶8−𝐶3𝐶7)𝜔0
2−𝐶4𝐶8

𝐶5
2𝜔0

6+(𝐶6
2−2𝐶5𝐶7)𝜔0

4+(𝐶7
2−2𝐶5𝐶8)𝜔0

2+𝐶8
2 +

                                                                                                                                                 
2𝑘𝜋

𝜔0
                                       

                                                                                                                                                     (53) 



IHJPAS. 2025, 38 (1) 
 

457 
 

where, k=0,1,2,…. Hence we get the corresponding 𝜏𝑘 > 0 for which system (3) has two wholly 

imaginary roots ±𝜔0. 

Let 𝜆(𝜏) = 𝜇(𝜏) + 𝑖𝜔(𝜏) be a root of C.E.(47) near 𝜏 = 𝜏𝑘with  𝜇(𝜏𝑘) = 0 and 𝜔(𝜏𝑘) = 𝜔0. Then 

comes the next theorem. 

 

Theorem3. The C.E. (47), roots satisfy the following transversality requirement; 

[
𝑑(𝑅𝑒𝜆(𝜏))

𝑑𝜏
]
𝜏=𝜏𝑘

> 0                                                                                                                    (54a) 

Provided that  

4𝜔0
6 + 2𝐶1

2𝜔0
5 + 2(𝐶2

2 + 2𝐶4)𝜔0
2 + 2𝐶2𝐶3𝜔0 >            

                                                   6𝐶2𝜔0
4 + 𝐶3(𝐶1 + 4) + 2𝐶2𝐶4 ,

                                                    (54b) 

2𝐶5
2𝜔0

5 + 2𝐶6
2𝜔0

2 > 𝐶5𝐶7𝜔0
3 + 2𝐶6(𝐶7 + 𝐶8).                                                                          (54c)                                        

 

Proof. By substituting 𝜆(𝜏) in Equation (47) and differentiating the resulting equation in 𝜏,we can 

get  

[4𝜆3 + 3𝐶1𝜆
2 + 2𝐶2𝜆 + 𝐶3 + (3𝐶5𝜆

2 + 2𝐶6𝜆 + 𝐶7)𝑒
−𝜆𝜏 −                    

                                          𝜏(𝐶5𝜆
3 + 𝐶6𝜆

2 + 𝐶7𝜆 + 𝐶8)𝑒
−𝜆𝜏]

𝑑𝜆

𝑑𝜆
=        

                                                                    𝜆(𝐶5𝜆
3 + 𝐶6𝜆

2 + 𝐶7𝜆 + 𝐶8)𝑒
−𝜆𝜏.

 

                                                                                                                                         (55) 

Thus,  

(
𝑑𝜆

𝑑𝜏
)
−1

=
4𝜆3+3𝐶1𝜆

2+2𝐶2𝜆+𝐶3

𝜆(𝐶5𝜆3+𝐶6𝜆2+𝐶7𝜆+𝐶8)𝑒−𝜆𝜏
+

3𝐶5𝜆
2+2𝐶6𝜆+𝐶7

𝜆(𝐶5𝜆3+𝐶6𝜆2+𝐶7𝜆+𝐶8)
−

𝜏

𝜆
                                                       (56) 

Since for = 𝜏0 , and 𝜆 = 𝑖𝜔0, we have got  

[
𝑑𝜆

𝑑𝜏
]
𝜏=𝜏0

−1

=
𝐶3 − 3𝐶1𝜔0

2 + 𝑖𝜔0(2𝐶2 − 4𝜔0
2)

𝐶1𝜔0
4 + 𝑖𝜔0(𝜔0

4 + 𝐶4 − 𝐶2𝜔0
2 − 𝐶3𝜔0)

+
−3𝐶5𝜔0

2 + 𝐶7 + 2𝑖𝐶6𝜔0

𝐶5𝜔0
4 + 𝑖𝜔0(𝐶7 + 𝐶8 − 𝐶6𝜔0

2)
−
𝑖𝜏0
𝜔0

 

  
Now, since  

𝑠𝑖𝑔𝑛 [
𝑑(𝑅𝑒𝜆)

𝑑𝜏
]
𝜏=𝜏0

= 𝑠𝑖𝑔𝑛 [𝑅𝑒 (
𝑑𝜆

𝑑𝜏
)
−1

]
𝜆=𝑖𝜔0

                                                                                (57) 

It is clear that : 

 

𝑅𝑒 [
𝐶3−3𝐶1𝜔0

2+𝑖𝜔0(2𝐶2−4𝜔0
2)

𝐶1𝜔0
4+𝑖𝜔0(𝜔0

4+𝐶4−𝐶2𝜔0
2−𝐶3𝜔0)

] =
𝐶1𝐶3𝜔0

3−3𝐶1
2𝜔0

5+(2𝐶2−4𝜔0
2)(𝜔0

4+𝐶4−𝐶2𝜔0
2−𝐶3𝜔0)

𝜔0(𝐶1
2𝜔0

6+(𝜔0
4+𝐶4−𝐶2𝜔0

2−𝐶3𝜔0)
2
)

,

𝑅𝑒 [
−3𝐶5𝜔0

2+𝐶7+2𝑖𝐶6𝜔0

𝐶5𝜔0
4+𝑖𝜔0(𝐶7+𝐶8−𝐶6𝜔0

2)
] =

𝐶5𝐶7𝜔0
3−3𝐶5

2𝜔0
5+2𝐶6(𝐶7+𝐶8−𝐶6𝜔0

2)

𝐶5
2𝜔0

6+(𝐶6+𝐶7−𝐶8𝜔0
2)
2 ,                                    

𝑅𝑒 [
𝑖𝜏0

𝜔0
] = 𝑧𝑒𝑟𝑜.                                                                                                                  

 

Hence, we have 

 

𝑅𝑒 [
𝑑𝜆

𝑑𝜏
]
𝜏=𝜏0

−1

= 
𝐶1𝐶3𝜔0

3−3𝐶1
2𝜔0

5+(2𝐶2−4𝜔0
2)(𝜔0

4+𝐶4−𝐶2𝜔0
2−𝐶3𝜔0)

Ψ0
+

                                                                                     
𝐶5𝐶7𝜔0

3−3𝐶5
2𝜔0

5+2𝐶6(𝐶7+𝐶8−𝐶6𝜔0
2)

Ψ1
.

  

Where  

 
Ψ0 = 𝜔0(𝐶1

2𝜔0
6 + (𝜔0

4 + 𝐶4 − 𝐶2𝜔0
2 − 𝐶3𝜔0)

2) > 0 ,

Ψ1 = 𝐶5
2𝜔0

6 + (𝐶6 + 𝐶7 − 𝐶8𝜔0
2)2 > 0 .                           

 



IHJPAS. 2025, 38 (1) 
 

458 
 

We obtain [
𝑑(𝑅𝑒𝜆(𝜏))

𝑑𝜏
]
𝜏=τ0

> 0    under (54a-54b). 

This outcome demonstrates how the roots of C.E. (47), as 𝜏 passes through τ0, traverse the 

imaginary axis from left to right. As a result, the system (3) experiences an H.B. at 𝜏 = τ0 and 

loses its stability.  

 

4. Numerical Simulation and Discussion 

     numerical simulation is employed in this section to show the findings of our investigation. In 

this section, the following hypothetical parameters have been chosen: 

Λ = 0.042 , 𝛽1 = 0.03 , 𝛽2 = 0.02 , 𝛽3 = 0.01 , 𝛼 = 0.3 ,
𝑛 = 50 ,   𝜃 = 0.1 , 𝛾 = 0.01 ,  𝜇1 = 0.03 , 𝜇 = 0.00015 , 𝜏 = 25.2225.

                                (58)    

Investigated is the dynamical behavior of system (3) near the SIEP point when the T.D. is 

increased. For the set of parameter values provided by (58), the system (3) is numerically solved, 

and Figure 1 shows the trajectory of the system (3). 

 
(a)                                                                             (b) 

 

 
         (c)                                                                                  (d)  

Figure 1.Periodic solution near (SIEP) of system (3) for (58) (a) Periodic solution near (SIEP). (b),(c) and (d) 3D-

periodic solution. 

 

Equation (58) is used to observe that for the provided data, with 𝜏 = 0  system (3) has a locally 

asymptotically stable (L.A.S.) to SIEP as shown in Figure 2. 



IHJPAS. 2025, 38 (1) 
 

459 
 

             

 Figure 2. System (3) trajectories using the information provided by (58) with 𝜏 = 0 approach to SIEP. 

 

Equation (58) is used to observe that for the provided data, with  𝛼 = 0 and 𝛽1 = 0.0003 system 

(3) has a locally asymptotically stable (L.A.S.) to FEP as shown in Figure 3. 

 

Figure 3. System (3) trajectories using the information provided by (58) with  𝛼 = 0 and 𝛽1 = 0.0003 approach to 

FEP. 

 

Equation (58) is used to observe that for the provided data, with  𝛾 = 0 , 𝛽1 = 0.0003 𝑎𝑛𝑑 𝛽2 =

0.0002  system (3) has a L.A.S. to SEP as shown in   Figure 4. 



IHJPAS. 2025, 38 (1) 
 

460 
 

 
Figure 4. System (3) trajectories using the information provided by equation (58) with  𝛾 = 0 , 𝛽1 = 0.0003 𝑎𝑛𝑑 𝛽2 =

0.0002  approach to SEP. 

 

We talk about how the time delay affects how the system behaves close to the TEP. 

For τ = 10 < τ0 = 13.5  and 𝛼 = 0  with the set of data in equation (58) TEP is still L.A.S.as 

shown in Figure 5.  

 

Figure 5. System (3) trajectories using the information provided by (58) with   𝛼 = 0 and τ = 10 approach to TEP. 

 

Now, for  τ0 = 13.5  and 𝛼 = 0  with the set of data in equation (58) a H.B. occurs at TEP as 

shown in Figure 6. 



IHJPAS. 2025, 38 (1) 
 

461 
 

 

(a)                                                              (b) 

 

Figure 6. System (3) trajectories using the information provided by (58) with  𝛼 = 0 and τ = 13.5. (a) Periodic  

solution  near TEP. (b) 3D- periodic solution. 

Equation (58) is used to observe that for the provided data, with 𝛽1 = 0.0003 𝑎𝑛𝑑 𝛽2 = 0.0001  

system (3) has a L.A.S. to FOEP as shown in Figure 7. 

 

Figure 7. System (3) trajectories using the information provided by equation (58) with   𝛽1 = 0.0003 𝑎𝑛𝑑 𝛽2 =

0.0001  approach to FOEP. 

 

We talk about how the time delay affects how the system behaves close to the FIEP.  

For τ = 10 < τ0 = 19  and 𝛾 = 0  with the set of data in (58) FIEP is still L.A.S. as shown in 

Figure 8. 



IHJPAS. 2025, 38 (1) 
 

462 
 

  

Figure 8. System (3) trajectories using the information provided by (58) with   𝛾 = 0 and τ = 10 approach to FIEP. 

 

Now, for  τ0 = 19  and 𝛾 = 0   with the set of data in equation (58) a H.B. occurs at FIEP as shown 

in Figure 9. 

 

(a)                                                             (b) 

Figure 9. System (3) trajectories using the information provided by (58) with  𝛾 = 0 and τ = 19. (a) Periodic  

solution  near FIEP. (b) 3D- periodic solution. 

 

3. Conclusion 

      A mathematical model was proposed and studied for the effect of two stages of the vaccine 

against the Coronavirus, which includes a time delay for the period of infection with the virus. The 

suggested system has six equilibrium points, namely FEP, SEP, TEP, FOEP, FIEP and SIEP. The 

FEP, SEP and FOEP are seen absolutely stable for all 𝜏 ≥ 0. The TEP, FIEP and SIEP is 

asymptotically stable for τ ∈ [0, τ0), but an H.B. occurs when 𝜏 = τ0. The TEP is still L.A.S. For 

τ = 10 < τ0 = 13.5  and 𝛼 = 0  with the set of data in equation (58) .While, for  τ = τ0 =

13.5  and 𝛼 = 0, a H.B. is demonstrated near  TEP. For τ = 10 < τ0 = 19  and 𝛾 = 0  with the 



IHJPAS. 2025, 38 (1) 
 

463 
 

set of data in equation (58), the FIEP is still L.A.S. . While, for  τ = τ0 = 19  and 𝛾 = 0, a H.B. 

is demonstrated near  FIEP . 

 

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. 

 

Ethical Clearance 

     Ethics of scientific research were carried out in accordance with international conditions. 

 

References 

1. Ferguson NM, Donnelly CA, Anderson RM. Transmission intensity and impact of control policies on 

the foot and mouth epidemic in Great Britain. Nature. 2001 Oct 4;413(6855):542-8. 

https://doi.org/10.1038/35097116 

2. 2.Cauchemez S, Fraser C, Van Kerkhove MD, Donnelly CA, Riley S, Rambaut A, Enouf V, van der 

Werf S, Ferguson NM. Middle East respiratory syndrome coronavirus: quantification of the extent of 

the epidemic, surveillance biases, and transmissibility. The Lancet infectious diseases. 2014 Jan 

1;14(1):50-6. https://doi.org/10.1016/S1473-3099(13)70304-9 

3. 3.Mohsen AA, Al-Husseiny HF, Zhou X, Hattaf K. Global stability of COVID-19 model involving the 

quarantine strategy and media coverage effects. AIMS public Health. 

2020;7(3):587. https://doi.org/10.3934/publichealth.2020047 

4. 4.AL-Husseiny HF, Mohsen AA, Zhou X. The Effect of Individuals Asymptomatic (Carrier) on The 

Dynamical Behavior Of a COVID-19 Virus. Ibn AL-Haitham Journal for Pure and Applied Sciences. 

2021 Apr 20;34(2):42-55. https://doi.org/10.30526/34.2.2624 

5. 5.Hattaf K, Mohsen AA, Harraq J, Achtaich N. Modeling the dynamics of COVID-19 with carrier 

effect and environmental contamination. International Journal of Modeling, Simulation, and Scientific 

Computing. 2021 Jun 19;12(03):2150048. 

https://doi.org/10.1142/S1793962321500483 

6. Abdulkadhim MM, Al-Husseiny HF. Global stability and bifurcation of a COVID-19 virus modeling 

with possible loss of the immunity. InAIP conference Proceedings 2020 Oct 27 (Vol. 2292, No. 1). AIP 

Publishing. https://doi.org/ 10.1063/5.0030669 

7. Zuo L, Liu M, Wang J. The impact of awareness programs with recruitment and delay on the spread of 

an epidemic. Mathematical problems in Engineering. 

2015;2015(1):235935.https://doi.org/10.1155/2015/235935 

8. Sirijampa A, Chinviriyasit S, Chinviriyasit W. Hopf bifurcation analysis of a delayed SEIR epidemic 

model with infectious force in latent and infected period. Advances in Difference Equations. 2018 

Dec;2018:1-24. https://doi.org/10.1186/s13662-018-1805-6 

https://doi.org/10.1038/35097116
https://doi.org/10.1016/s1473-3099(13)70304-9
https://doi.org/10.3934/publichealth.2020047
https://doi.org/10.30526/34.2.2624
https://doi.org/10.1142/S1793962321500483
http://dx.doi.org/10.1063/5.0030669
https://doi.org/10.1155/2015/235935
https://advancesindifferenceequations.springeropen.com/articles/10.1186/s13662-018-1805-6


IHJPAS. 2025, 38 (1) 
 

464 
 

9. Yaseen RM, Mohsen AA, Al-Husseiny HF, Hattaf K. Stability and Hopf bifurcation of an 

epidemiological model with effect of delay the awareness programs and vaccination: analysis and 

simulation. Commun. Math. Biol. Neurosci.. 2023 Mar 4;2023:Article-ID. 

10. Yin Z, Yu Y, Lu Z. Stability analysis of an age-structured SEIRS model with time delay. Mathematics. 

2020 Mar 23;8(3):455. https://doi.org/10.3390/math8030455 

11. Mohsen AA, AL-Husseiny HF, Naji RK. The dynamics of Coronavirus pandemic disease model in the 

existence of a curfew strategy. Journal of Interdisciplinary Mathematics. 2022 Aug 18;25(6):1777-97. 

https://doi.org/10.1080/09720502.2021.2001139 

12. Zhang Z, Upadhyay RK. Dynamical analysis for a deterministic SVIRS epidemic model with Holling 

type II incidence rate and multiple delays. Results in Physics. 2021 May 1;24:104181. 

https://doi.org/10.1016/j.rinp.2021.104181 

13. 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 Jul 

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

14. Shafeeq SK, Abdulkadhim MM, Mohsen AA, Al-Husseiny HF, Zeb A. Bifurcation analysis of a 

vaccination mathematical model with application to COVID-19 pandemic. Commun. Math. Biol. 

Neurosci.. 2022 Dec 9;2022:Article-ID. 

15. Hussien RM, Naji RK. The dynamics of the SEIR epidemic model under the influence of delay. 

Commun. Math. Biol. Neurosci.. 2022 Jun 27;2022:Article-ID. 

16. Mohsen AA, AL-Husseiny HF, Hattaf K, Boulfoul B. A mathematical Model for the Dynamics of 

COVID-19 Pandemy Involving the Infective Immigrants. Iraqi Journal of Science. 2021;62(1):295-

307.DOI: https://doi.org/10.24996/ijs.2021.62.1.28 

17. Naji R, Muhseen A. Stability analysis with bifurcation of an SVIR epidemic model involving 

immigrants. Iraqi journal of Science. 2013;54(2):397-408. 

18. Ahmed LS, AL-Husseiny HF. Dynamical Behavior of an eco-epidemiological Model involving 

Disease in predator and stage structure in prey. Iraqi Journal of Science. 2019 Aug 26:1766-

82.DOI: https://doi.org/10.24996/ijs.2019.60.8.14 

19. Helal MM, Yaseen RM, Mohsen AA, AL-Husseiny HF, Sabbar Y. Dynamics of a Social Model for 

Marriage and Divorce Relationship with Fear Effect. Malaysian Journal of Mathematical Sciences. 

2024 Jun 1;18(2). https://doi.org/10.47836/mjms.18.2.04 

20. Saadi RR, Al-Husseiny HF. A Mathematical Modeling of the Vaccination Effect on the SARS-CoV-2 

Transmission: Analysis and Simulation. Iraqi Journal of Science. 2024 Mar 29;65(3):1548-70.  

https://doi.org/ 10.24996/ijs.2024.65.3.31 

https://doi.org/10.3390/math8030455
https://doi.org/10.1080/09720502.2021.2001139
https://doi.org/10.1080/09720502.2021.2001139
http://dx.doi.org/10.1016/j.rinp.2021.104181
https://doi.org/10.3390/math10132344
https://doi.org/10.3390/math10132344
https://doi.org/10.24996/ijs.2021.62.1.28
https://doi.org/10.24996/ijs.2019.60.8.14
http://dx.doi.org/10.47836/mjms.18.2.04
https://doi.org/%2010.24996/ijs.2024.65.3.31