Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

2684 
https://internationalpubls.com 

An Analytical and Numerical simulation of Fuzzy Delay Differential 

Equation Model of HIV infection of CD4+T-Cells 

T.Muthukumar1, R.Ramesh2∗, L.Chitra3, C. Kothai Andal4, K. Kalaiselvi5 

1,2,5Department of Mathematics 

Dr.Mahalingam College of Engineering and Technology, Pollachi-642 003, Tamilnadu, India. 

vtmuthukumar@gmail.com 

3Department of Electrical and Electronics Engineering, 

Dr.Mahalingam College of Engineering and Technology,Pollachi-642 003, Tamilnadu, India. 

4Department of Electrical and Electronics Engineering, AMC Engineering College, Bengaluru, Karnataka, India. 

*Corresponds: rameshwaran141@gmail.com 

 

Article History: 

Received: 12-01-2025 

Revised: 15-02-2025 

Accepted: 01-03-2025 

Abstract:  

           A new fuzzy mathematical model of HIV infection of CD4+ T-Cells consisting 

of fuzzy Delay Differential Equations (FDDEs) has been suggested with the appropriate 

theory and analysis. In this fuzzy dynamical system, we consider three equations having 

parameters uninfected, viral and infected cells. Therefore, utilizing a fuzzy dynamical 

system and 𝛼− cuts, we prepare to suppress the HIV epidemic in this approach. We 

transform this system into a fuzzy linear system of differential equations, and then we 

endeavor to solve these FDDEs using the Runge-Kutta method at the order of five. In 

addition, we presented the stability of fuzzy differential equations before the numerical 

analysis 

Keywords: HIV infection; Fuzzy Delay Differential equation; Stability; Numerical 

solution; Runge-Kutta Method of order five. 

     Mathematics Subject Classification: 74H15, 34A07 

1. Introduction 

          HIV infection of CD4+ cells has been a subject of discussion for many studies.  Perelson [14] 

provided a simple explanation of how to form a model for interaction between HIV (Human immune-

deficiency virus). Even though the interactions between the human immune system and HIV are 

extremely complicated, we still don’t understand the pathogenicity. They developed the HIV model 

and evaluated the model’s behavior. Uninfected CD4+ T-cells, latently infected CD4+ T-cells, infected 

CD4+ T-cells, and the free virus was listed as the four components. The fluctuations of a density of 

uninfected target cells (𝑈𝑐), infected target cells (𝐼𝑐), and free virus (𝐹𝑣) were explored in the basic 

differential equation model developed by [14] to analyze viral dynamics. Numerous researchers, 

including [3], [5], [6], [10], [18], [19], [22], detailed the model utilising ODE techniques (ODE). That 

model does not include an intracellular time delay between a cell becoming affected and a virus-

producing a new infection. Herz et al., [8], assumed that virus generation lags behind cell infection by 

a delay of τ . Finally, there are two techniques to solve the system of equations. This reveals that the 

density of cells that were freshly infected at time 𝑡 −  𝜏 

and are still alive at time t determines the recruitment of virus-producing cells at time 𝑡. The first 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

2685 
https://internationalpubls.com 

one is analytical when it is possible, and the second is numerical methods since they enable the 

prediction of virus concentrations and target cell behaviour. 

             We will discuss a fuzzy differential equation model of HIV infection of CD4+ T-cells with 

delay using the fifth order Runge-Kutta method in this study. Rebecca et al. [4] have described the 

most basic model of HIV infection of CD4+ T-cells. We must turn the nonlinear fuzzy delay model 

into a fuzzy linear model in this study. Then it will provide a strategy for determining the analytical 

solution of the fuzzy linearized system. 

Variables and parameters for viral spread[4] 

Terms Description of variables and constants Values 

 Dependent variables  

𝑈𝑐 Uninfected CD4+ T-cell population size 1000mm−3 

𝐼𝑐 

𝐹𝑣 

Infected CD4+ T-cell density 

Initial density of HIV RNA 

0 

10−3 

  Parameters and Constants  
𝜇𝑢𝑐  

 
  Rate of natural death for CD4+T cells 0.02day−1 

µ

I 
 
c 

Infectious people die in droves for CD4+ T-cells 0.26day−1 

µb  Rate of lytic death in infected cells 0.24day−1 

µ

F 
 
v 

Rate of free virus death 2.4day−1 

k1  Rate CD4+ T-cells become infected with virus 2.4 × 10−5mm3day−1 
′ 

k  The rate at which infected cells activate 2 × 10−5mm3day−1 

r  Growth rate of CD+ T-cell population 0.03day−1 

N  The number of virions that infected CD4+ T cells produce Varies 

U

c 
 
max 

Maximum number of CD4+ T-cells in a population 1500mm−3 

s  Uninfected CD4+ T-cells are the source term 10(day)−1(mm−3) 

  derived quantities  

T0  The population of CD4+ T cells for HIV-negative patients 1000mm−3 

 

2. HIV Infection of CD4+cells with a Fuzzy Delay Model 

Consider the FDDEs model of CD4+T-cell HIV infection, 

