# Dynamical Analysis of Posttreatment HIV-1 Infection Model

Authors: M. Pradeesh; A. Manivannan; S. Lakshmanan; F. A. Rihan; Prakash Mani
Journal: Complexity
 (2022)

---

## Abstract

This paper aims to explore the dynamic characteristics of the post treatment human immunodeficiency virus (HIV) type-1 model by proposing the theoretical frameworks. Distinct from the previous works, this study explores the effect of effector cells, loss of functional effector cells, and two types of anti-retroviral therapies such as reverse transcriptase inhibitors (RTIs) and protease inhibitors (PIs) and also the effect of intracellular time delay. Based on the Routh—Hurwitz criterion and eigenvalue analysis, the stability of the proposed HIV-1 model is analyzed. To reveal the significance of time delay, the Hopf-type bifurcation analysis is performed. The optimal control algorithm is designed by choosing the antiviral therapies such as RTI and PI as control parameters. Numerical simulations are performed to validate the effectiveness of the proposed theoretical frameworks.


---

## Body

## 1. Introduction

According to the World Health Organization (WHO), 680,000 individuals died from HIV-1-related diseases worldwide in 2020, ranging from 480,000 to 1.0 million, while 1.5 million people were newly infected with HIV-1. Antiretroviral therapy (ART) has drastically reduced the number of people infected with HIV during the 1990s, with 27.5 million (approximately) people undergoing treatment in 2020. Because of antiretroviral medication, the HIV infection rate has decreased by 49% over the last two decades, from 2000 to 2020. Despite the researcher’s valiant efforts in terms of treatment options and drugs, a cure for HIV-1 remains a pipe dream, necessitating lifelong treatment. Mathematical models have been demonstrated to be a useful tool for comprehending the dynamics of disease progression, identifying key determinants, and evaluating the efficacy of antiretroviral therapy.During the 1990s, HIV was thought to be a lethal disease, similar to other lentiviruses, because HIV remains within the host without causing symptoms and progresses to a chronic stage known as acquired immunodeficiency virus (AIDS), with a nearly ten-year delay between HIV and AIDS. In this case, mathematical modeling of HIV can help in estimating the lifespan of infected cells, evaluating therapeutic efficacy, and realizing that new virions require a host with deoxyribonucleic acid (DNA) to replicate. Currently, a variety of therapy options are available to help people with HIV get better. However, controlling the virus is the only option available and curing the disease still requires seamless efforts in the research domain.Mathematical models have been created to investigate the dynamic properties of cell populations using parameters including latent reservoirs, immunological responses, total carrying capacity, and time-delayed fractional differential operators based on literature reviews (see [1–5]). Immature infected cells also known as latent stage are those that have been infected but are not yet infectious, which is considered in the present study. Immune responses such as CD4+ T-cells/CD8+ T-cells are activated when a foreign agent enters the body, causing the body’s alarm system to go into overdrive [6, 7].When systems are in motion, there will always be a degree of lag time. Time delays are inevitable because it has an ability to cause a significant impact on cell populations. When modeling the kinetics of HIV-1 infection, two types of delays are taken into account. One situation is that the uninfected cells interact with infected/free virions and there is an intracellular temporal delay. Besides, after a foreign agent has been ingested, immune response cells must be activated, which results in a delay in the immunological response. This study aims to explore the effect of intracellular time delay with respect to infected cell population [8, 9]. Initially, ART is given to every primarily infected HIV person, but based on the stage of infection the level of drugs may be redefined. However, it is a challenging task with respect to the immune boosters which may vary in the individual. Besides, it can be seen that, if there is a change in the period of drugs provided, then it will reflect in the virion populations. Hence, this study explores the effect of antiretroviral therapies which are suggested as posttreatment for a long period of time. The proposed model’s stability analysis will provide some insights into disease progression in relation to the system parameters such as infection rate and time delays. Followed by, bifurcation analysis is used to determine the threshold value of the significant parameter that has the potential to cause fluctuations in the cell population. The optimal control algorithm is designed by selecting cell populations, which aids in better understanding and extraction of system parameters, resulting in cell populations with stable staining [10–13]. Recently, the models were proposed on Zika virus, HIV, SARS-CoV-2, and other viral dynamics models such as maize streak virus in maize spread by leaf hopper mostly in Africa, canine distemper virus, and rabies epidemics in red fox with respect to significant factors such as vaccination parameter involved models and optimal control strategies (for more details, refer [14–22]). Distinct from the existing models, the present study focuses on considering the effect of time delays and two kinds of antiretroviral therapies. The overall contribution is listed in the following:(1)
This paper models the dynamics of HIV-1 infection by considering the factors such as healthy CD4+ cells, latent reservoirs, infected cells, free virions, and immune responses also considering the effect of time delay.(2)
Positivity and boundedness of solutions of the differential model are proved. The reproduction number is determined through the next-generation matrix, which helps to identify the community spread.(3)
Conditions for the existence of Hopf bifurcation are proved by choosing the intracellular time delay as a bifurcation parameter.(4)
Optimal control algorithm is designed to ensure the stabilization of the proposed model.(5)
Numerical simulations are performed to validate the proposed theoretical frameworks.


## 2. Model of Posttreatment HIV-1 Viral Dynamics

The schematic representation of our model is given below.(1)x˙t=s−γxt−1−ε1tβxtyt−τ,l˙t=αL1−ε1tβxtyt+ρ−a−dLlt,y˙t=1−αL1−ε1tβxtyt+alt−δyt−mytzt,u˙t=1−ε2tpyt−cut,z˙t=bzytKB+ytzt−dztKD+ytzt−μzt.In model (1), the first equation represents the rate of change in the susceptible cell populations, where s denotes the rate of production, γ is decay rate and ɛ1,ɛ2 denotes the efficacy of antiretroviral therapies. Consider that εit, where i=1,2, in 0,1; if εi=1, where i=1,2., then cent percent the treatment is effective which makes the infection zero, and if εit=0, i=1,2, then, no progress in the therapy. β represents the rate of infection between the infected and uninfected cell population. τ stands for intracellular time delay. The second equation describes the state of latent infection; that is, target cells are infected but not yet infectious. Suppose if the infected cells are matured enough to infect the susceptible cells. ρ, αL, and dL are scalars. The third equation explains the rate of change in the infectious cell population with the death rate δ, migrated from latent infection a, removing the infection by immune responses at the rate of m. The fourth equation explores the rate of change in the free virions; the rate of proliferation from the infection is given by p and the decay rate is c. Finally, the production of effector cells can have the maximum proliferation in an infected cell with a maximum rate bz and is given by the term bzyt/yt+KBzt. The loss of functional effector cells is defined by dzyt/yt+KDzt, with an assumption KD>KB.


### 2.1. Positivity and Boundedness

In order to ensure the convergence of the solutions of the model, it becomes necessary to ensure that solutions of the state variables in model (1) are positive and ultimately bounded.Theorem 1.
Assume thatxt, lt, yt, ut, zt) be a solution of the proposed model (1) with the initial conditions xt>0, lt>0, yt>0, ut>0, zt>0. It is positive and ultimately bounded for t>0.Proof.
Consider the proof regardingxt>0 for all t>0. In this regard, assume that there exists t1>0, which implies xt1=0, xt>0, t∈0,t1 such that x˙t≤0. From (1), it is clear that x˙t1=s, s>0, which is a contradiction and leads to the proof that xt>0,∀t>0.
Similarly, the proof can be extended for all the remaining state equations.(2)lt=l0e−a+dL−ρt+∫0tαL1−ϵ1tβxξyξea+dL−ρξ−tdξ,zt=z0e∫0tbzyξ/KB+yξzξ−dzyξ/KD+yξzξ−μdξ,yt=y0e−δt+∫0t1−αL1−ϵ1tβxξyξ+alξ−myξzξe−δξ−tdξ,ut=u0e−ct+∫0t1−ϵ2tyecξ−tdξ.
Now, to prove {lt, yt, ut, zt>0, ∀t>0. And consider that t2>0 and define(3)minlt2,yt2,ut2,zt2=0.
Iflt2=0, lt>0 for t∈0,t2 and yt>0, ut>0, zt>0 for t∈0,t2,l˙t2 = 0, then we have l˙t2=αL1−ϵ1tβxt2yt2>0, which is a contradiction.
Ifyt2=0, yt>0 for t∈0,t2 and {lt, ut, zt>0 for t∈0,t2 with y˙t2≤0. However, from (1), one can have y˙t2=alt2>0.
Similarly, forut2=0 and zt2=0 is also a contradiction. Thus, xt,lt,yt,ut,zt>0, ∀t>0. To proceed with ensuring the boundedness of the solutions, we extend the results of the positivity of the solution for model (1).(4)x˙t≤s−γxt.
Taking the limits will lead to(5)limt⟶∞supxt≤sγ.
LetHt=xt+lt+yt(6)H˙t=x˙t+l˙t+y˙tH˙t=s−γxt−1−ε1tβxtyt+αL1−ε1tβxtyt+ρ−a−dLlt+1−αL1−ε1tβxtyt+alt−δyt−mytzt≤s−ϒxt+ρ−dLlt−δyt≤s+σHt.
Define thatσ = minγ,ρ−dL,δ and(7)limt⟶∞supHt≤sσ,
To prove the ultimate boundedness of free virions, the same approach is followed:(8)u˙t=pyt−cut=psσ1−ϵ2t−cutlimt→∞suput≤ps1−ϵ2tσc.
Sinceut cannot be negative, if ϵ2t=0, then suput=0. Similarly, for zt, we get(9)z˙t=bzytKB+ytzt−dztKD+ytzt−μzt,z˙tz=bzytKB+yt−dzytKD+yt−μ.
Applying integration on both sides(10)limt⟶∞supzt≤bzΩ1−dzΩ2−μ,where(11)Ω1=bzs/σKB+s/σ,Ω2=dzs/σKD+s/σ.
Hence, it is proved that all the state variables solutions are ultimately bounded.
Consider(12)Δ=x,l,y,u,z∈C+5:xt≤sσ,lt≤sσ,yt≤sσ,ut≤ps1−ε2tσc,zt≤bzΩ1−dzΩ2−μ,where σ, Ω1, and Ω2 are given in the above equation. From the given theorem, it is clear that, within the region, Δ is a positive invariant.


## 2.1. Positivity and Boundedness

In order to ensure the convergence of the solutions of the model, it becomes necessary to ensure that solutions of the state variables in model (1) are positive and ultimately bounded.Theorem 1.
Assume thatxt, lt, yt, ut, zt) be a solution of the proposed model (1) with the initial conditions xt>0, lt>0, yt>0, ut>0, zt>0. It is positive and ultimately bounded for t>0.Proof.
Consider the proof regardingxt>0 for all t>0. In this regard, assume that there exists t1>0, which implies xt1=0, xt>0, t∈0,t1 such that x˙t≤0. From (1), it is clear that x˙t1=s, s>0, which is a contradiction and leads to the proof that xt>0,∀t>0.
Similarly, the proof can be extended for all the remaining state equations.(2)lt=l0e−a+dL−ρt+∫0tαL1−ϵ1tβxξyξea+dL−ρξ−tdξ,zt=z0e∫0tbzyξ/KB+yξzξ−dzyξ/KD+yξzξ−μdξ,yt=y0e−δt+∫0t1−αL1−ϵ1tβxξyξ+alξ−myξzξe−δξ−tdξ,ut=u0e−ct+∫0t1−ϵ2tyecξ−tdξ.
Now, to prove {lt, yt, ut, zt>0, ∀t>0. And consider that t2>0 and define(3)minlt2,yt2,ut2,zt2=0.
Iflt2=0, lt>0 for t∈0,t2 and yt>0, ut>0, zt>0 for t∈0,t2,l˙t2 = 0, then we have l˙t2=αL1−ϵ1tβxt2yt2>0, which is a contradiction.
Ifyt2=0, yt>0 for t∈0,t2 and {lt, ut, zt>0 for t∈0,t2 with y˙t2≤0. However, from (1), one can have y˙t2=alt2>0.
Similarly, forut2=0 and zt2=0 is also a contradiction. Thus, xt,lt,yt,ut,zt>0, ∀t>0. To proceed with ensuring the boundedness of the solutions, we extend the results of the positivity of the solution for model (1).(4)x˙t≤s−γxt.
Taking the limits will lead to(5)limt⟶∞supxt≤sγ.
LetHt=xt+lt+yt(6)H˙t=x˙t+l˙t+y˙tH˙t=s−γxt−1−ε1tβxtyt+αL1−ε1tβxtyt+ρ−a−dLlt+1−αL1−ε1tβxtyt+alt−δyt−mytzt≤s−ϒxt+ρ−dLlt−δyt≤s+σHt.
Define thatσ = minγ,ρ−dL,δ and(7)limt⟶∞supHt≤sσ,
To prove the ultimate boundedness of free virions, the same approach is followed:(8)u˙t=pyt−cut=psσ1−ϵ2t−cutlimt→∞suput≤ps1−ϵ2tσc.
Sinceut cannot be negative, if ϵ2t=0, then suput=0. Similarly, for zt, we get(9)z˙t=bzytKB+ytzt−dztKD+ytzt−μzt,z˙tz=bzytKB+yt−dzytKD+yt−μ.
Applying integration on both sides(10)limt⟶∞supzt≤bzΩ1−dzΩ2−μ,where(11)Ω1=bzs/σKB+s/σ,Ω2=dzs/σKD+s/σ.
Hence, it is proved that all the state variables solutions are ultimately bounded.
Consider(12)Δ=x,l,y,u,z∈C+5:xt≤sσ,lt≤sσ,yt≤sσ,ut≤ps1−ε2tσc,zt≤bzΩ1−dzΩ2−μ,where σ, Ω1, and Ω2 are given in the above equation. From the given theorem, it is clear that, within the region, Δ is a positive invariant.


## 3. Equilibria

This section describes the derivation of equilibria of PTC HIV model based on three cases such as disease-free, immune-free equilibrium, and endemic equilibrium.(i)
Model (1) without infection exhibits disease-free equilibrium E0=x0,0,0,0,0, where x0=s/d.(ii)
IfR0>1, model (1) has an immune-free equilibrium E1=x1,l1,y1,u1,0 with the coefficients(13)x1=δa+dL−ρβ1−ϵ1taαL+1−αLa+dL−ρ,l1=αLβ1−ϵ1tx1y1a+dL−ρ,u1=p1−ϵ2t1c,y1=dR0−1β1−ϵ1t.(iii)
The endemic equilibrium of the systemE2x+, l+, y+, u+, z+ is given by(14)x+=sγ+β1−ϵ1ty+,l+=αL1−ϵ1tβx+y+a+dL−ρy+=−Ky+Ky2−4IC2I,u+=p+1−ϵ2tc,z+=a+dL−ρ1−αL1−ϵ1tβx++aαL1−ϵ1tβx+y+−δa+dL−ρma+dL−ρ,with Ky=KDbz−dzKB−μKB−μKD and I=bz−dz−μ,C=−μKBKD.


### 3.1. Basic Reproduction Number

In epidemiology, a basic reproduction number is the average number of persons an affected person can transmit the secondary infection. The basic reproduction number is an indicator that helps to determine the community spread of infection, which can be calculated using the next generation matrix. For the proposed model, the basic reproduction number, say,R0, is determined for various situations and described in the following sections.


#### 3.1.1. For Disease-Free Equilibrium

The basic reproduction number of the model without a viral latent reservoir is(15)1−ϵ1tβxtyt−δyt>0,1−ϵ1tβxt−δyt>0.Sinceyt≠0, then(16)R0=1−ϵ1tβ⋅sγ⋅1δ.


#### 3.1.2. For Immune-Free Equilibrium

The basic reproduction number of the model with immune-free equilibrium is(17)ℱ=αL1−ε1tβxtyt1−αL1−ε1tβxtyt,V=ρ−a−dLltalt−δyt.Then,ℱ and V help us to find the next generation matrix FV−1 as calculated:(18)F=0αL1−ε1tβxt01−αL1−ε1tβxt,V=ρ−a−dL0a−δ,V−1=1−δρ−a−dL−δ0−aρ−a−dL,FV−1=0αL1−ε1tβxt01−αL1−ε1tβxt⋅1−δρ−a−dL−δ0−aρ−a−dL,=−aαL1−ε1tβxtδa+dL−ρ−αL1−ε1tβxtδ−a1−αL1−ε1tβxtδa+dL−ρ−1−ε1t1−αLβxtδ.Now, the characteristic equation for the above matrix is(19)λ2−−aαL1−ε1tβxtδa+dL−ρ−1−ε1t1−αLβxtδ,λ+−αL1−ε1tβxtδ−a1−αL1−ε1tβxtδa+dL−ρ=0,where x0=s/ϒ, then the basic reproduction number is given by(20)R0=1−ε1t⋅β⋅1−αL+aαLa+dL−ρ⋅sϒ⋅1δ.


#### 3.1.3. For Endemic Equilibrium

The basic reproduction number of the model with endemic equilibrium is given by(21)ℱ=αL1−ε1tβxtyt1−αL1−ε1tβxtyt,V=ρ−a−dLltalt−δyt.Then,ℱ and V help us to find the next generation matrix FV−1 as calculated below:(22)F=0αL1−ε1tβxt01−αL1−ε1tβxt,V=ρ−a−dL0a−δ,V−1=1−δρ−a−dL−δ0−aρ−a−dL,FV−1=0αL1−ε1tβxt01−αL1−ε1tβxt⋅1−δρ−a−dL−δ0−aρ−a−dL,=−aαL1−ε1tβxtδa+dL−ρ−αL1−ε1tβxtδ−a1−αL1−ε1tβxtδa+dL−ρ−1−ε1t1−αLβxtδ.Now, the characteristic equation for the above matrix is(23)λ2−−aαL1−ε1tβxtδa+dL−ρ−1−ε1t1−αLβxtδ,λ+−αL1−ε1tβxtδ−a1−αL1−ε1tβxtδa+dL−ρ=0.The basic reproduction number is calculated from the next generation matrix is given by(24)λ=2Isβ1−ε1taαL+a+dL−ρ1−αLδa+dL−ρ2Iγ+β1−ε1t−Ky+Ky2−4IC,where Ky, I, and C are described above. The spectral radius or the largest eigenvalue of the next generation matrix is called the basic reproduction number.


## 3.1. Basic Reproduction Number

In epidemiology, a basic reproduction number is the average number of persons an affected person can transmit the secondary infection. The basic reproduction number is an indicator that helps to determine the community spread of infection, which can be calculated using the next generation matrix. For the proposed model, the basic reproduction number, say,R0, is determined for various situations and described in the following sections.


### 3.1.1. For Disease-Free Equilibrium

The basic reproduction number of the model without a viral latent reservoir is(15)1−ϵ1tβxtyt−δyt>0,1−ϵ1tβxt−δyt>0.Sinceyt≠0, then(16)R0=1−ϵ1tβ⋅sγ⋅1δ.


### 3.1.2. For Immune-Free Equilibrium

The basic reproduction number of the model with immune-free equilibrium is(17)ℱ=αL1−ε1tβxtyt1−αL1−ε1tβxtyt,V=ρ−a−dLltalt−δyt.Then,ℱ and V help us to find the next generation matrix FV−1 as calculated:(18)F=0αL1−ε1tβxt01−αL1−ε1tβxt,V=ρ−a−dL0a−δ,V−1=1−δρ−a−dL−δ0−aρ−a−dL,FV−1=0αL1−ε1tβxt01−αL1−ε1tβxt⋅1−δρ−a−dL−δ0−aρ−a−dL,=−aαL1−ε1tβxtδa+dL−ρ−αL1−ε1tβxtδ−a1−αL1−ε1tβxtδa+dL−ρ−1−ε1t1−αLβxtδ.Now, the characteristic equation for the above matrix is(19)λ2−−aαL1−ε1tβxtδa+dL−ρ−1−ε1t1−αLβxtδ,λ+−αL1−ε1tβxtδ−a1−αL1−ε1tβxtδa+dL−ρ=0,where x0=s/ϒ, then the basic reproduction number is given by(20)R0=1−ε1t⋅β⋅1−αL+aαLa+dL−ρ⋅sϒ⋅1δ.


### 3.1.3. For Endemic Equilibrium

The basic reproduction number of the model with endemic equilibrium is given by(21)ℱ=αL1−ε1tβxtyt1−αL1−ε1tβxtyt,V=ρ−a−dLltalt−δyt.Then,ℱ and V help us to find the next generation matrix FV−1 as calculated below:(22)F=0αL1−ε1tβxt01−αL1−ε1tβxt,V=ρ−a−dL0a−δ,V−1=1−δρ−a−dL−δ0−aρ−a−dL,FV−1=0αL1−ε1tβxt01−αL1−ε1tβxt⋅1−δρ−a−dL−δ0−aρ−a−dL,=−aαL1−ε1tβxtδa+dL−ρ−αL1−ε1tβxtδ−a1−αL1−ε1tβxtδa+dL−ρ−1−ε1t1−αLβxtδ.Now, the characteristic equation for the above matrix is(23)λ2−−aαL1−ε1tβxtδa+dL−ρ−1−ε1t1−αLβxtδ,λ+−αL1−ε1tβxtδ−a1−αL1−ε1tβxtδa+dL−ρ=0.The basic reproduction number is calculated from the next generation matrix is given by(24)λ=2Isβ1−ε1taαL+a+dL−ρ1−αLδa+dL−ρ2Iγ+β1−ε1t−Ky+Ky2−4IC,where Ky, I, and C are described above. The spectral radius or the largest eigenvalue of the next generation matrix is called the basic reproduction number.


## 3.1.1. For Disease-Free Equilibrium

The basic reproduction number of the model without a viral latent reservoir is(15)1−ϵ1tβxtyt−δyt>0,1−ϵ1tβxt−δyt>0.Sinceyt≠0, then(16)R0=1−ϵ1tβ⋅sγ⋅1δ.


## 3.1.2. For Immune-Free Equilibrium

The basic reproduction number of the model with immune-free equilibrium is(17)ℱ=αL1−ε1tβxtyt1−αL1−ε1tβxtyt,V=ρ−a−dLltalt−δyt.Then,ℱ and V help us to find the next generation matrix FV−1 as calculated:(18)F=0αL1−ε1tβxt01−αL1−ε1tβxt,V=ρ−a−dL0a−δ,V−1=1−δρ−a−dL−δ0−aρ−a−dL,FV−1=0αL1−ε1tβxt01−αL1−ε1tβxt⋅1−δρ−a−dL−δ0−aρ−a−dL,=−aαL1−ε1tβxtδa+dL−ρ−αL1−ε1tβxtδ−a1−αL1−ε1tβxtδa+dL−ρ−1−ε1t1−αLβxtδ.Now, the characteristic equation for the above matrix is(19)λ2−−aαL1−ε1tβxtδa+dL−ρ−1−ε1t1−αLβxtδ,λ+−αL1−ε1tβxtδ−a1−αL1−ε1tβxtδa+dL−ρ=0,where x0=s/ϒ, then the basic reproduction number is given by(20)R0=1−ε1t⋅β⋅1−αL+aαLa+dL−ρ⋅sϒ⋅1δ.


## 3.1.3. For Endemic Equilibrium

The basic reproduction number of the model with endemic equilibrium is given by(21)ℱ=αL1−ε1tβxtyt1−αL1−ε1tβxtyt,V=ρ−a−dLltalt−δyt.Then,ℱ and V help us to find the next generation matrix FV−1 as calculated below:(22)F=0αL1−ε1tβxt01−αL1−ε1tβxt,V=ρ−a−dL0a−δ,V−1=1−δρ−a−dL−δ0−aρ−a−dL,FV−1=0αL1−ε1tβxt01−αL1−ε1tβxt⋅1−δρ−a−dL−δ0−aρ−a−dL,=−aαL1−ε1tβxtδa+dL−ρ−αL1−ε1tβxtδ−a1−αL1−ε1tβxtδa+dL−ρ−1−ε1t1−αLβxtδ.Now, the characteristic equation for the above matrix is(23)λ2−−aαL1−ε1tβxtδa+dL−ρ−1−ε1t1−αLβxtδ,λ+−αL1−ε1tβxtδ−a1−αL1−ε1tβxtδa+dL−ρ=0.The basic reproduction number is calculated from the next generation matrix is given by(24)λ=2Isβ1−ε1taαL+a+dL−ρ1−αLδa+dL−ρ2Iγ+β1−ε1t−Ky+Ky2−4IC,where Ky, I, and C are described above. The spectral radius or the largest eigenvalue of the next generation matrix is called the basic reproduction number.


## 4. Stability Analysis

LetE2x+,l+,y+,u+,z+ be any arbitrary positive equilibrium of the system (1). Then, the Jacobian matrix was evaluated at positive equilibrium’s E2 in the view of biological aspects. This leads us to the following characteristic polynomial. The plus signs are ignored for the convenience of calculations.(25)J=−γ−1−ε1tβy0−1−ε1tβxe−λτ00αL1−ε1tβyρ−a−dLαL1−ε1tβx00βy1−αL1−ε1taa330−my00p1−ε2t−c000a530a55,where a33=βx1−αL1−ε1t−mz−δ,a53=bzz/KB+y−dzz/KD+y−bzyz/KB+y2+dzyz/KD+y2 and a55=bzy/KB+y−μ−dzy/KD+y


### 4.1. Stability Analysis of the Model without Time Delay

The characteristic polynomial without delay is given by(26)λ5+P1λ4+P2λ3+P3λ2+P4λ+P5=0.Routh–Hurwitz criterion: definenth root-based Routh–Hurwitz matrix as follows:(27)Hn=P1100⋯0P3P2P11⋯0P5P4P3P2⋯1⋮⋮⋮0000⋯Pn,where Pj=0, if j>n. When n=5, the matrix is simplified into(28)H5=P11000P3P2P110P5P4P3P2100P5P4P30000P5.Suppose all the roots of the characteristic polynomial with negative real part, then the determinant of Routh–Hurwitz matrices are positive and vice versa. That is,detHI>0,1=1,2,…,5. This can be employed to verify the proof of Theorem 1. The necessary and sufficient condition to exist for the negative real part for equation (6) is P1>0,P1P2−P3>0,P1P2P3−P1P4+P32−P1P5>0,P1P2P3P4−P2P5−P1P42+P4P5−P3P3P4−P2P5+P1P4P5+P52>0, and detH5>0.The expansion of the coefficientsP1,P2,P3,P4,P5 is given in Appendix.


### 4.2. Stability Analysis of the Model with Time Delay

The characteristic polynomial of the Jacobian is(29)λ5+Q1λ4+Q2λ3+Q3λ2+Q4λ+Q5+e−λτR1λ3+R2λ2+R3λ+R4=0.The coefficientsQi, i=1,2,3,4,5, and Ri, i=1,2,3,4, are given in Appendix. We rewrite the above equation as(30)Qλ+e−λτRλ=0.Suppose some of the eigenvalues are purely imaginary, that is,λ=iω, then the characteristic equation becomes(31)iω5+Q1iω4+Q2iω3+Q3iω2+Q4iω+Q5+e−iωτR1iω3+R2iω2+R3iω+R4=0,iω5+Q1ω4−iQ2ω3−Q3ω2+iQ4ω+Q5+e−iωτ−R1iω3−R2ω2+iR3ω+R4=0,iω5−Q2ω3+Q4ω+Q1ω4−Q3ω2+Q5+cosωτ−isinωτ−iR1ω3−R2ω2+iR3ω+R4=0,iω5−Q2ω3+Q4ω+Q1ω4−Q3ω2+Q5+cosωτ−isinωτi−R1ω3+R3ω−R2ω2+R4=0.Equating real and imaginary parts in the above equation:(32)ω5−Q2ω3+Q4ω+cosωτ−R1ω3−R3ω−isinωτ−R2ω2+R4=0,(33)Q1ω4−Q3ω2+Q5+cosωτ−R2ω2+R4+isinωτ−R1ω3+R3ω=0.Squaring and adding equations (32) and (33), we get(34)ω5−Q2ω3+Q4ω2+Q1ω4−Q3ω2+Q52=cos2ωτ−R2ω2+R42+sin2ωτ−R1ω3+R3ω2+cos2ωτ−R1ω3+R3ω2+sin2ωτ−R2ω2+R42,ω5−Q2ω3+Q4ω2+Q1ω4−Q3ω2+Q52=−R1ω3+R3ω2+−R2ω2+R42.Simplifying the above equations leads to the following:−ω10+−Q12+2Q2ω8+−Q22+R12−2Q4+2Q1Q3ω6+−Q32+R22−2Q1Q5+2Q2Q4−2R1R3ω4+−Q42+R32+2Q3Q5−2R2R4ω2−Q52+R42=0.Since this is the differential equation of order 10, we get almost “10” roots such asω1=z1,ω2=z2,ω3=z3,ω4=z4,ω5=z5,ω6=z6,ω7=z7,ω8=z8,ω9=z9,andω10=z10.Eliminatingcosωτ from equations (32) and (33), we get(35)sinωτ=ω5−Q2ω3+Q4ωA−Q1ω4−Q3ω2+Q5BA2+B2,ωτ=sin−1ω5−Q2ω3+Q4ωA−Q1ω4−Q3ω2+Q5BA2+B2,τi=1ωsin−1ω5−Q2ω3+Q4ωA−Q1ω4−Q3ω2+Q5BA2+B2,where i=0,1,2,…, A=R4−R2ω2 and B=−R1ω3+R3ω.Hence, forτ=0, E2 is asymptotically stable by Routh–Hurwitz criterion [1], E2 remains stable for τ<τ0; we choose(36)τ0=minτj,which completes the proof.


## 4.1. Stability Analysis of the Model without Time Delay

The characteristic polynomial without delay is given by(26)λ5+P1λ4+P2λ3+P3λ2+P4λ+P5=0.Routh–Hurwitz criterion: definenth root-based Routh–Hurwitz matrix as follows:(27)Hn=P1100⋯0P3P2P11⋯0P5P4P3P2⋯1⋮⋮⋮0000⋯Pn,where Pj=0, if j>n. When n=5, the matrix is simplified into(28)H5=P11000P3P2P110P5P4P3P2100P5P4P30000P5.Suppose all the roots of the characteristic polynomial with negative real part, then the determinant of Routh–Hurwitz matrices are positive and vice versa. That is,detHI>0,1=1,2,…,5. This can be employed to verify the proof of Theorem 1. The necessary and sufficient condition to exist for the negative real part for equation (6) is P1>0,P1P2−P3>0,P1P2P3−P1P4+P32−P1P5>0,P1P2P3P4−P2P5−P1P42+P4P5−P3P3P4−P2P5+P1P4P5+P52>0, and detH5>0.The expansion of the coefficientsP1,P2,P3,P4,P5 is given in Appendix.


## 4.2. Stability Analysis of the Model with Time Delay

The characteristic polynomial of the Jacobian is(29)λ5+Q1λ4+Q2λ3+Q3λ2+Q4λ+Q5+e−λτR1λ3+R2λ2+R3λ+R4=0.The coefficientsQi, i=1,2,3,4,5, and Ri, i=1,2,3,4, are given in Appendix. We rewrite the above equation as(30)Qλ+e−λτRλ=0.Suppose some of the eigenvalues are purely imaginary, that is,λ=iω, then the characteristic equation becomes(31)iω5+Q1iω4+Q2iω3+Q3iω2+Q4iω+Q5+e−iωτR1iω3+R2iω2+R3iω+R4=0,iω5+Q1ω4−iQ2ω3−Q3ω2+iQ4ω+Q5+e−iωτ−R1iω3−R2ω2+iR3ω+R4=0,iω5−Q2ω3+Q4ω+Q1ω4−Q3ω2+Q5+cosωτ−isinωτ−iR1ω3−R2ω2+iR3ω+R4=0,iω5−Q2ω3+Q4ω+Q1ω4−Q3ω2+Q5+cosωτ−isinωτi−R1ω3+R3ω−R2ω2+R4=0.Equating real and imaginary parts in the above equation:(32)ω5−Q2ω3+Q4ω+cosωτ−R1ω3−R3ω−isinωτ−R2ω2+R4=0,(33)Q1ω4−Q3ω2+Q5+cosωτ−R2ω2+R4+isinωτ−R1ω3+R3ω=0.Squaring and adding equations (32) and (33), we get(34)ω5−Q2ω3+Q4ω2+Q1ω4−Q3ω2+Q52=cos2ωτ−R2ω2+R42+sin2ωτ−R1ω3+R3ω2+cos2ωτ−R1ω3+R3ω2+sin2ωτ−R2ω2+R42,ω5−Q2ω3+Q4ω2+Q1ω4−Q3ω2+Q52=−R1ω3+R3ω2+−R2ω2+R42.Simplifying the above equations leads to the following:−ω10+−Q12+2Q2ω8+−Q22+R12−2Q4+2Q1Q3ω6+−Q32+R22−2Q1Q5+2Q2Q4−2R1R3ω4+−Q42+R32+2Q3Q5−2R2R4ω2−Q52+R42=0.Since this is the differential equation of order 10, we get almost “10” roots such asω1=z1,ω2=z2,ω3=z3,ω4=z4,ω5=z5,ω6=z6,ω7=z7,ω8=z8,ω9=z9,andω10=z10.Eliminatingcosωτ from equations (32) and (33), we get(35)sinωτ=ω5−Q2ω3+Q4ωA−Q1ω4−Q3ω2+Q5BA2+B2,ωτ=sin−1ω5−Q2ω3+Q4ωA−Q1ω4−Q3ω2+Q5BA2+B2,τi=1ωsin−1ω5−Q2ω3+Q4ωA−Q1ω4−Q3ω2+Q5BA2+B2,where i=0,1,2,…, A=R4−R2ω2 and B=−R1ω3+R3ω.Hence, forτ=0, E2 is asymptotically stable by Routh–Hurwitz criterion [1], E2 remains stable for τ<τ0; we choose(36)τ0=minτj,which completes the proof.


## 5. Hopf Bifurcation Analysis

In general, any physical system will reflect the changes in the qualitative behavior subject to states and parameters changes; for instance, changes in the rate of infection will reflect in the cell populations. Hence, determining the significant parameters that have an ability to affect the stability of the systems is considered bifurcation parameters. By choosing bifurcation parameters, various kinds of solution nature can be realized; among that Hopf-type bifurcation explores the point where solution trajectories cross the origin, say, from negative to positive, in which the system has purely imaginary eigenvalues. The present model possesses the Hopf-type bifurcation while the bifurcation parameter exceeds the threshold value. The process of deriving the stability conditions and proving the existence of Hopf bifurcation are given below.Theorem 2.
Consider that the characteristic polynomial is of the form(37)f0λ+f1λe−λτ=0,where f0 and f1 are continuously differentiable with respect to λ. One of the roots is λτ=ατ+iωτ, where λτ is continuously differentiable with respect to τ, and satisfies ατ0=0 and ωτ0=ω0 for a positive real number t0. Denote(38)ϕω=f0iω2−f1iω2,which results in(39)signdReλdτ|τ=τ0=sign12ωdϕdω|ω=ω0.Proof.
Consider equation (38), and by taking the derivative of f0iω2 with respect to ω, will lead to(40)ddωf0iω2=ddωRef0iω2+Imf0iω2=2Ref0iω.Ref0′iωi+2Imf0f0iω.Imf0′iωi=2Ref0iω¯f0′iωi=−2Imf0iω¯f0′iω.
Then,(41)12ω.dϕdω=ddωf0iω2−f1iω2=1ωImf1iω¯f1′iω−f0iω¯f0′iω=Imf1iω2f1′iωωf1iω−f0iω2f0′iωωf0iω.
Sincef0iω02=f1iω02, we have(42)12ω.dϕdω|ω=ω0=f0iω02Imf1′iω0ω0f1iω0−f0′iω0ω0f0iω0.
Now, we turn to the left side of (38), calculating the derivative of both sides of f0λ+f1e−λτ=0 with respect to τ, we obtain(43)f0′λdλdτ+f1′λdλdτe−λτ−λ+τdλdτf1λe−λτ=0.
Thus,(44)dλdτ−1=f0′λ+f1λe−λτ−τf1λe−λτλf1λe−λτ=f0′λeλτ+f1′λλf1λ−τλ.
Sincef0iω0+f1iω0e−iω0τ0=0, we have(45)Redλdτ|τ=τ0−1=Ref0′iω0eiω0τ0+f1′iω0ω0f1iω0=Ref0′iω0ω0f0iω0i+Re−f1′iω0ω0f1iω0i=Imf1′iω0ω0f1iω0−f0′iω0ω0f0iω0.
Therefore,(46)signdReλdτ|τ=τ0=signRedλdτ|τ=τ0=signRedλdτ|τ=τ0−1=sign12ω.dϕdω|ω=ω0.
The proof is completed.


## 6. Optimal Control Design

Considering (1), based on two variables such as RTI and PI as controls, namely, ε1 and ε2. Here, ε1 represents the drug reverse transcriptase, and it safeguards the healthy CD4+ from infection, so that the healthy immune cells are maintained in the right proportion. Also, ε2 is represented the drug protease inhibitors, which maintains the release of the free virions to burst which are active and fully infected. In general, the treatment is initiated with an antiretroviral drug or the combination of two or more drugs, which creates some side effects while using a regular basis. Since the decision of choosing the drug combination is complex and, in this regard, an optimal strategy is a useful tool to understand the situation and make decisions. Now, optimal control provides different options with respect to estimating the costs of the drugs used for the therapies and observing the drug’s effectiveness for the disease. The therapy should last long as the optimal control helps us to find the suitable drug combination for the disease. Formulating the optimization problem based on optimal pair, existence is discussed in the following sections.


### 6.1. The Optimization Problem

In order to state the optimization problem, we first considerε1 and ε2 vary with time.(47)x˙t=s−γxt−1−ε1tβxtyt,l˙t=αL1−ε1tβxtyt+ρ−a−dLlt,y˙t=1−αL1−ε1tβxtyt+alt−δyt−mytzt,u˙t=p1−ε2tyt−cut,z˙t=bzytKB+ytzt−dzytKD+ytzt−μzt.The optimization problem is designed in terms of maximizing the following objective function constructed from the model parameters:(48)Jε1,ε2=∫0tfxt+zt+ut−A12ε12t+A22ε22tdt.Here, the upper boundtf denotes the period of the treatment and assumptions A1>0 and A2>0, respectively, stand for benefit and treatment costs. The scalars ε1t and ε2t are bounded and Lebesgue integrable. The objective of the control is to increase the uninfected cell populations through immune cells and decrease the cell count of free virions and infected cells. Hence, ε1∗,ε2∗ is the control pair that needs to be investigated. We assume that the control pair is nonempty, convex, and closed and it is integrable in the objective functional.(49)Jε1∗,ε2∗=maxJε1,ε2:ε1,ε2∈U,where U is the control set defined by U=ε1t,ε2t:εit measurable, 0≤εi≤1,t∈0,tf,i=1,2.


### 6.2. Optimal System

In order to investigate the properties of the optimal system, the Pontryagin’s minimum principle given in [13] is utilized and it provides necessary stability conditions for the designed optimal control problem. The advantage of the principle is that it handles (47)–(49) in terms of maximizing a Hamiltonian H through ɛ1 and ɛ2(50)Ht,x,l,y,u,z,yτ,ε1,ε2,λ=A12ε12+A22ε22−x−u−z+∑i=05λifi,with(51)f1=s−γx−1−ε1βxyτf2=αL1−ε1βxy+ρ−a−dLl,f3=1−αL1−ε1βxy+al−δy−myz,f4=p1−ε2y−cu,f5=bzyKB+yz−dzyKD+yz−μz.Based on the above discussion, the following theorem can be derived.Theorem 3.
For any optimal controlε1∗, ε2∗, and solutions x∗,l∗,y∗,u∗,z∗ of the corresponding state system (1), there exists adjoint variables λ1, λ2, λ3, λ4, and λ5 satisfying the equations(52)λ1′t=1+λ1tϒ+β1−ε1t+τyt+λ2tαLε1−1βyt+λ3t1−αLε1−1βytλ2′t=λ2t−ρ+a+dL−λ3taλ3′t=λ1tβxt1−ε1t+τ+λ2tαlε1t−1βxt+λ3t1−αLε1t−1βxt+δ+mzt+λ4tρε2t−1+λ5t−bzztkB+yt+dzztkD+yt+bzytztkB+yt2−dzytztkD+yt2λ4′t=1+cλ4tλ5′t=1+λ3tmyt+λ5tdzytKD+yt−bzytKB+yt+μ,with the transversality conditions(53)λitf=0,i=1,…,5.
Moreover, the optimal control is given by(54)ε1∗=min1,max0,−βA1λ1x∗tyτ∗,ε2∗=min1,max0,pA2λ4y∗t.Proof.
The adjoint equations and transversality conditions can be obtained by using Pontryagin’s minimum principle with delay, such that suppose if consider the equation in the vector forX=λ1,λ2,λ3,λ4,λ5.(55)λ1′t=−∂x˙t∂xt∂l˙t∂xt∂y˙t∂xt∂u˙t∂xt∂z˙t∂xtXT−∂H∂x,λ2′t=−∂x˙t∂lt∂l˙t∂lt∂y˙t∂lt∂u˙t∂lt∂z˙t∂ltXT−∂H∂l,λ3′t=−∂x˙t∂yt∂l˙t∂yt∂y˙t∂yt∂u˙t∂yt∂z˙t∂ytXT−∂H∂y,λ4′t=−∂x˙t∂ut∂l˙t∂ut∂y˙t∂ut∂u˙t∂ut∂z˙t∂utXT−∂H∂u,λ5′t=−∂x˙t∂zt∂l˙t∂zt∂y˙t∂zt∂u˙t∂zt∂z˙t∂ztXT−∂H∂z.
The optimal controlε1∗ and ε2∗ can be solved from the optimality conditions(56)∂H∂ε1t=0,∂H∂ε2t=0.
That is,(57)∂H∂ε1t=A1ε1+βxtyτλ1=0,∂H∂ε2t=A2ε2−λ4pyt=0.
Consider the optimal controlsε1∗ and ε2∗ defined in (54) are bounded. If we substitute ε1∗ and ε2∗ in systems (47) and (52), we obtain the following optimality system:(58)x˙t=s−ϒx∗t−1−ε1∗tβx∗ty∗t−Γ,l˙∗t=αL1−ε1∗tβx∗ty∗t+ρ−a−dLl∗t,y˙∗t=1−αL1−ε1∗tβx∗ty∗t+al∗t−δy∗t−my∗tz∗t,u˙∗t=p1−ε2∗ty∗t−cu∗t,z˙∗t=bzy∗tKB+y∗tz∗t−dzy∗tKD+y∗tz∗t−μz∗t,λ1′t=1+λ1tϒ+β1−ε1t+τyt+λ2tαLε1−1βyt+λ3t1−αLε1−1βytλ2′t=λ2t−ρ+a+dL−λ3taλ3′t=λ1tβxt1−ε1t+τ+λ2tαlε1t−1βxt+λ3t1−αLε1t−1βxt+δ+mzt+λ4tρε2t−1+λ5t−bzztkB+yt+dzztkD+yt+bzytztkB+yt2−dzytztkD+yt2λ4′t=1+cλ4t,λ5′t=1+λ3tmyt+λ5tdzytKD+yt−bzytKB+yt+μ,ε1∗=min1,max0,−βA1λ1x∗tyτ∗,ε2∗=min1,max0,pA2λ4y∗t,λitf=0,i=1,…,5.
The significance of stability conditions, validating the existence of Hopf bifurcation properties, and effect of the optimal controller are explored in the following section.


## 6.1. The Optimization Problem

In order to state the optimization problem, we first considerε1 and ε2 vary with time.(47)x˙t=s−γxt−1−ε1tβxtyt,l˙t=αL1−ε1tβxtyt+ρ−a−dLlt,y˙t=1−αL1−ε1tβxtyt+alt−δyt−mytzt,u˙t=p1−ε2tyt−cut,z˙t=bzytKB+ytzt−dzytKD+ytzt−μzt.The optimization problem is designed in terms of maximizing the following objective function constructed from the model parameters:(48)Jε1,ε2=∫0tfxt+zt+ut−A12ε12t+A22ε22tdt.Here, the upper boundtf denotes the period of the treatment and assumptions A1>0 and A2>0, respectively, stand for benefit and treatment costs. The scalars ε1t and ε2t are bounded and Lebesgue integrable. The objective of the control is to increase the uninfected cell populations through immune cells and decrease the cell count of free virions and infected cells. Hence, ε1∗,ε2∗ is the control pair that needs to be investigated. We assume that the control pair is nonempty, convex, and closed and it is integrable in the objective functional.(49)Jε1∗,ε2∗=maxJε1,ε2:ε1,ε2∈U,where U is the control set defined by U=ε1t,ε2t:εit measurable, 0≤εi≤1,t∈0,tf,i=1,2.


## 6.2. Optimal System

In order to investigate the properties of the optimal system, the Pontryagin’s minimum principle given in [13] is utilized and it provides necessary stability conditions for the designed optimal control problem. The advantage of the principle is that it handles (47)–(49) in terms of maximizing a Hamiltonian H through ɛ1 and ɛ2(50)Ht,x,l,y,u,z,yτ,ε1,ε2,λ=A12ε12+A22ε22−x−u−z+∑i=05λifi,with(51)f1=s−γx−1−ε1βxyτf2=αL1−ε1βxy+ρ−a−dLl,f3=1−αL1−ε1βxy+al−δy−myz,f4=p1−ε2y−cu,f5=bzyKB+yz−dzyKD+yz−μz.Based on the above discussion, the following theorem can be derived.Theorem 3.
For any optimal controlε1∗, ε2∗, and solutions x∗,l∗,y∗,u∗,z∗ of the corresponding state system (1), there exists adjoint variables λ1, λ2, λ3, λ4, and λ5 satisfying the equations(52)λ1′t=1+λ1tϒ+β1−ε1t+τyt+λ2tαLε1−1βyt+λ3t1−αLε1−1βytλ2′t=λ2t−ρ+a+dL−λ3taλ3′t=λ1tβxt1−ε1t+τ+λ2tαlε1t−1βxt+λ3t1−αLε1t−1βxt+δ+mzt+λ4tρε2t−1+λ5t−bzztkB+yt+dzztkD+yt+bzytztkB+yt2−dzytztkD+yt2λ4′t=1+cλ4tλ5′t=1+λ3tmyt+λ5tdzytKD+yt−bzytKB+yt+μ,with the transversality conditions(53)λitf=0,i=1,…,5.
Moreover, the optimal control is given by(54)ε1∗=min1,max0,−βA1λ1x∗tyτ∗,ε2∗=min1,max0,pA2λ4y∗t.Proof.
The adjoint equations and transversality conditions can be obtained by using Pontryagin’s minimum principle with delay, such that suppose if consider the equation in the vector forX=λ1,λ2,λ3,λ4,λ5.(55)λ1′t=−∂x˙t∂xt∂l˙t∂xt∂y˙t∂xt∂u˙t∂xt∂z˙t∂xtXT−∂H∂x,λ2′t=−∂x˙t∂lt∂l˙t∂lt∂y˙t∂lt∂u˙t∂lt∂z˙t∂ltXT−∂H∂l,λ3′t=−∂x˙t∂yt∂l˙t∂yt∂y˙t∂yt∂u˙t∂yt∂z˙t∂ytXT−∂H∂y,λ4′t=−∂x˙t∂ut∂l˙t∂ut∂y˙t∂ut∂u˙t∂ut∂z˙t∂utXT−∂H∂u,λ5′t=−∂x˙t∂zt∂l˙t∂zt∂y˙t∂zt∂u˙t∂zt∂z˙t∂ztXT−∂H∂z.
The optimal controlε1∗ and ε2∗ can be solved from the optimality conditions(56)∂H∂ε1t=0,∂H∂ε2t=0.
That is,(57)∂H∂ε1t=A1ε1+βxtyτλ1=0,∂H∂ε2t=A2ε2−λ4pyt=0.
Consider the optimal controlsε1∗ and ε2∗ defined in (54) are bounded. If we substitute ε1∗ and ε2∗ in systems (47) and (52), we obtain the following optimality system:(58)x˙t=s−ϒx∗t−1−ε1∗tβx∗ty∗t−Γ,l˙∗t=αL1−ε1∗tβx∗ty∗t+ρ−a−dLl∗t,y˙∗t=1−αL1−ε1∗tβx∗ty∗t+al∗t−δy∗t−my∗tz∗t,u˙∗t=p1−ε2∗ty∗t−cu∗t,z˙∗t=bzy∗tKB+y∗tz∗t−dzy∗tKD+y∗tz∗t−μz∗t,λ1′t=1+λ1tϒ+β1−ε1t+τyt+λ2tαLε1−1βyt+λ3t1−αLε1−1βytλ2′t=λ2t−ρ+a+dL−λ3taλ3′t=λ1tβxt1−ε1t+τ+λ2tαlε1t−1βxt+λ3t1−αLε1t−1βxt+δ+mzt+λ4tρε2t−1+λ5t−bzztkB+yt+dzztkD+yt+bzytztkB+yt2−dzytztkD+yt2λ4′t=1+cλ4t,λ5′t=1+λ3tmyt+λ5tdzytKD+yt−bzytKB+yt+μ,ε1∗=min1,max0,−βA1λ1x∗tyτ∗,ε2∗=min1,max0,pA2λ4y∗t,λitf=0,i=1,…,5.
The significance of stability conditions, validating the existence of Hopf bifurcation properties, and effect of the optimal controller are explored in the following section.


## 7. Numerical Simulation

This section describes the numerical evaluations of the proposed model (1) by choosing the experimental range of parameter values provided in Table 1. The simulations are performed through Runge–Kutta fourth-order numerical approximation scheme. The outcomes of the cell simulations corresponding to cell populations are demonstrated through phase-space diagrams. Figures 1 and 2 illustrate the solution behavior of cell populations within the threshold and exceeding the threshold of time delay, respectively. Also, Figure 3 explores the effect of optimal control design, in which the objective is achieved by increasing the number of uninfected cell populations and diminishing the number of virions. In addition, Figures 4–7 provide insights by comparing the controlled and uncontrolled system states such as uninfected, infected, latent, and free virions, respectively. Figures 2 shows the interpretation of variables of the equilibrium when τ<0.93. The top left corner in the panel of Figure 2 evidence the rate of change in healthy CD4+ and infected cells and the bottom left depicts the uninfected cells and the free virions. Similarly, the phase portrait of all the cell populations such as uninfected and infected cells with free virions and immune cells is accordingly picturized. The drug parameter ε1t and the transmission rate β have a high influence on the stability of the equilibrium. That is, for 0≤ε1t≤0.6, the equilibrium is stable. For 1.5×10−7≤β≤1.5×10−6, the equilibrium remains stable. When the transmission rate gets higher or lesser, it loses the stability of equilibrium, where αL, r, μ, and a also affect the stability when it increases rapidly. The rest of the parameters are having less significance in terms of stability when compared with the remaining parameters.Table 1
Parameter values and their sources.
NotationsValuesReferencess105 cells/mL/day[23]γ0.01—β1.5×10−6mL/day[24, 25]δ1day−1[26]N2000[27, 28]p2000day−1—C23day−1[29]ϵ0.001—A0.001day−1—dL0.004day−1[30]t1/244 months[31]ρ0.0045day−1[32]αL10−6—sz1cells/mL/day—bz1day−1[33]KB0.1cells/mL[34]dz2day−1[35]KD5cells/mL[34]μ0.0002day−1—M0.42mL/cells/day—Figure 1
Unstable responses for the time delayτ>0.93.Figure 2
Stable responses for the time delayτ<0.93.Figure 3
Optimal control-based state responses with respect to time.Figure 4
Comparison of controlled and uncontrolled susceptible cell populations.Figure 5
Comparison of controlled and uncontrolled infected cell populations.Figure 6
Comparison of controlled and uncontrolled latent cell populations.Figure 7
Comparison of controlled and uncontrolled free virions.Figure1 represents the state trajectories where the equilibrium point loses its stability due to time delay that makes the transmission rate slow. If the drug is already administered into the body, the drug starts to function among the cells; there is a period for the drug to make progress and with respect to time delays, it becomes necessary to have the effect of time response, and if the drug efficacy is not up to expectation, then the endemic equilibrium will loose the stability.Figure8 shows that the uninfected cells start to bifurcate as τ crosses the value 0.93.Figure 8
Bifurcation behavior for time delay with respect to uninfected cell populations.Figures3 represents the uninfected cell population after the optimal control is established. The top figures show the increasing cell population of the uninfected cells and the reducing rate of the infected cells; it is visible after a certain stage the immune cells maintained in proportion, with no increasing number of cells. Similarly, the free virions are depleted in the process.Figures4–7 provide the comparisons between controlled and uncontrolled susceptible cell populations, infected cell populations, latent cell populations, and free virions. The graphs exhibit the controlled cell population after the implementation of optimal control.


## 8. Conclusion

The aim of the paper is to explore the dynamical analysis of posttreatment HIV-1 infection with respect to various significant parameters such as the effect of time delays, two different kinds of antiretroviral therapies, and loss of functional effector cells. By employing the Routh–Hurwitz criterion, the stability properties of the model with respect to discrete-type constant time-delay which is chosen as a bifurcation parameter have been presented. To ensure the effect of time delays, the existence of Hopf-type bifurcation in the behavior of solutions has been proved through proving the corresponding transversality conditions. To reveal the effect of combination of drug therapies, along with those parameters, the situation is modeled as an optimal control problem in terms of the objective function. Through objective function, the results explore the maximization in the number of uninfected cell population and minimization in the number of infected cells and corresponding results are picturized in the numerical section. In the future direction, the model can be extended in terms of fractional order differential operator, considering continuous-type time-delays, stochastic disturbances, and impulse in the antiviral therapy, which also plays a significant role in the system dynamics.



---

*Source: 9752628-2022-10-08.xml*