{
 
 

 
 
𝑑

𝑑𝑡
[𝑈𝑐(𝑡)]

𝛼 = [𝑠]𝛼 − [𝜇𝑡𝑈𝑐(𝑡)]
𝛼 + [𝑟𝑈𝑐(𝑡)]

𝛼 (1 −
[𝑈𝑐(𝑡)]

𝛼+[𝐼𝑐(𝑡)]
𝛼

𝑈𝑐𝑚𝑎𝑥
) − 𝑘[𝐹𝑣(𝑡)]

𝛼[𝑈𝑐(𝑡)]
𝛼,

𝑑

𝑑𝑡
[𝐼𝑐(𝑡)]

𝛼 = 𝑘′1[𝐹𝑣(𝑡 − 𝜏)]
𝛼[𝑈𝑐(𝑡 − 𝜏)]

𝛼 − 𝜇𝐼[𝐹𝑣(𝑡)]
𝛼,

𝑑

𝑑𝑡
[𝐹𝑣(𝑡)]

𝛼 = [𝑁𝜇𝑏𝐼𝑐(𝑡)]
𝛼 − 𝑘1[𝐹𝑣(𝑡)]

𝛼[𝑈𝑐(𝑡)]
𝛼 − 𝜇𝐹𝑣[𝐹𝑐(𝑡)]

𝛼,

(2.1)  

 

with the fuzzy initial function is given by 

   

[𝑈𝑐(𝑡)]
𝛼 = [1001 − 𝑒𝑡]𝐹1(𝛼), 

[𝐼𝑐(𝑡)]
𝛼 = [0]𝐹2(𝛼), 

[𝐹𝑣(𝑡)]
𝛼 = [1.001 − 𝑒𝑡]𝐹3(𝛼), 𝑡 ∈ [−𝜏, 0], 

where ,  

𝐹1(𝛼) = [1000𝛼 + 850(1 − 𝛼),1000𝛼 + 1150(1 − 𝛼)] 

𝐹2(𝛼) = [5𝛼 + 4(1 − 𝛼), 5𝛼 + 6(1 − 𝛼)], 

𝐹3(𝛼) = [7000𝛼 + 6750(1 − 𝛼), 7000𝛼 + 7250(1 − 𝛼)], 𝛼 ∈ [0,1]. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

2686 
https://internationalpubls.com 

            Where [𝑈𝑐(𝑡)]
𝛼 symbolizes the number of healthy CD4+T-cells at time 𝑡. [𝐼𝑐(𝑡)]

𝛼 reflects 

the accumulation of carrying the virus CD4+T-cells at time t, [𝐹𝑣(𝑡)] corresponds to the number 

of free HIV at time t, and the constant of proportionality 𝜏, which reflects the time frame of the 

delay, is also present. The parameters are as follows: s represents the precursors source of CD4+T-

cells, µ𝑇 indicates the natural death rate of CD4+T-cells, r represents their rate of growth (therefore, 

𝑟 > µ𝑇 generally), and 𝑈𝑐𝑚𝑎𝑥 defines their carrying capacity. Given as a loss term for both healthy 

cells and virus since they are both lost by binding to one another and as the source term for infected 

cells, the parameter 𝑘1 describes the rate of infection of T-cells with free virus. Since both healthy 

cells and viruses are lost via attaching to one another, issued as a loss term for both, while also 

serving as the parameter for infected cells. The ratio 
𝑘1
′

𝑘1
 represents the proportion of T-cells that 

ever become actively infected. 𝑘1
′  indicates the rate at which infected cells become actively 

infected. To reflect the presumption that we initially do-not know whether the cells die normally 

or via bursting, µ𝐼 is a general word for the death of infected cells. To describe the dynamics of 

healthy cells, a time delay is added to the system.  

         Only healthy cells that the virus infected τ time units ago(i.e.,at time(t−τ) become infectious 

at time t. Thus, the incidence term of healthy cells is changed from 𝑘1
′𝑈𝑐(𝑡)𝐹𝑣(𝑡) to 

𝑘1
′𝑈𝑐(𝑡 − 𝜏)𝐹𝑣(𝑡 −  𝜏).  

          Additionally, the catalytic mortality rate for virus particles is µ𝑏.Each lysing proved to 

possess N viral particles, hence this term is multiplied by N to indicate the source of free virus 

(assuming a one-time initial infection). Finally, µ𝑉 is the viral loss rate. 

 

3. Linearization Technique 

        (i)To find the equilibrium points        

              Determining the values of equilibrium is required in order to adequately comprehend the 

dynamics of the three-component model. A system’s equilibrium point is a stable solution, (2.1) 

indicating that if the system starts out at that value, it will stay there forever. In other words, the 

inhabitants are static, and as a result, each population’s rate of change is zero. Two steady states 

exist for the system: 

 

(i) The uninfected steady state 𝐸0 = ([𝑈𝑐0(𝑡)]
𝛼
, 0,0), where [𝑈𝑐0(𝑡)]

𝛼
 is given by 

[𝑈𝑐0(𝑡)]
𝛼
= [

𝑟 − 𝜇𝑇[(𝑟 − 𝜇𝑇)
2 + 4𝑟𝑠𝑈𝑐

−1
𝑚𝑎𝑥

]
1/2

2𝑟𝑈𝑐𝑚𝑎𝑥−1
] [𝐹1(𝛼)]. 

(ii) The positively infected steady state [�̅�(𝑡)]𝛼 = ([𝑈𝑐̅̅ ̅(𝑡)]
𝛼, [𝐼�̅�(𝑡)]

𝛼, [𝐹�̅�(𝑡)]
𝛼), where 

[𝑈𝑐̅̅ ̅(𝑡)]
𝛼, [𝐼�̅�(𝑡)]

𝛼, [𝐹�̅�(𝑡)]
𝛼   are given by 

[𝑈𝑐̅̅ ̅(𝑡)]
𝛼 = (

𝜇𝐹𝑣𝜇𝐼𝑐
𝑘1
′𝑁𝜇𝑏 − 𝑘1𝜇𝐼𝑐

) [𝐹1(𝛼)],  [𝐼�̅�(𝑡)]
𝛼 = (

𝑘1
′𝑈𝑐𝐹𝑣̅̅ ̅̅ ̅̅

𝜇𝐼𝑐
) [𝐹2(𝛼)], 

[𝐹�̅�(𝑡)]
𝛼 = (

𝜇𝐼𝑐 [(𝑠 + (𝑟 − 𝜇𝑈𝑐)𝑈𝑐̅̅ ̅)𝑈𝑐𝑚𝑎𝑥 − 𝑟𝑈𝑐̅̅ ̅
2
]

𝑈𝑐̅̅ ̅[𝑘1
′𝑟𝑈𝑐̅̅ ̅ + 𝑘1𝜇1𝑈𝑐𝑚𝑎𝑥]

) [𝐹3(𝛼)]. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

2687 
https://internationalpubls.com 

Using the parameter values in Table, we find the three fixed(or) equilibrium points of the system (2.1) 

are, 

(i) Initially: [𝐸0]
𝛼 = ([𝑈𝑐0]

𝛼, 0, 0) = (1000, 0, 0), 

(ii) For N=500; 

[�̅�(𝑡)]𝛼 = ([𝑈𝑐̅̅ ̅(𝑡)]
𝛼, [𝐼�̅�(𝑡)]

𝛼, [𝐹�̅�(𝑡)]
𝛼) = (260.7[𝐹1(𝛼)], 35.5[𝐹2(𝛼)], 1768.2[𝐹3(𝛼)]) 

(iii) For N=1000; 

 

[�̅�(𝑡)]𝛼 = ([𝑈𝑐̅̅ ̅(𝑡)]
𝛼, [𝐼�̅�(𝑡)]

𝛼, [𝐹�̅�(𝑡)]
𝛼) = (130.2[𝐹1(𝛼)], 34.9[𝐹2(𝛼)], 3480.1[𝐹3(𝛼)]) 

iii)  To find the Jacobian matrix at the equilibrium points 

            The nonlinear fuzzy system (2.1) can be written as 

{
 
 
 
 

 
 
 
 𝑑

𝑑𝑡
([𝑈𝑐(𝑡)]

𝛼 − [𝑈𝑐̅̅ ̅(𝑡)]
𝛼) = 𝑀([𝑈𝑐(𝑡)]

𝛼 − [𝑈𝑐̅̅ ̅(𝑡)]
𝛼) −

𝑟𝑈𝑐̅̅ ̅
𝛼

𝑈𝑐𝑚𝑎𝑥
([𝐼𝑐(𝑡)]

𝛼 − [𝐼�̅�(𝑡)]
𝛼)

−𝑘1[𝑈𝑐̅̅ ̅(𝑡)]
𝛼([𝐹𝑣(𝑡)]

𝛼 − [𝐹�̅�(𝑡)]
𝛼),

𝑑

𝑑𝑡
([𝐼𝑐(𝑡)]

𝛼 − [𝐼�̅�(𝑡)]
𝛼) = 𝑘1

′ [𝐹�̅�(𝑡)]
𝛼([𝑈𝑐(𝑡 − 𝜏)]

𝛼 − [𝑈𝑐̅̅ ̅(𝑡)]
𝛼) − 𝜇1([𝐼𝑐(𝑡)]

𝛼 − [𝐼�̅�(𝑡)]
𝛼)

+𝑘1
′ [𝑈𝑐̅̅ ̅(𝑡)]

𝛼([𝐹𝑣(𝑡 − 𝜏)]
𝛼 − [𝐹�̅�(𝑡)]

𝛼),

𝑑

𝑑𝑡
([𝐹𝑣(𝑡)]

𝛼 − [𝐹�̅�(𝑡)]
𝛼) = −𝑘1[𝐹�̅�(𝑡)]

𝛼([𝑈𝑐(𝑡)]
𝛼 − [𝑈𝑐̅̅ ̅(𝑡)]

𝛼) + 𝑁𝜇𝑏([𝐼𝑐(𝑡)]
𝛼 − [𝐼�̅�(𝑡)]

𝛼)

−(𝑘1[𝑈𝑐̅̅ ̅(𝑡)]
𝛼 + 𝜇𝐹𝑣)([𝐹𝑣(𝑡)]

𝛼 − [𝐹�̅�(𝑡)]
𝛼).

            (3.1) 

 

Where 

𝑀 = µ𝑇
𝛼 + (

𝑟(2𝑈𝑐̅̅ ̅(𝑡) + 𝐼�̅�(𝑡))

𝑈𝑐𝑚𝑎𝑥
)

𝛼

+ (𝑘1𝐹𝑣(𝑡) − 𝑟)
𝛼 

System(3.1)is a linearized fuzzy system. Then, using a fuzzy matrix, we express the system 

(3.1) as follows.: 

[ 𝑋′(𝑡)]𝛼 = 𝐴1[𝑋(𝑡)]
𝛼 + 𝐴2[𝑋(𝑡 − 𝜏)]

𝛼                                                                             (3.2) 

where [𝑋(𝑡)]𝛼 = [

 [𝑤1(𝑡)]
𝛼

 [𝑤2(𝑡)]
𝛼

 [𝑤2(𝑡)]
𝛼

] and 𝐴1, 𝐴2 are 3x3fuzzy matrices are given by, 

𝐴1 = [

−𝑀
−𝑟𝑈𝑐̅̅̅̅

𝑈𝑐𝑚𝑎𝑥
−𝑘1𝑈𝑐̅̅ ̅

0 −𝜇1 0

−𝑘1𝐹�̅� 𝑁𝜇𝑏 −(𝑘1𝑈𝑐 + 𝜇𝐹𝑣̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅)

] [

𝐹1(𝛼)
𝐹2(𝛼)
𝐹3(𝛼)

] and  

 

𝐴2 = [
0 0 0
𝑘1𝐹�̅� 0 𝑘1′𝑈𝑐̅̅ ̅

0 0 0

] [

𝐹1(𝛼)
𝐹2(𝛼)

𝐹3(𝛼)
]                                                                                       (3.3) 

Using the parameter values from the table and the equilibrium point for N=500, 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

2688 
https://internationalpubls.com 

[�̅�]𝛼 = (260.7[𝐹1(𝛼)], 35.5[𝐹2(𝛼)], 1768.2[𝐹3(𝛼)]) 

And 𝜏 = 1, (3.1) can be written as, 

{

 [𝑤′1(𝑡)]
𝛼 = −0.0436[𝑤1(𝑡)]

𝛼 − 0.0052[𝑤2(𝑡)]
𝛼 − 0.00626[𝑤3(𝑡)]

𝛼,

 [𝑤′2(𝑡)]
𝛼 = 0.0354[𝑤1(𝑡 − 1)]

𝛼 − 0.26[𝑤2(𝑡)]
𝛼 + 0.0052[𝑤3(𝑡 − 1)]

𝛼,

 [𝑤′3(𝑡)]
𝛼 = −0.0424[𝑤1(𝑡)]

𝛼 + 120[𝑤2(𝑡)]
𝛼 − 2.24063[𝑤3(𝑡)]

𝛼.

                (3.4) 

At the equilibrium point for N=1000, 

[�̅�]𝛼 = (130.2[𝐹1(𝛼)], 34.9[𝐹2(𝛼)], 3480.1[𝐹3(𝛼)]). 

And 𝜏 = 1, (3.1) can be written as, 

{

 [𝑤′1(𝑡)]
𝛼 = −0.07942[𝑤1(𝑡)]

𝛼 − 0.0026[𝑤2(𝑡)]
𝛼 − 0.003124[𝑤3(𝑡)]

𝛼,

 [𝑤′2(𝑡)]
𝛼 = 0.0696[𝑤1(𝑡 − 1)]

𝛼 − 0.26[𝑤2(𝑡)]
𝛼 + 0.0026[𝑤3(𝑡 − 1)]

𝛼,

 [𝑤′3(𝑡)]
𝛼 = −0.0835[𝑤1(𝑡)]

𝛼 + 240[𝑤2(𝑡)]
𝛼 − 2.20312[𝑤3(𝑡)]

𝛼.

                (3.5) 

 [𝑤1(𝑡)]
𝛼  =  [1001 − 𝑒𝑡][𝐹1(𝛼)], 

 [𝑤2(𝑡)]
𝛼 = [0][𝐹2(𝛼)],  

[𝑤3(𝑡)]
𝛼  =  [1.001 − 𝑒𝑡][𝐹3(𝛼)],   𝑓𝑜𝑟 𝑡 ∈ [−1,0]. 

With initial functions. 

The analytical solution is  

[𝑋(𝑡)]𝛼  =  [𝛷(𝑡)]𝛼𝐶 + [𝛷(𝑡)]𝛼∫ [𝛷−1(𝑠)]𝛼𝐵(𝑠, 𝛼)𝑑𝑠, 𝑡 ∈ 𝐼.                                               (3.6)
𝑡

𝑡0

 

The fundamental matrix of the system(3.5) is given by, 

[𝛷(𝑡)]𝛼 = [

𝑎1𝑒
∝𝑡 𝑏1𝑒

∝𝑡 𝑐1𝑒
∝𝑡

𝑎2𝑒
∝𝑡 𝑏2𝑒

∝𝑡 𝑐2𝑒
∝𝑡

𝑎3𝑒
∝𝑡 𝑏3𝑒

∝𝑡 𝑐3𝑒
∝𝑡

] [

𝐹1(𝛼)

𝐹2(𝛼)
𝐹3(𝛼)

]                                                                              (3.7) 

where 𝛼 = −0.0434876, 𝛽 = −0.25, 𝛾 = −2.40641, and 

 

{
a1 = 0.999839, a2 = 0, a3 = −0.0179587,

b1 = 0.0293406,   b2 = 0.0178856, b3 = 0.9994090,

d1 = 0.000264938, d2 = 0, d3 = 0.9999960
 

Let, 𝑏 = [
0

35.4406 − (0.0406𝑒𝑡−1)
0

] and Φ−1(𝑡)𝐵 =
1

𝑀
[

𝑚1𝑒
−𝛼𝑡

𝑚2𝑒
−𝛽𝑡

𝑚3𝑒
−𝛾𝑡

] where,  

𝑀 = 𝑎1ℎ2 − 𝑏1ℎ1 + 𝑑1ℎ3,  

𝑚1 = 𝑛1ℎ1 − 𝑛2ℎ4 + 𝑛3ℎ7, 

 𝑚2 = −𝑛1ℎ2 + 𝑛2ℎ5 − 𝑛3ℎ8,  

𝑚3 = 𝑛1ℎ3 − 𝑛2ℎ6 + 𝑛3ℎ9,  

𝑛1 = 0,          𝑛2 = 35.4406 − (0.0406𝑒
𝑡−1)), 𝑛3 = 0 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

2689 
https://internationalpubls.com 

and  

{

h1 = γ3a2 − γ2a3, h2 = γ3b2 − γ2b3,   h3 = a2b3 − b2a3,

h4 = γ
1
b3 − γ3b1, h5 = γ

3
a1 − γ1a3,   h6 = a1b3 − a3b1,

h7 = γ1b2 − γ2b1, h8 = γ2a1 − γ1a2,   h9 = a1b2 − b1a2.

  

 

The linear system (3.5)’s analytical solution is, 

𝑥(𝑡) = 𝐶1𝑎1𝑒
𝛼𝑡 + 𝐶2𝑏1𝑒

𝛽𝑡 + 𝐶3𝑑1𝑒
𝛾𝑡 + 1/𝑀 [

𝑎1𝑚1

𝛼
(𝑒𝛼 − 1) +

𝑏1𝑚1

𝛽
(𝑒𝛽 − 1) +

𝑑1𝑚1

𝛾
(𝑒𝛾 − 1)

] 

𝑦(𝑡) = 𝐶1𝑎2𝑒
𝛼𝑡 + 𝐶2𝑏2𝑒

𝛽𝑡 + 𝐶3𝑑2𝑒
𝛾𝑡 + 1/𝑀 [

𝑎2𝑚1

𝛼
(𝑒𝛼 − 1) +

𝑏2𝑚1

𝛽
(𝑒𝛽 − 1) +

𝑑2𝑚1

𝛾
(𝑒𝛾 − 1)

]’ 

𝑧(𝑡) = 𝐶1𝑎3𝑒
𝛼𝑡 + 𝐶2𝑏3𝑒

𝛽𝑡 + 𝐶3𝑑3𝑒
𝛾𝑡 + 1/𝑀 [

𝑎3𝑚1

𝛼
(𝑒𝛼 − 1) +

𝑏3𝑚1

𝛽
(𝑒𝛽 − 1) +

𝑑3𝑚1

𝛾
(𝑒𝛾 − 1)

]’ 

where 𝐶1 = 10001.4, 𝐶2 = 0, 𝐶3 = 179.793. 

4.  Stability Analysis 

       The linearized system’s characteristic equation, (3.2), is generated by, 

∆(𝜆)  =  |𝜆𝐼 + 𝐴1 − 𝑒
−𝜆𝜏𝐴2| = 0 

That is,  

𝜆3 + 𝑐1𝜆2 + 𝑐2𝜆 + 𝑐3𝑒
−𝜆𝜏 + 𝑐4𝜆𝑒

−𝜆𝜏 + 𝑐5  =  0,                                                              (4.1) 

Compared to the ODE model, it is difficult to find the eigenvalues of (4.1), since it is a transcen 

dental equation and it has infinitely many eigenvalues. Let 𝜆 = 𝜂(𝜏) + 𝑖𝜔(𝜏)(𝜔 > 0), be the 

eigenvalue off the characteristic equation (4.1), where 𝜂(𝜏) and 𝜔(𝜏) depend only on the delay 

𝜏. Suppose 𝑖𝜔(𝜔 > 0) is a root of (4.1) if and only if  ℎ(𝑧) = 𝑧3 + µ𝑧2  + 𝑣𝑍 + 𝜌 = 0, where 

𝑧 = 𝜔2, µ = 𝑐1
2 − 2𝑐2, 𝑣 = 𝑐2

2 − 2𝑐1𝑐5 − 𝑐4
2, 𝜌 = 𝑐5

2 − 𝑐3
2. The following proposition 

demonstrates that the steady state of the delay model (2.1) is asymptotically stable for all delay 

values if the parameter values in Table 1 satisfy the specified conditions.  

            4.1. Proposition  

             If (i)𝑐1 > 0, 𝑐3 + 𝑐5 > 0, 𝑐1(𝑐2 + 𝑐4)  −  (𝑐3 + 𝑐5)  >  0, (ii) 𝜌 ≥  0, and 𝜈 >  0, then 

the fuzzy delay model’s infected steady state (2.1) is completely stable. that is, ”E” is stable 

asymptotically across all 𝜏 ≥  0.  

Case (i): For N = 500, then𝑐1 = 3.07, 𝑐2 = 1.70173, 𝑐 3 = −0.2256, 𝑐4 = −0.6255, 𝑐5 =

0.2525, 𝑐3 + 𝑐5 = 0.027 > 0, 𝑐1(𝑐2 + 𝑐4) − (𝑐3 + 𝑐5) = 3.27703 > 0, and ρ = 0.01286 ≥ 0,v 

= 0.95428 > 0. The infected steady state 𝐸 =  (260.7[𝐹1(𝛼)], 35.5[𝐹2(𝛼)],1768.2[𝐹3(𝛼)]) is 

asymptotically stable, as demonstrated by the numerical simulations and all of the parameter 

values in Table I.  

Case (ii): There will be a significant increase in the number of infected cells and a decrease in 

the number of uninfected CD4+ T-cells as the N value is raised, but the steady state will remain 

stable. If N equals 1000, then 𝑐1 = 2.7422, 𝑐2 = 0.8360, 𝑐3 = 0.00299, 𝑐4 = −0.62477, 𝑐5 =

0.04955, 𝑐3 + 𝑐5 = 0.0525535 > 0, 𝑐1(𝑐2 + 𝑐4) − (𝑐3 + 𝑐5) = 0.526945 > 0, and 𝜌 =



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

2690 
https://internationalpubls.com 

 0.00244 ≥  0, 𝑣 =  0.036806 >  0. As a result, the infected steady state 𝐸 =

(130.2[𝐹1(𝛼)], 34.9[𝐹2(𝛼)], 3480.1[𝐹3(𝛼)])is asymptotically stable based on all of the 

parameter values in Table 1. Therefore, for any delay 𝜏 ≥  0, the infected steady state E is 

statistically stable. If proposition 4.1’s conditions (i) and (ii) are not met, the steady state 

stability is dependent on the delay value, and the delay may even cause oscillations. When 𝜏 

passes through a critical value 𝜏𝑗, a Hopf-bifurcation [9] occurs. 

𝜏𝑗 =
1

𝜔0
𝑎𝑟𝑐 𝑐𝑜𝑠 (

𝑐4𝜔0
4 + (𝑐1𝑐3 − 𝑐2𝑐4)𝜔0

2 − 𝑐3𝑐5
𝑐3
2 + 𝑐4

2𝜔0
2 ) +

2𝑗𝜋

𝜔0
, 𝑗 = 0,1,2, . .. 

Such that 𝜂(𝑡0)=0, 𝜔(𝜏0) = 𝜔0. The transversality condition is: 

𝑑

𝑑𝜏
𝑅𝑒(𝜆(𝜏))|𝜏=𝜏0 =

𝑑

𝑑𝜏
𝑅𝑒(𝜂(𝜏))|𝜏=𝜏0 > 0. 

4.2. Proposition  

      If (i)𝑐1 > 0, 𝑐3 + 𝑐5 > 0, 𝑐1(𝑐2 + 𝑐4) − (𝑐3 + 𝑐5) > 0 (ii) 𝜌 <  0, (iii) 𝜌 ≥  0, and 𝑣0, are 

satisfied, then the infected steady state of the fuzzy delay model (2.1) is asymptotically stable 

when 𝜏, 𝜏0, and unstable when 𝜏 > 𝜏0, where 

𝜏0 =
1

𝜔0
𝑎𝑟𝑐 𝑐𝑜𝑠 (

𝑐4𝜔0
4 + (𝑐1𝑐3 − 𝑐2𝑐4)𝜔0

2 − 𝑐3𝑐5
𝑐3
2 + 𝑐4

2𝜔0
2 ). 

When 𝜏 = 𝜏0, a Hopf bifurcation occurs; that is, a family of periodic solutions bifurcates from 

E as 𝜏 passes through the critical value 𝜏0. 

4.3. Proposition  

The characteristic equation of the linearized system is given by |𝐴 = 𝜆𝐼| = 𝜆3  

+𝑐2𝜆
2 + 𝑐1𝜆 + 𝑐0  =  0, where for 𝑁 = 500, 𝑐2 = 2.7099, 𝑐1 = 0.741621, 𝑐0  =

 0.0272083, for 𝑁 = 1000, 𝑐2 = 2.74225, 𝑐1 = 0.836078, 𝑐0  =  0.049557. The eigenvalues 

of the matrix A are for 𝑁 =  500, 𝛼 =  −0.0434876, 𝛽 =  −0.26, 𝛾 =  −2.40641, for 𝑁 =

 1000, 𝛼 =  −0.0793157, 𝛽 =  −0.26, 𝛾 =  −2.40324. If and only if all of the eigenvalues 

have negative real parts, the system is said to be stable. If any one of the eigenvalues has a 

positive real part, the system is said to be unstable. We have three eigenvalues that are linearly 

independent. Negative real parts are present in all eigenvalues. Our system is therefore stable.  

5.  Numerical Simulation of SIR model 

              In this section, we are using the RKM-5. We finding value of 𝑈𝑐(𝑡), 𝐼𝑐(𝑡) and 𝐹𝑣(𝑡) 

at ℎ = 0.1 for the best approximation. For 0 ≤  𝛼 ≤  1. To evaluate 𝑈𝑐(𝑡), 𝐼𝑐(𝑡) and 𝐹𝑣(𝑡), 

consider, 

[𝑈𝑐(𝑡 + 1)]
𝛼 = [𝑈𝑐(𝑡)]

𝛼  + 
1

90
(7[𝐾1]

𝛼 + 32[𝐾3]
𝛼 + 12[𝐾4]

𝛼 + 32[𝐾5]
𝛼 + 7[𝐾6]

𝛼) , 

            [𝐼𝑐(𝑡 + 1)]
𝛼 = [𝐼𝑐(𝑡)]

𝛼  +  
1

90
(7[𝐾1]

𝛼 + 32[𝐾3]
𝛼 + 12[𝐾4]

𝛼 + 32[𝐾5]
𝛼 + 7[𝐾6]

𝛼) ,      (5.1)             

[𝐹𝑉(𝑡 + 1)]
𝛼 = [𝐹𝑣(𝑡)]

𝛼  +  
1

90
(7[𝐾1]

𝛼 + 32[𝐾3]
𝛼 + 12[𝐾4]

𝛼 + 32[𝐾5]
𝛼 + 7[𝐾6]

𝛼) ,  

In estimate (5.1), for 1 ≤  𝑝 ≤  6 and 0 ≤  𝛼 ≤  1, we use  

[𝐾𝑝]
𝛼
= [𝐾𝑝(𝑡; 𝛼), 𝐾𝑝(𝑡; 𝛼)], [𝐿𝑝]

𝛼
= [𝐿𝑝(𝑡; 𝛼), 𝐿𝑝(𝑡; 𝛼)] [𝑀𝑝]

𝛼
 =  [𝑀𝑝(𝑡; 𝛼),𝑀𝑝(𝑡; 𝛼)].  



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

2691 
https://internationalpubls.com 

For 0 ≤  𝑡 ≤ 𝑛, 𝑛 = 1,2,3, . . ., and for  𝑞 = 𝑡 + 1, 𝑡 = 0,1,2,3, . ..  

[𝐾(𝑞)]𝛼 = [𝐾𝑞(𝑡; 𝛼), 𝐾𝑞(𝑡; 𝛼)], [𝐿(𝑞)]
𝛼 = [𝐿𝑝(𝑡; 𝛼), 𝐿𝑝(𝑡; 𝛼)] [𝑀(𝑞)]

𝛼  =  [𝑀𝑝(𝑡; 𝛼),𝑀𝑝(𝑡; 𝛼)]. 

Where, 

[𝐾1]
𝛼 = ℎ(𝑠 − µ𝑡[𝑈𝑐(𝑡)]

𝛼 + 𝑟[𝑈𝑐(𝑡)]
𝛼(1 − ([𝑈𝑐(𝑡)]

𝛼 + [𝐼𝑐(𝑡)]
𝛼)/𝑈𝑐𝑚𝑎𝑥) − 𝑘1[𝐹𝑣(𝑡)]

𝛼[𝑈𝑐(𝑡)]
𝛼), 

[𝐿1]
𝛼  =  ℎ(𝑘1

′ [𝐹𝑣(𝑡 − 𝜏)]
𝛼[𝑈𝑐(𝑡 − 𝜏)]

𝛼 − µ𝐼[𝐼𝑐(𝑡)]
𝛼),  

[𝑀1]
𝛼  =  ℎ(𝑁µ𝑏[𝐼𝑐(𝑡)]

𝛼 − 𝑘1[𝐹𝑣(𝑡)]
𝛼[𝑈𝑐(𝑡)]

𝛼 − µ𝑣[𝐹𝑣(𝑡)]
𝛼). 

[𝐾2]
𝛼  =  ℎ(𝑠 − µ𝑡([𝑈𝑐(𝑡)]

𝛼 + 1/2[𝐾1]
𝛼) + 𝑟([𝑈𝑐(𝑡)]

𝛼 + 1/2[𝐾1]
𝛼) 

                   (1 − (([𝑈𝑐(𝑡)]
𝛼 + 1/2[𝐾1]

𝛼) + ([𝐼𝑐(𝑡)]
𝛼 + 1/2[𝐿1]

𝛼))/𝑈𝑐𝑚𝑎𝑥) − 𝑘1([𝐹𝑣(𝑡)]
𝛼

+ 1/2[𝑀1]
𝛼)([𝑈𝑐(𝑡)]

𝛼 + 1/2[𝐾1]
𝛼)), 

[𝐿2]
𝛼  =  ℎ(𝑘1

′ ([𝐹𝑣(𝑡 − 𝜏)]
𝛼 + 1/2[𝑀1]

𝛼)([𝑈𝑐(𝑡 − 𝜏)]
𝛼 + 1/2[𝐾1]

𝛼) − µ𝐼 ([𝐼𝑐(𝑡)]
𝛼 + 1/2[𝐿1]

𝛼)),  

[𝑀2]
𝛼  =  ℎ(𝑁µ𝑏([𝐼𝑐(𝑡)]

𝛼 + 1/2[𝐿1]
𝛼) − 𝑘1([𝐹𝑣(𝑡)]

𝛼 + 1/2[𝑀1]
𝛼)([𝑈𝑐(𝑡)]

𝛼 + 1/2[𝐾1]
𝛼)

− µ𝑣([𝐹𝑣(𝑡)]
𝛼 + 1/2[𝑀1]

𝛼)). 

 [𝐾3]
𝛼  =  ℎ(𝑠 − µ𝑡([𝑈𝑐(𝑡)]

𝛼 + 3/16[𝐾1]
𝛼 + 1/16[𝐾2]

𝛼) + 𝑟([𝑈𝑐(𝑡)]
𝛼 + 3/16[𝐾1]

𝛼 +

                  1/16[𝐾2]
𝛼)(1 − (([𝑈𝑐(𝑡)]

𝛼 + 3/16[𝐾1]
𝛼 + 1/16[𝐾2]

𝛼) + ([𝐼𝑐(𝑡)]
𝛼 + 3/16[𝐿1]

𝛼 +

                   1/16[𝐿2]
𝛼))/𝑈𝑐𝑚𝑎𝑥) − 𝑘1([𝐹𝑣(𝑡)]

𝛼 + 3/16[𝑀1]
𝛼 + 1/16[𝑀2]

𝛼)([𝑈𝑐(𝑡)]
𝛼 +

                    3/16[𝐾1]
𝛼 + 1/16[𝐾2]

𝛼)), 

[𝐿3]
𝛼  =  ℎ(𝑘1

′ ([𝐹𝑣(𝑡 − 𝜏)]
𝛼 + 3/16[𝑀1]

𝛼 + 1/16[𝑀2]
𝛼)([𝑈𝑐(𝑡 − 𝜏)]

𝛼 + 3/16[𝐾1]
𝛼   +

                     1/16[𝐾2]
𝛼) − µ𝐼𝑐 ([𝐼𝑐(𝑡)]

𝛼 + 3/16[𝐿1]
𝛼 + 1/16[𝐿2]

𝛼)),  

[𝑀3]𝛼 =  ℎ(𝑁µ𝑏([𝐼𝑐(𝑡)]𝛼 + 3/16[𝐿1]𝛼 + 1/16[𝐿2]𝛼)  − 𝑘1([𝐹𝑣(𝑡)]𝛼 +  3/16[𝑀1]𝛼 

+  1/16[𝑀2]𝛼)([𝑈𝑐(𝑡)]𝛼 +  3/16[𝐾1]𝛼 +  1/16[𝐾2]𝛼) − µ𝑣([𝐹𝑣(𝑡)]𝛼 

+  3/16[𝑀1]𝛼 +  1/16[𝑀2]𝛼)).  

[𝐾4]
𝛼  =  ℎ(𝑠 − µ𝑡([𝑈𝑐(𝑡)]

𝛼 + 1/2[𝐾3]
𝛼) + 𝑟([𝑈𝑐(𝑡)]

𝛼 + 1/2[𝐾3]
𝛼)(1 − (([𝑈𝑐(𝑡)]

𝛼 +

                   1/2[𝐾3]
𝛼) + ([𝐼𝑐(𝑡)]

𝛼 + 1/2[𝐿3]
𝛼))/𝑈𝑐𝑚𝑎𝑥) − 𝑘1([𝐹𝑣(𝑡)]

𝛼 + 1/2[𝑀3]
𝛼)([𝑈𝑐(𝑡)]

𝛼 +

                    1/2[𝐾3]
𝛼)),  

[𝐿4]
𝛼  =  ℎ(𝑘1

′ ([𝐹𝑣(𝑡 − 𝜏)]
𝛼 + 1/2[𝑀3]

𝛼)([𝑈𝑐(𝑡 − 𝜏)]
𝛼 + 1/2[𝐾3]

𝛼) − µ𝐼𝑐 ([𝐼(𝑡)]
𝛼 + 1/2[𝐿3]

𝛼)),  

[𝑀4]
𝛼  =  ℎ(𝑁µ𝑏([𝐼𝑐(𝑡)]

𝛼 + 1/2[𝐿3]
𝛼) − 𝑘1([𝐹𝑣(𝑡)]

𝛼 + 1/2[𝑀3]
𝛼)([𝑈𝑐(𝑡)]

𝛼 + 1/2[𝐾3]
𝛼) −

                   µ𝑣([𝐹𝑣(𝑡)]
𝛼 + 1/2[𝑀3]

𝛼)).  

[𝐾5]
𝛼  =  ℎ(𝑠 − µ𝑡([𝑈𝑐(𝑡)]

𝛼 − 3/16[𝐾2]
𝛼 + 6/16[𝐾3]

𝛼 + 9/16[𝐾4]
𝛼)

+ 𝑟([𝑈𝑐(𝑡)]
𝛼 − 3/16[𝐾2]

𝛼 + 6/16[𝐾3]
𝛼 + 9/16[𝐾4]

𝛼)(1

− (([𝑈𝑐(𝑡)]
𝛼 − 3/16[𝐾2]

𝛼 + 6/16[𝐾3]
𝛼 + 9/16[𝐾4]

𝛼)

+ ([𝐼𝑐(𝑡)]
𝛼 − 3/16[𝐿2]

𝛼 + 6/16[𝐿3]
𝛼 + 9/16[𝐿4]

𝛼))/𝑈𝑐𝑚𝑎𝑥) − 𝑘1([𝐹𝑣(𝑡)]
𝛼

− 3/16[𝑀2]
𝛼 + 6/16[𝑀3]

𝛼 + 9/16[𝑀4]
𝛼)([𝑈𝑐(𝑡)]

𝛼 − 3/16[𝐾2]
𝛼 + 6/16[𝐾3]

𝛼

+ 9/16[𝐾4]
𝛼)), 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

2692 
https://internationalpubls.com 

[𝐿5]
𝛼 = ℎ(𝑘1

′ ([𝐹𝑣(𝑡 − 𝜏)]
𝛼 − 3/16[𝑀2]

𝛼 + 6/16[𝑀3]
𝛼 + 9/16[𝑀4]

𝛼)([𝑈𝑐(𝑡 − 𝜏)]
𝛼 − 3/16[𝐾2]

𝛼

+ 6/16[𝐾3]
𝛼 + 9/16[𝐾4]

𝛼)

− µ𝐼𝑐([𝐼𝑐(𝑡)]
𝛼 − 3/16[𝐿2]

𝛼 + 6/16[𝐿3]
𝛼 + 9/16[𝐿4]

𝛼)), 

[𝑀5]
𝛼 =  ℎ(𝑁µ𝑏([𝐼𝑐(𝑡)]

𝛼 − 3/16[𝐿2]
𝛼 + 6/16[𝐿3]

𝛼 + 9/16[𝐿4]
𝛼)

− 𝑘1([𝐹𝑣(𝑡)]
𝛼 − 3/16[𝑀2]

𝛼 + 6/16[𝑀3]
𝛼 + 9/16[𝑀4]

𝛼)([𝑈𝑐(𝑡)]
𝛼 − 3/16[𝐾2]

𝛼

+ 6/16[𝐾3]
𝛼 + 9/16[𝐾4]

𝛼) − µ𝑣([𝐹𝑣(𝑡)]
𝛼 − 3/16[𝑀2]

𝛼 + 6/16[𝑀3]𝛼

+ 9/16[𝑀4]
𝛼)). 

[𝐾6]
𝛼  =  ℎ(𝑠 − µ𝑡([𝑈𝑐(𝑡)]

𝛼 + 1/7[𝐾1]
𝛼 + 4/7[𝐾2]

𝛼 + 6/7[𝐾3]
𝛼 − 12/7[𝐾4]

𝛼 + 8/7[𝐾5]
𝛼)  

+ 𝑟([𝑈𝑐(𝑡)]
𝛼 + 1/7[𝐾1]

𝛼 + 4/7[𝐾2]
𝛼 + 6/7[𝐾3]

𝛼 − 12/7[𝐾4]
𝛼 + 8/7[𝐾5]

𝛼) (1 

− (([𝑈𝑐(𝑡)]
𝛼 + 1/7[𝐾1]

𝛼 + 4/7[𝐾2]
𝛼 + 6/7[𝐾3]

𝛼 

                             −12/7[𝐾4]
𝛼 + 8/7[𝐾5]

𝛼)  + ([𝐼𝑐(𝑡)]
𝛼  + 1/7[𝐿1]

𝛼 + 4/7[𝐿2]
𝛼 + 6/7[𝐿3]

𝛼 −

                              12/7[𝐿4]
𝛼 + 8/7[𝐿5]

𝛼))/𝑈 𝑐𝑚𝑎𝑥)  − 𝑘1([𝐹𝑣(𝑡)]
𝛼  +  1/7[𝑀1]

𝛼  +  4/7[𝑀2]
𝛼 +

                             6/7[𝑀3]
𝛼 − 12/7[𝑀4]

𝛼  + 8/7[𝑀5]
𝛼) ([𝑈𝑐(𝑡)]

𝛼  +  1/7[𝐾1]
𝛼 +  4/7[𝐾2]

𝛼  +

                             6/7[𝐾3]
𝛼  −  12/7[𝐾4]

𝛼  +  8/7[𝐾5]
𝛼)),  

[𝐿6]
𝛼  =  ℎ(𝑘1

′ ([𝐹𝑣(𝑡 − 𝜏)]
𝛼 + 1/7[𝑀1]

𝛼 + 4/7[𝑀2]
𝛼 + 6/7[𝑀3]

𝛼 − 12/7[𝑀4]
𝛼

+ 8/7[𝑀5]
𝛼)([𝑈𝑐(𝑡 − 𝜏)]

𝛼 + 1/7[𝐾1]
𝛼 + 4/7[𝐾2]

𝛼 + 6/7[𝐾3]
𝛼 − 12/7[𝐾4]𝛼

+ 8/7[𝐾5]
𝛼) − µ𝐼𝑐  ([𝐼𝑐(𝑡)]

𝛼 + 1/7[𝐿1]
𝛼 + 4/7[𝐿2]

𝛼 + 6/7[𝐿3]
𝛼 − 12/7[𝐿4]

𝛼

+ 8/7[𝐿5]
𝛼)), 

[𝑀6]
𝛼 =  ℎ(𝑁µ𝑏([𝐼𝑐(𝑡)]

𝛼 − 3/16[𝐿2]
𝛼 + 6/16[𝐿3]

𝛼 + 9/16[𝐿4]
𝛼) − 𝑘1([𝐹𝑣(𝑡)]

𝛼 + 1/7[𝑀1]
𝛼

+ 4/7[𝑀2]
𝛼 + 6/7[𝑀3]

𝛼 − 12/7[𝑀4]
𝛼 + 8/7[𝑀5]

𝛼)([𝑈𝑐(𝑡)]
𝛼 + 1/7[𝐾1]

𝛼

+ 4/7[𝐾2]
𝛼 + 6/7[𝐾3]

𝛼 − 12/7[𝐾4]
𝛼 + 8/7[𝐾5]𝛼)  − µ𝑣([𝐹𝑣(𝑡)]𝛼 + 1/7[𝑀1]𝛼

+ 4/7[𝑀2]
𝛼 + 6/7[𝑀3]

𝛼 − 12/7[𝑀4]
𝛼 + 8/7[𝑀5]

𝛼)). 

 

 

 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

2693 
https://internationalpubls.com 

 

6. Result and Discussion 

                   First, we reconstructed Rebecca et al.’s[4] DDE model of HIV infection in CD4 + T 

cells, into a three system of linear equations. In order to extend the infection, we discover a limit 

on the number of viral particles released by each infectious cell. Under this constraint, the system 

is in an infected steady state with a positive equilibrium. We established sufficient parameters for 

the steady state stability of the infected virus by employing stability analysis. Numerical results 

demonstrated that the analysis was correct, and all of our stability conditions are met. Despite the 

fact that our estimation of the number of viral particles released from each infectious cell is lower 

than Perelson et al.’s [14], It has no effect on the existence or stability of the infected steady state. 

There is a solution to every nonlinear differential equation. Analytical solutions to the majority of 

nonlinear differential equations are extremely challenging to find. We must therefore use numerical 

approaches. It is simple to solve all nonlinear differential equations using numerical techniques. 

To solve the delay differential equation numerically, Runge-Kutta’s fifth-order method is the most 

effective approach. 

Refrences 

[1] J.C. Butcher, A history of Runge-Kutta methods, Applied Numerical Mathematics, 20 (1996), 

247-260. 

[2] Capistran, M.A., and Solis F.J., On the modeling of long-term HIV-I infection dynamics, 

Mathematical and computer modeling,(2009), 777-782 

[3] S.M. Ciupe, B.L. de Bivort and P.W. Nelson, Estimates of kinetics parameters from HIV 

patient data during primary infection through the eyes of three different models, Math. Biosci. 

in press. 

[4] R.V. Culshaw, S. Ruan, A delay differential equation model of HIV infection of CD4+ T-cells, 

Math. Biosci. 165 (2000), 27-39. 

[5] R.V. Culshaw, S. Ruan, G. Webb, A mathematical model of cell-to-cell spread of HIV-1 that 

includes a time delay, Journal of Mathematical Biology, 46 (2003), 425-444. 

[6] R.P. Duffin, R.H. Tullis, Mathematical models of the complete course of HIV infection and 

AIDS, Journal of Theoretical Medicine, Vol.4 (2002), pp. 215-221. 

[7] J. Hussain, R. Lalawmpuii, Modelling the dynamics of CD4+ T-cells with and without delay, 

MIZO Academy of Sciences, Vol 13, ISSN 0975-6175. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

2694 
https://internationalpubls.com 

[8] A.V.M. Herz, S.Bonhoeffer, R.M.Anderson, R.M.May and M.A.Nowak, Viral dynamics in 

vivo: Limitations on estimates of intracellular delay and virus delay. In.Proc.Natl.Acad.Sci 

(Medical Sciences), 7247-7251, USA, 1996. 

[9] B.D.Hasssard, N.D.Kazarinoff, Y.H.Wan, Theory and Applications of Hopf Bifurcation, Cam- 

bridge University, Cambridge,1981. 

[10] P.W. Nelson, A.S. Perelson, Mathematical analysis of delay differential equation models of 

HIV-1 infection, Math. Biosci. 179 (2002) 73. 

[11] M.A.Nowak, and R.M.May, Mathematical biology of HIV infections, Antigenic variation and 

diversity Threshold, Mathematical Biosciences, 106(1) pp.1-21,1991. 

[12] T.Muthukumar and T.Jayakuamr, Numerical Solution for Hybrid Fuzzy Differential Equations 

by Fifth Order Runge-Kutta Nystrom Method, Journal of Mathematical Sciences and 

Modelling, 2(1) (2019) 39-50. 

[13] T.Muthukumar, T.Jayakumar and D.Prasantha Bharathi, Numerical Solution of Fuzzy Delay 

Differential Equations by Runge-Kutta Fehlberge Method, GIS Science journal, (7)(12)2020, 

507-516. 

[14] A.S. Perelson, Modelling viral and immune system dynamics, Nature Reviews Immunology, 2 

(2002), 28-36. 

[15] L.F. Shampine, Numerical Solution of Ordinary Differential Equations, Chapman and Hall, 

New York, 1994.  

[16] L.F. Shampine, I. Gladwell, S. Thomson, Solving ODEs with Matlab, Cambridge University 

Press, UK, 2003. 

[17] R.W. Shonkwiler, J. Herod, Mathematical Biology, An introduction with Maple and Matlab, 

Second edition. 

[18] X. Song, S. Cheng, A Delay differential equation model of HIV infection of CD4+ T-cells, 

Journal of korean mathematical society, 42 (2005), No. 5, pp. 1071-1086. 

[19] L. Wang, M.Y. Li, Mathematical analysis of the global dynamics of a model for HIV infection 

of CD4+ T-cells, Math. Bio. 200 (2006), 44-57. 

[20] M. Yan, Z. Xiang, A Delay-Differential Equation Model of HIV infection of CD4+ T-cells 

with Cure Rate, International Mathematical Forum, Vol. 7 (2012), No. 30, 1475-1481. 

[21] S. Yuzbasi, A numerical approach to solve the model for HIV infection of CD4+ T-cells, 

Applied Mathematical Modelling 36 (2012), 5876-5890. 

[22] X. Zhou, X. Song, X. Shi, A differential equation model of HIV infection of CD4+ T-cells with 

cure rate, Journal of Mathematical Analysis, 342 (2008), 1342-1355.