# Bifurcation and Stabillity Analysis of HIV Transmission Model with Optimal Control

Authors: Kumama Regassa Cheneke; Koya Purnachandra Rao; Geremew Kenassa Edessa
Journal: Journal of Mathematics
 (2021)

---

## Abstract

A mathematical model of HIV transmission is built and studied in this paper. The system’s equilibrium is calculated. A next-generation matrix is used to calculate the reproduction number. The novel method is used to examine the developed model’s bifurcation and equilibrium stability. The stability analysis result shows that the disease-free equilibrium is locally asymptotically stable if0<R0<1 but unstable if R0>1. However, the endemic equilibrium is locally and globally asymptotically stable if R0>1 and unstable otherwise. The sensitivity analysis shows that the most sensitive parameter that contributes to increasing of the reproduction number is the transmission rate β2 of HIV transmission from HIV individuals to susceptible individuals and the parameter that contributes to the decreasing of the reproduction number is identified as progression rate η of HIV-infected individuals to AIDS individuals. Furthermore, it is observed that as we change η from 0.1 to 1, the reproduction number value decreases from 1.205 to 1.189, where the constant value of β2=0.1. On the other hand, as we change the value of β2 from 0.1 to 1, the value of the reproduction number increases from 0.205 to 1.347, where the constant value of η=0.1. Further, the developed model is extended to the optimal control model of HIV/AIDS transmission, and the cost-effectiveness of the control strategy is analyzed. Pontraygin’s Maximum Principle (PMP) is applied in the construction of the Hamiltonian function. Moreover, the optimal system is solved using forward and backward Runge–Kutta fourth-order methods. The numerical simulation depicts the number of newly infected HIV individuals and the number of individuals at the AIDS stage reduced as a result of taking control measures. The cost-effectiveness study demonstrates that when combined and used, the preventative and treatment control measures are effective. MATLAB is used to run numerical simulations.


---

## Body

## 1. Introduction

Human immunodeficiency virus (HIV) is a retrovirus that attacks the human immune system and causes a highly killing disease called acquired immunodeficiency syndrome (AIDS) [1–6]. HIV was discovered in the early 1980, and it has been persisting in the population [2]. HIV is transmitted through unsafe sex, blood transfusion, breast feeding, materials exposed to the virus, and mother to child during pregnancy [3]. Currently, there is no curing treatment for HIV-infected individuals [4]. In 2018, the number of human individuals living with HIV is estimated to be 37.9 million and the number of dead individuals with AIDS-related disease is 1.2 million. Among HIV-infected individuals, about 62% are tested and taking antiretro therapy (ART) [4]. The data indicate that Africa is the continent that is highly exposed to the human immunodeficiency virus in the world. Particularly in 2018, Ethiopia has about 690,000 peoples living with HIV, 23,000 new people are infected with HIV, and 11,000 individuals dead with AIDS-related disease [5]. To control the transmission of the human immunodeficiency virus, different protective and treatment strategies are used [6]. Some of the strategies used in controlling the transmission and progression of HIV are using condom, be faithful, abstaining, and ART [7]. Even though, different control strategies are used to eradicate and combat the transmission of HIV in human population, the viruses still become one of the global issues that need attention to save human population from this merciless killing disease.In this study, we extended the classical SIA model of HIV to the SWIA model of HIV with optimal control problem to identify the best control measures that reduces the transmission and progression of the human immunodeficiency virus along with minimum cost. We have incorporated two control measures in the model and analyzed the cost-effectiveness of the controls.Mathematical models are efficient tools to describe and predict the transmission dynamics of disease in the population [8, 9].


## 2. Mathematical Model Formulation

In this study, we extended the SWIA model of HIV into optimal control problem which is significantly different from our previous work in approach and objective of the study. Similar to our previous work, the SWIA model of HIV is described as follows. (i) Susceptible individualsS: they are individuals who are free of the human immunodeficiency virus but have a chance of possible exposure to the virus having unsafe sexual practices. (ii) Window stage individuals I: they are new individuals infected with the human immunodeficiency virus but the presence of virus in the blood cannot be verified by the laboratory test. (iii) HIV stage individuals I: they are pre-AIDS human individuals whose blood test gives the positive results. (iv) AIDS stage individuals A: they are humans at advanced stage of infection caused by the human immunodeficiency virus and resistant to treatment as a result of weakening the human immune system.Moreover, the following assumptions are considered in the development of the model:(i)
Human population size is assumed to be nonconstant.(ii)
Individuals are recruited into susceptible population at the recruitment rateλ.(iii)
Human immunodeficiency virus transmitted at the constant transmission rate ofβ1,β2,andβ3 from W,I,andA, respectively.(iv)
The total population size at timet is denoted by Nt and given by(1)Nt=St+Wt+It+At.(v)
Window stage individuals transfer to HIV stage at the progression rateρ.(vi)
HIV stage individuals transfer to AIDS stage at the rate ofη.(vii)
All humans die naturally at the constant rateμ.(viii)
AIDS stage individuals die at the constant rateδ.(2)dSdt=λ−SNβ1W+β2I+β3A−μS,(3)dWdt=SNβ1W+β2I+β3A−ρ+μW,(4)dIdt=ρW−η+μI,(5)dAdt=ηI−δ+μA.withS(0)>0, W(0)>0, I(0)>0, and A(0)>0.


## 3. Mathematical Analysis of the Model

### 3.1. Well-Posedness of the Model

#### 3.1.1. Invariant Region

Theorem 1.
The setΩ⊂ℜ+4 is the invariant region of boundedness if for all initial solutions S0,W0,I0,A0∈Ω and for all t≥0, then St,Wt,It,At∈Ω. That is,(6)Ω=St,Wt,It,At:St,Wt,It,At∈ℜ+4,∀t≥0.Proof.
Considering model (4) and adding all corresponding terms on the left of equality and terms on right of equality, we get(7)dNdt≤λ−μN+δA.
Solving the preceding inequality and considering the expression as timet gets larger, we obtain(8)Nt≤λμ.
The preceding inequality shows that all solution variables are bounded inℜ+4. Therefore, using equation (4), Ω⊂ℜ+4 is the invariant set of solutions such that(9)Ω=St,Wt,It,At:St,+Wt+It+At≤λμ,∀t≥0.


#### 3.1.2. Positivity of Solutions

Theorem 2.
All solutions of model (4) are nonnegative for all time t provided that initial conditions are nonnegative.Proof.
Taking the first equation of model (4), we have(10)dSdt=λ−Sβ1W+β2I+β3AN−μS.
Neglecting the termλ, the foregoing equation is reduced to the next inequality(11)dSS≥−β1W+β2I+β3AN+μdt.
Solving the foregoing inequality over time interval0,t, we get(12)St≥S0e−μt−∫0tβ1Wξ+β2Iξ+β3Aξ/Nξdξ.
Since the initial conditionS0 and the exponential expression e−μt−∫0tβ1Wξ+β2Iξ+β3Aξ/Nξdξ are nonnegative, the solution variable St is nonnegative for all time t. Similarly, other solution variables are nonnegative.


#### 3.1.3. Existence and Uniqueness of Solutions

Theorem 3.
The solutions of model (4) exist and are unique in ℜ+4.Proof.
Since the expression on the right-hand side of model (4) is bounded and continuously differentiable, the proof directly follows from Cauchy–Lipschitz theorem (see [10]).
Therefore, the formulated model (4) is mathematically well-posed and biologically acceptable.


### 3.2. Equilibriums of the Model

#### 3.2.1. Disease-Free Equilibrium (DFE)

The disease-free equilibriumE0 is a point in the system where there is no disease in the population. The computed disease-free equilibrium of the model one is given by(13)E0=S0,0,0,0,where S0=λ/μ.


#### 3.2.2. Endemic Equilibrium (EE)

The endemic equilibrium of the model is the point where the disease persists in the population. The computed endemic equilibrium of model (4) is given by(14)E1=S∗,W∗,I∗,A∗,where S∗=a+b+cλ/R0+μa+b+c−1,W∗=λ/ρ+μR0−1/R0+μa+b+c−1,I∗=ρW∗/η+μ,A∗=ηI∗/δ+μ,a=1/ρ+μ,b=ρ/ρ+μη+μ,andc=ηρ/ρ+μη+μδ+μ.


### 3.3. Basic Reproduction NumberR0

The basic reproduction number is the average number of infected individuals produced by one infectious individual in the susceptible population during entire period of infection [11]. We have computed the basic reproduction number using the next-generation matrix employed in [12]. Accordingly, from model (4), let f=Sβ1W+β2I+β3A/N00 represent the newly infected individuals in the infectious compartments and v=ρ+μW−ρW+η+μI−ηI+δ+μA be the remaining terms of infected compartments. So that, the Jacobian matrices obtained from f and v at the disease-free equilibrium are denoted by F and V, respectively. So that,(15)F=β1β2β3000000,V=ρ+μ00−ρη+μ00−ηδ+μ.The next-generation matrixFV−1 is computed as follows:(16)FV−1=β1ρ+μ+β2ρρ+μη+μ+β3ηρρ+μη+μδ+μβ2η+μ+β3ηη+μδ+μβ3δ+μ000000.The computed eigenvalues of the foregoing next-generation matrix are given by(17)λ1=β1ρ+μ+β2ρρ+μη+μ+β3ηρρ+μη+μδ+μ,λ2=0,λ3=0.The basic reproduction numberR0=ρFV−1 and given by(18)R0=β1ρ+μ+β2ρρ+μη+μ+β3ηρρ+μη+μδ+μ.


### 3.4. Stability Analysis of Disease-Free Equilibrium

Theorem 4.
The disease-free equilibrium becomes locally asymptotically stable ifR0<1 and becomes unstable if R0>1.Proof.
To determine the local stability of the disease-free equilibrium, we constructed a JacobianJ from model (2) at the disease-free equilibrium as follows:(19)J=−μ−β1−β2−β30β1−ρ+μβ2β30ρ−η+μ000η−δ+μ.
The eigenvalues of the above matrix can be obtained by solving the characteristic equation as follows:(20)detJ−λI=0.
The characteristic equation can be written as(21)λ+μa0λ3+a1λ2+a2λ+a3=0,where a0=1,a1=η+ρ+δ+3μ−β1,a2=ηδ+2ημ+2δμ+ηη+δ+2μ+3μ2−β1η+δ+2μ−β2,anda3=ημ2+δμ2+ρμ2+μ3+ηδμ+ηδρ+ημρ+δμρ−β1μ2+ηδ+ημ+δμ−β2δρ+μρ−β3ηρ.
By the Hurwitz–Routh principle, all eigenvalues of the characteristic equation are negative if and only if the following conditions are satisfied:(22)a0>0,a1a2>a3.
Therefore, by stability theory of differential equations, the disease-free equilibrium becomes locally asymptotically stable if0<R0<1 and becomes unstable if R0>1.Theorem 5.
The disease-free equilibrium of the constructed model (4) is globally asymptotically unstable if 0<R0<1 and unstable otherwise.Proof.
To prove global stability of the disease-free equilibrium, we apply the methods applied in [13] as follows.
Accordingly, letX∈ℝ1 be uninfected individuals in class S and Y∈ℝ3 be individuals in infected class I,U,A. Hence, we write model (4) as follows:(23)dXdt=HX,Y,dYdt=GX,Y,GX,0=0.
The disease-free equilibrium of the foregoing system is computed as(24)E01=X0,0.
To guarantee global stability, the method we applied must meet conditions H1 and H2 stated as follows:H1. For dX/dt=HX,0,X0 is globally asymptotically stableH2. GX,Y=PY−G^X,Y,G^X,Y≥0 for X,Y∈Ω
Here,P=DYGX,0 be the Metzler matrix and Ω is a region where solutions are acceptable.Theorem 6.
The disease-free equilibrium of model (4) is globally asymptotically stable in a region Ω if R0<1 and unstable whenever R0>1 provided that the stated two conditions H1 and H2 are satisfied.Proof.
Considering model (4), we have(25)HX,0=λ−μS=HS,0.
PuttingHX,0=0 and solving, we obtain S=λ/μ. Hence, X0=λ/μ,0.
We observe thatX0 is the globally asymptotically stable equilibrium of equation as follows:(26)dXdt=HX,0.
Further, considering the infected compartments of model (27), we have(27)dSdt=λ−SNβ1W+β2I+β3A−μS,dWdt=SNβ1W+β2I+β3A−ρ+μW,dIdt=ρW−η+μI,dAdt=ηI−δ+μA,GX,Y=β1W+β2I+β3A−ρ+μ00ρ−η+μ00η−δ+μWIA−β1W+β2I+β3A−SNβ1W+β2I+β3A00.
Further, at the disease-free equilibrium, the preceding equation is reduced to the form as follows:(28)GX,Y=−ρ+μ00ρ−η+μ00η−δ+μWIA−000.
Now, comparing the required condition and the preceding equation, we obtain(29)P=−ρ+μ00ρ−η+μ00η−δ+μ,Y=WIA,G^X,Y=000.
Now, we observe thatGX,Y=PY−G^X,Y, where G^X,Y≥0,∀X,Y.
Therefore, the disease-free equilibriumE0 is globally asymptotically stable if R0<1.


### 3.5. Bifurcation and Stability Analysis of Endemic Equilibrium

In this subsection, we simulate the endemic equilibrium versus reproduction number to identify the kind of bifurcation that occurs at the pointR0=1 and to determine the stability behavior of equilibriums if R0<1 and R0>1. The simulation of the endemic equilibrium S∗,W∗,I∗,A∗ versus R0 is given in Figure 1.Figure 1
Simulation of the endemic equilibrium versus reproduction number.Theorem 7.
The endemic equilibrium of model (2) undergoes forward bifurcation at R0=1.Proof.
(see [3, 14] and Figure 1).Theorem 8.
The endemic equilibriumE1 is globally asymptotically positively stable if R0>1.Proof.
(see [3, 14] and Figure 1).Theorem 9.
The endemic equilibrium of model (27) is locally asymptotically stable if R0>1 and unstable if 0<R0<1.Proof.
We have already computed the endemic equilibriumE1 as(30)E1=S∗,W∗,I∗,A∗,where S∗=a+b+cλ/R0+μa+b+c−1,W∗=λ/ρ+μR0−1/R0+μa+b+c−1,I∗=ρW∗/η+μ,A∗=ηI∗/δ+μ,a=1/ρ+μ,b=ρ/ρ+μη+μ,andc=ρη/ρ+μη+μδ+μ.
Clearly, the endemic equilibrium exists only ifR0>1. Therefore, using center manifold theory stated in [15, 16], Figure 1, and stability theory of differential equations, the endemic equilibrium of model (27) is locally asymptotically stable if R0>1 and unstable if R0<1.


### 3.6. Sensitivity Analysis

In this subsection, we determine the most sensitive parameter that contributes in increasing (decreasing) the reproduction number using the normalized forward sensitivity index stated in [17]. That is, the normalized forward sensitive index of R0 with respect to parameter p is denoted by YpR0 and defined as(31)YpR0=∂R0∂ppR0.From our earlier computations,R0 is given by(32)R0=β1ρ+μ+β2ρρ+μη+μ+β3ηρρ+μη+μδ+μ.To determine the most sensitive parameter in the model, we perform the following computations and results are given in Table1:(33)Yβ1R0=∂R0∂β1β1R0=β1ρ+μR0,Yβ2R0=∂R0∂β2β2R0=ρβ2ρ+μη+μR0,Yβ3R0=∂R0∂β3β3R0=ηρβ3ρ+μη+μδ+μR0,YρR0=∂R0∂ρρR0=−β1ρ+μ2+β2μρ+μ2η+μ+β3ημρ+μ2η+μδ+μρR0,YηR0=∂R0∂ηηR0=−β2ρη+μ2ρ+μ+β3ρμη+μ2ρ+μδ+μηR0,YμR0=∂R0∂μμR0=−β1ρ+μ2−β2ρη+ρ+2μρ+μ2η+μ2−β3ηρη+δ+2μρ+μ+η+μδ+μρ+μ2η+μ2δ+μ2μR0,YδR0=∂R0∂δδR0=−β3ηρδρ+μη+μδ+μ2R0.Table 1
Computed sensitivity index value.
ParameterValueSensitivity indexβ10.2+0.23β20.16+0.768β30.0024+0.001ρ0.5−0.19η0.1−1.28μ0.02−0.011δ1−0.001


## 3.1. Well-Posedness of the Model

### 3.1.1. Invariant Region

Theorem 1.
The setΩ⊂ℜ+4 is the invariant region of boundedness if for all initial solutions S0,W0,I0,A0∈Ω and for all t≥0, then St,Wt,It,At∈Ω. That is,(6)Ω=St,Wt,It,At:St,Wt,It,At∈ℜ+4,∀t≥0.Proof.
Considering model (4) and adding all corresponding terms on the left of equality and terms on right of equality, we get(7)dNdt≤λ−μN+δA.
Solving the preceding inequality and considering the expression as timet gets larger, we obtain(8)Nt≤λμ.
The preceding inequality shows that all solution variables are bounded inℜ+4. Therefore, using equation (4), Ω⊂ℜ+4 is the invariant set of solutions such that(9)Ω=St,Wt,It,At:St,+Wt+It+At≤λμ,∀t≥0.


### 3.1.2. Positivity of Solutions

Theorem 2.
All solutions of model (4) are nonnegative for all time t provided that initial conditions are nonnegative.Proof.
Taking the first equation of model (4), we have(10)dSdt=λ−Sβ1W+β2I+β3AN−μS.
Neglecting the termλ, the foregoing equation is reduced to the next inequality(11)dSS≥−β1W+β2I+β3AN+μdt.
Solving the foregoing inequality over time interval0,t, we get(12)St≥S0e−μt−∫0tβ1Wξ+β2Iξ+β3Aξ/Nξdξ.
Since the initial conditionS0 and the exponential expression e−μt−∫0tβ1Wξ+β2Iξ+β3Aξ/Nξdξ are nonnegative, the solution variable St is nonnegative for all time t. Similarly, other solution variables are nonnegative.


### 3.1.3. Existence and Uniqueness of Solutions

Theorem 3.
The solutions of model (4) exist and are unique in ℜ+4.Proof.
Since the expression on the right-hand side of model (4) is bounded and continuously differentiable, the proof directly follows from Cauchy–Lipschitz theorem (see [10]).
Therefore, the formulated model (4) is mathematically well-posed and biologically acceptable.


## 3.1.1. Invariant Region

Theorem 1.
The setΩ⊂ℜ+4 is the invariant region of boundedness if for all initial solutions S0,W0,I0,A0∈Ω and for all t≥0, then St,Wt,It,At∈Ω. That is,(6)Ω=St,Wt,It,At:St,Wt,It,At∈ℜ+4,∀t≥0.Proof.
Considering model (4) and adding all corresponding terms on the left of equality and terms on right of equality, we get(7)dNdt≤λ−μN+δA.
Solving the preceding inequality and considering the expression as timet gets larger, we obtain(8)Nt≤λμ.
The preceding inequality shows that all solution variables are bounded inℜ+4. Therefore, using equation (4), Ω⊂ℜ+4 is the invariant set of solutions such that(9)Ω=St,Wt,It,At:St,+Wt+It+At≤λμ,∀t≥0.


## 3.1.2. Positivity of Solutions

Theorem 2.
All solutions of model (4) are nonnegative for all time t provided that initial conditions are nonnegative.Proof.
Taking the first equation of model (4), we have(10)dSdt=λ−Sβ1W+β2I+β3AN−μS.
Neglecting the termλ, the foregoing equation is reduced to the next inequality(11)dSS≥−β1W+β2I+β3AN+μdt.
Solving the foregoing inequality over time interval0,t, we get(12)St≥S0e−μt−∫0tβ1Wξ+β2Iξ+β3Aξ/Nξdξ.
Since the initial conditionS0 and the exponential expression e−μt−∫0tβ1Wξ+β2Iξ+β3Aξ/Nξdξ are nonnegative, the solution variable St is nonnegative for all time t. Similarly, other solution variables are nonnegative.


## 3.1.3. Existence and Uniqueness of Solutions

Theorem 3.
The solutions of model (4) exist and are unique in ℜ+4.Proof.
Since the expression on the right-hand side of model (4) is bounded and continuously differentiable, the proof directly follows from Cauchy–Lipschitz theorem (see [10]).
Therefore, the formulated model (4) is mathematically well-posed and biologically acceptable.


## 3.2. Equilibriums of the Model

### 3.2.1. Disease-Free Equilibrium (DFE)

The disease-free equilibriumE0 is a point in the system where there is no disease in the population. The computed disease-free equilibrium of the model one is given by(13)E0=S0,0,0,0,where S0=λ/μ.


### 3.2.2. Endemic Equilibrium (EE)

The endemic equilibrium of the model is the point where the disease persists in the population. The computed endemic equilibrium of model (4) is given by(14)E1=S∗,W∗,I∗,A∗,where S∗=a+b+cλ/R0+μa+b+c−1,W∗=λ/ρ+μR0−1/R0+μa+b+c−1,I∗=ρW∗/η+μ,A∗=ηI∗/δ+μ,a=1/ρ+μ,b=ρ/ρ+μη+μ,andc=ηρ/ρ+μη+μδ+μ.


## 3.2.1. Disease-Free Equilibrium (DFE)

The disease-free equilibriumE0 is a point in the system where there is no disease in the population. The computed disease-free equilibrium of the model one is given by(13)E0=S0,0,0,0,where S0=λ/μ.


## 3.2.2. Endemic Equilibrium (EE)

The endemic equilibrium of the model is the point where the disease persists in the population. The computed endemic equilibrium of model (4) is given by(14)E1=S∗,W∗,I∗,A∗,where S∗=a+b+cλ/R0+μa+b+c−1,W∗=λ/ρ+μR0−1/R0+μa+b+c−1,I∗=ρW∗/η+μ,A∗=ηI∗/δ+μ,a=1/ρ+μ,b=ρ/ρ+μη+μ,andc=ηρ/ρ+μη+μδ+μ.


## 3.3. Basic Reproduction NumberR0

The basic reproduction number is the average number of infected individuals produced by one infectious individual in the susceptible population during entire period of infection [11]. We have computed the basic reproduction number using the next-generation matrix employed in [12]. Accordingly, from model (4), let f=Sβ1W+β2I+β3A/N00 represent the newly infected individuals in the infectious compartments and v=ρ+μW−ρW+η+μI−ηI+δ+μA be the remaining terms of infected compartments. So that, the Jacobian matrices obtained from f and v at the disease-free equilibrium are denoted by F and V, respectively. So that,(15)F=β1β2β3000000,V=ρ+μ00−ρη+μ00−ηδ+μ.The next-generation matrixFV−1 is computed as follows:(16)FV−1=β1ρ+μ+β2ρρ+μη+μ+β3ηρρ+μη+μδ+μβ2η+μ+β3ηη+μδ+μβ3δ+μ000000.The computed eigenvalues of the foregoing next-generation matrix are given by(17)λ1=β1ρ+μ+β2ρρ+μη+μ+β3ηρρ+μη+μδ+μ,λ2=0,λ3=0.The basic reproduction numberR0=ρFV−1 and given by(18)R0=β1ρ+μ+β2ρρ+μη+μ+β3ηρρ+μη+μδ+μ.


## 3.4. Stability Analysis of Disease-Free Equilibrium

Theorem 4.
The disease-free equilibrium becomes locally asymptotically stable ifR0<1 and becomes unstable if R0>1.Proof.
To determine the local stability of the disease-free equilibrium, we constructed a JacobianJ from model (2) at the disease-free equilibrium as follows:(19)J=−μ−β1−β2−β30β1−ρ+μβ2β30ρ−η+μ000η−δ+μ.
The eigenvalues of the above matrix can be obtained by solving the characteristic equation as follows:(20)detJ−λI=0.
The characteristic equation can be written as(21)λ+μa0λ3+a1λ2+a2λ+a3=0,where a0=1,a1=η+ρ+δ+3μ−β1,a2=ηδ+2ημ+2δμ+ηη+δ+2μ+3μ2−β1η+δ+2μ−β2,anda3=ημ2+δμ2+ρμ2+μ3+ηδμ+ηδρ+ημρ+δμρ−β1μ2+ηδ+ημ+δμ−β2δρ+μρ−β3ηρ.
By the Hurwitz–Routh principle, all eigenvalues of the characteristic equation are negative if and only if the following conditions are satisfied:(22)a0>0,a1a2>a3.
Therefore, by stability theory of differential equations, the disease-free equilibrium becomes locally asymptotically stable if0<R0<1 and becomes unstable if R0>1.Theorem 5.
The disease-free equilibrium of the constructed model (4) is globally asymptotically unstable if 0<R0<1 and unstable otherwise.Proof.
To prove global stability of the disease-free equilibrium, we apply the methods applied in [13] as follows.
Accordingly, letX∈ℝ1 be uninfected individuals in class S and Y∈ℝ3 be individuals in infected class I,U,A. Hence, we write model (4) as follows:(23)dXdt=HX,Y,dYdt=GX,Y,GX,0=0.
The disease-free equilibrium of the foregoing system is computed as(24)E01=X0,0.
To guarantee global stability, the method we applied must meet conditions H1 and H2 stated as follows:H1. For dX/dt=HX,0,X0 is globally asymptotically stableH2. GX,Y=PY−G^X,Y,G^X,Y≥0 for X,Y∈Ω
Here,P=DYGX,0 be the Metzler matrix and Ω is a region where solutions are acceptable.Theorem 6.
The disease-free equilibrium of model (4) is globally asymptotically stable in a region Ω if R0<1 and unstable whenever R0>1 provided that the stated two conditions H1 and H2 are satisfied.Proof.
Considering model (4), we have(25)HX,0=λ−μS=HS,0.
PuttingHX,0=0 and solving, we obtain S=λ/μ. Hence, X0=λ/μ,0.
We observe thatX0 is the globally asymptotically stable equilibrium of equation as follows:(26)dXdt=HX,0.
Further, considering the infected compartments of model (27), we have(27)dSdt=λ−SNβ1W+β2I+β3A−μS,dWdt=SNβ1W+β2I+β3A−ρ+μW,dIdt=ρW−η+μI,dAdt=ηI−δ+μA,GX,Y=β1W+β2I+β3A−ρ+μ00ρ−η+μ00η−δ+μWIA−β1W+β2I+β3A−SNβ1W+β2I+β3A00.
Further, at the disease-free equilibrium, the preceding equation is reduced to the form as follows:(28)GX,Y=−ρ+μ00ρ−η+μ00η−δ+μWIA−000.
Now, comparing the required condition and the preceding equation, we obtain(29)P=−ρ+μ00ρ−η+μ00η−δ+μ,Y=WIA,G^X,Y=000.
Now, we observe thatGX,Y=PY−G^X,Y, where G^X,Y≥0,∀X,Y.
Therefore, the disease-free equilibriumE0 is globally asymptotically stable if R0<1.


## 3.5. Bifurcation and Stability Analysis of Endemic Equilibrium

In this subsection, we simulate the endemic equilibrium versus reproduction number to identify the kind of bifurcation that occurs at the pointR0=1 and to determine the stability behavior of equilibriums if R0<1 and R0>1. The simulation of the endemic equilibrium S∗,W∗,I∗,A∗ versus R0 is given in Figure 1.Figure 1
Simulation of the endemic equilibrium versus reproduction number.Theorem 7.
The endemic equilibrium of model (2) undergoes forward bifurcation at R0=1.Proof.
(see [3, 14] and Figure 1).Theorem 8.
The endemic equilibriumE1 is globally asymptotically positively stable if R0>1.Proof.
(see [3, 14] and Figure 1).Theorem 9.
The endemic equilibrium of model (27) is locally asymptotically stable if R0>1 and unstable if 0<R0<1.Proof.
We have already computed the endemic equilibriumE1 as(30)E1=S∗,W∗,I∗,A∗,where S∗=a+b+cλ/R0+μa+b+c−1,W∗=λ/ρ+μR0−1/R0+μa+b+c−1,I∗=ρW∗/η+μ,A∗=ηI∗/δ+μ,a=1/ρ+μ,b=ρ/ρ+μη+μ,andc=ρη/ρ+μη+μδ+μ.
Clearly, the endemic equilibrium exists only ifR0>1. Therefore, using center manifold theory stated in [15, 16], Figure 1, and stability theory of differential equations, the endemic equilibrium of model (27) is locally asymptotically stable if R0>1 and unstable if R0<1.


## 3.6. Sensitivity Analysis

In this subsection, we determine the most sensitive parameter that contributes in increasing (decreasing) the reproduction number using the normalized forward sensitivity index stated in [17]. That is, the normalized forward sensitive index of R0 with respect to parameter p is denoted by YpR0 and defined as(31)YpR0=∂R0∂ppR0.From our earlier computations,R0 is given by(32)R0=β1ρ+μ+β2ρρ+μη+μ+β3ηρρ+μη+μδ+μ.To determine the most sensitive parameter in the model, we perform the following computations and results are given in Table1:(33)Yβ1R0=∂R0∂β1β1R0=β1ρ+μR0,Yβ2R0=∂R0∂β2β2R0=ρβ2ρ+μη+μR0,Yβ3R0=∂R0∂β3β3R0=ηρβ3ρ+μη+μδ+μR0,YρR0=∂R0∂ρρR0=−β1ρ+μ2+β2μρ+μ2η+μ+β3ημρ+μ2η+μδ+μρR0,YηR0=∂R0∂ηηR0=−β2ρη+μ2ρ+μ+β3ρμη+μ2ρ+μδ+μηR0,YμR0=∂R0∂μμR0=−β1ρ+μ2−β2ρη+ρ+2μρ+μ2η+μ2−β3ηρη+δ+2μρ+μ+η+μδ+μρ+μ2η+μ2δ+μ2μR0,YδR0=∂R0∂δδR0=−β3ηρδρ+μη+μδ+μ2R0.Table 1
Computed sensitivity index value.
ParameterValueSensitivity indexβ10.2+0.23β20.16+0.768β30.0024+0.001ρ0.5−0.19η0.1−1.28μ0.02−0.011δ1−0.001


## 4. Optimal Control of HIV Model

In this section, we analyze the optimal control model (2) to identify the best control strategy that reduces the number of AIDS individuals and minimizes the total cost used in during the interventions. The incorporated control strategies are described as follows:(i)
Preventive Controlu1. HIV prevention measures include HIV-AIDS education, condom usage, abstinence, and faithfulness.(ii)
Treatment Controlu2. ART treatment reduces the number of virus in human blood, and as a result, the number of individuals that progress to AIDS stage will be reduced. Our goal is to minimize the number of AIDS individuals by effectively using the control measures along with minimum cost. The normalized optimal control model is given by(34)dsdt=b1−s−1−u1sβ1w+β2i+β3a+sδa,dwdt=1−u1sβ1w+β2i+β3a−ρ+bw+wδa,didt=ρw−1−u2η+bi+iδa,dadt=1−u2ηi−b+δa+δa2.To characterize the optimal levels of the controls, we define the lebesgue measurable control setU as (35)U=u1,u2:0≤u1≤1,0≤u2≤1,0≤t≤T.The 0 value of control measure means that there are no control measures taken against HIV/AIDS prevention or treatment [18]. The value of control measure indicates that full control measure is taken against the HIV/AIDS prevention or treatment. The objective functional or cost functional that describes the consumed total cost to minimize the number of window individuals and AIDS individuals is defined as(36)J=min∫0tfc1w+c2a+12w1u12+w2u22dt,where the coefficients c1 and c2 of state variables are constants and coefficients w1 and w2 are the measure of the consumed costs of intervention associated to controls u1 and u2, respectively [7]. All constant coefficients are positive [18]. Our aim is to obtain the optimal controls that minimize new infected individuals and reduces the number of AIDS individuals along with minimum cost of interventions. Since cost is nonlinear, we assume quadratic expression 1/2wiui2 where i=1,2.


### 4.1. Existence and Description of Optimal Control Solution

Theorem 10 (existence of optimal solution).
There exist optimal controlsu1∗andu2∗ and state variables s∗,w∗,i∗, and a∗ for the initial value problems (34) and (39) that minimize Ju1,u2 over U.Proof.
The admissible control setU is determined using Fleming and Rishel’s theorem.(i)
The solution set of (34) and (39) with corresponding control function be nonempty(ii)
The state system is a linear function of control variables and coefficients depending on state variables and time(iii)
The integrandL=c1w+c2a+1/2w1u12+w2u22 of objective functional J is convex on U and L≥αu1,u2β−γ, where α>0 and β>1
To verify condition (i), we use [7]. If the functions of state equations are continuous, bounded, and Lipschitz in state variables, then to every admissible control U, there is unique solution. The total population size is bounded, the state variables are bounded, and the partial derivatives of functions constructed from state equations are bounded. Therefore, condition (i) is satisfied. Condition (ii) holds by observing the linearity of the state equation in controls u1 and u2. Condition (iii) follows from definition that any quadratic function, linear, and constant functions are convex. Hence, L is convex on U. Next, we show boundedness of L.
Sincew1u22≤w2 and 1/2w1u22≤1/2w22 for u2∈0,1, then L=c1w+c2a+1/2w1u12+w2u22≥1/2w1u12+w2u22−1/2w22.(37)⇒L≥minw12,w22u12+u22−12w22,⇒L≥minw12,w22u1,u22−12w22.
Thus,L≥αu1,u2β−γ, where α=minw1/2,w2/2,β=2,andγ=1/2w22.
Applying Pontryagin’s Maximum Principle (PMP) as in [7, 19], we obtained the necessary conditions to be fulfilled by optimal control pairs. Therefore, in order to construct the optimal system that minimize the cost functional, we defined the Hamiltonian function as follows:(38)H=c1w+c2a+12w1u12+w2u22+λ1b1−s−1−u1sβ1w+β2i+β3a+sδa+λ21−u1sβ1w+β2i+β3a−ρ+bw+wδa+λ3ρw−1−u2η+bi+iδa+λ41−u2ηi−b+δa+δa2,where λi,i=1,2,3,4 are adjoint variable functions corresponding to state variables s,w,i, and a determined as follows [4] for existence of optimal control pairs.
Variablesλi,i=1,2,3,4 such that(39)dλ1dt=−∂H∂s=1−u1β1w+β2i+β3aλ1−λ2+b−δaλ1,dλ2dt=−∂H∂w=−c1+1−u1sβ1λ1−λ2+ρλ2−λ3+b−δaλ2,dλ3dt=−∂H∂i=1−u1sβ2λ1−λ2+1−u2ηλ3−λ4+b−δaλ3,dλ4dt=−∂H∂a=−c2+1−u1sβ3λ1−λ2−sδλ1−wδλ2+iδλ3+b+δ−2δaλ4.
With transversality conditionsλitf=0,i=1,2,3,4. Further, we obtained the control set u1∗t,u2∗t such that u1∗t=max0,min1,u1∗,u2∗t=max0,min1,u2∗, where u1∗=sβ1w+β2i+β3aλ2−λ1/w1andu2∗=ηiλ4−λ3/w2.


## 4.1. Existence and Description of Optimal Control Solution

Theorem 10 (existence of optimal solution).
There exist optimal controlsu1∗andu2∗ and state variables s∗,w∗,i∗, and a∗ for the initial value problems (34) and (39) that minimize Ju1,u2 over U.Proof.
The admissible control setU is determined using Fleming and Rishel’s theorem.(i)
The solution set of (34) and (39) with corresponding control function be nonempty(ii)
The state system is a linear function of control variables and coefficients depending on state variables and time(iii)
The integrandL=c1w+c2a+1/2w1u12+w2u22 of objective functional J is convex on U and L≥αu1,u2β−γ, where α>0 and β>1
To verify condition (i), we use [7]. If the functions of state equations are continuous, bounded, and Lipschitz in state variables, then to every admissible control U, there is unique solution. The total population size is bounded, the state variables are bounded, and the partial derivatives of functions constructed from state equations are bounded. Therefore, condition (i) is satisfied. Condition (ii) holds by observing the linearity of the state equation in controls u1 and u2. Condition (iii) follows from definition that any quadratic function, linear, and constant functions are convex. Hence, L is convex on U. Next, we show boundedness of L.
Sincew1u22≤w2 and 1/2w1u22≤1/2w22 for u2∈0,1, then L=c1w+c2a+1/2w1u12+w2u22≥1/2w1u12+w2u22−1/2w22.(37)⇒L≥minw12,w22u12+u22−12w22,⇒L≥minw12,w22u1,u22−12w22.
Thus,L≥αu1,u2β−γ, where α=minw1/2,w2/2,β=2,andγ=1/2w22.
Applying Pontryagin’s Maximum Principle (PMP) as in [7, 19], we obtained the necessary conditions to be fulfilled by optimal control pairs. Therefore, in order to construct the optimal system that minimize the cost functional, we defined the Hamiltonian function as follows:(38)H=c1w+c2a+12w1u12+w2u22+λ1b1−s−1−u1sβ1w+β2i+β3a+sδa+λ21−u1sβ1w+β2i+β3a−ρ+bw+wδa+λ3ρw−1−u2η+bi+iδa+λ41−u2ηi−b+δa+δa2,where λi,i=1,2,3,4 are adjoint variable functions corresponding to state variables s,w,i, and a determined as follows [4] for existence of optimal control pairs.
Variablesλi,i=1,2,3,4 such that(39)dλ1dt=−∂H∂s=1−u1β1w+β2i+β3aλ1−λ2+b−δaλ1,dλ2dt=−∂H∂w=−c1+1−u1sβ1λ1−λ2+ρλ2−λ3+b−δaλ2,dλ3dt=−∂H∂i=1−u1sβ2λ1−λ2+1−u2ηλ3−λ4+b−δaλ3,dλ4dt=−∂H∂a=−c2+1−u1sβ3λ1−λ2−sδλ1−wδλ2+iδλ3+b+δ−2δaλ4.
With transversality conditionsλitf=0,i=1,2,3,4. Further, we obtained the control set u1∗t,u2∗t such that u1∗t=max0,min1,u1∗,u2∗t=max0,min1,u2∗, where u1∗=sβ1w+β2i+β3aλ2−λ1/w1andu2∗=ηiλ4−λ3/w2.


## 5. Numerical Simulations of Optimal Control Solutions

Numerical simulations are very important to describe the behavior of populations qualitatively. Numerical solution of optimal control problem is simulated by assuming parameter values logically as described in Table2.Table 2
Parameter and value assigned.
ParameterValueSourceλ800000Assumedβ10.2Assumedβ20.16Assumedβ30.0024Assumedη0.1[18, 9]δ1[18, 9]ρ0.5Assumedb0.03Assumedμ0.02[9]So far, we have discussed the optimal control problem analytically. In this section, to validate our analytical findings, we perform numerical simulations to analyze the minimum cost required during application of control measures. So, we solve control systems (34) and (39) numerically by using the value of parameters assigned in Table 2 and initial population sizes along with fixed final time tf=10 years. Assume the initial population size in the proportionality form as s=0.6,w=0.1,i=0.2,anda=0.1. Also, we assumed constants in objective functional as c1=2,c2=3,w1=10,andw2=20.


### 5.1. Cost-Effectiveness Analysis

In this section, we analyze cost-effectiveness of control strategies to rank them in terms of their cost. To identify the best control strategies along with minimum coast, we apply the method applied in other researchers work in [20]. They had stated a formula as follows:(40)incremental cost−effectiveness ratioICER=difference in costs between strategiesdifference in health effects between strategies.Based on the description of controls, infectious averted, and total cost given in Table3, we compute the incremental cost-effectiveness ratio (ICER) as(41)ICER1=2.64700.00975=271.4872,(42)ICER2=2.6931−2.64700.41184−0.00975=0.1147.Table 3
Description of controls, infectious averted, and cost.
StrategiesDescriptionTotal infectious avertedTotal cost (USD)1Prevention control0.009752.64702Treatment control0.411842.69313Prevention and treatment controls0.008792.644ComparingICER2<ICER1, thus select strategy (2) and ignore strategy (1) as it is expensive and less effective. Again, we compare ICER (2) and ICER (3) by first computing ICER (2) and ICER (3) as follows:(43)ICER2=2.69310.41184=6.5392,(44)ICER3=2.644−2.69310.00879−0.41184=0.1218.Comparing the computed valuesICER3<ICER2 shows that strategy (2) is more expensive and less effective than strategy (3). Therefore, we conclude that using both preventive and treatment control strategies is more effective than using controls separately toward prevention and controlling of the human immunodeficiency virus.


## 5.1. Cost-Effectiveness Analysis

In this section, we analyze cost-effectiveness of control strategies to rank them in terms of their cost. To identify the best control strategies along with minimum coast, we apply the method applied in other researchers work in [20]. They had stated a formula as follows:(40)incremental cost−effectiveness ratioICER=difference in costs between strategiesdifference in health effects between strategies.Based on the description of controls, infectious averted, and total cost given in Table3, we compute the incremental cost-effectiveness ratio (ICER) as(41)ICER1=2.64700.00975=271.4872,(42)ICER2=2.6931−2.64700.41184−0.00975=0.1147.Table 3
Description of controls, infectious averted, and cost.
StrategiesDescriptionTotal infectious avertedTotal cost (USD)1Prevention control0.009752.64702Treatment control0.411842.69313Prevention and treatment controls0.008792.644ComparingICER2<ICER1, thus select strategy (2) and ignore strategy (1) as it is expensive and less effective. Again, we compare ICER (2) and ICER (3) by first computing ICER (2) and ICER (3) as follows:(43)ICER2=2.69310.41184=6.5392,(44)ICER3=2.644−2.69310.00879−0.41184=0.1218.Comparing the computed valuesICER3<ICER2 shows that strategy (2) is more expensive and less effective than strategy (3). Therefore, we conclude that using both preventive and treatment control strategies is more effective than using controls separately toward prevention and controlling of the human immunodeficiency virus.


## 6. Results and Discussion

The optimum control model of HIV/AIDS is numerically solved in this paper using Runge–Kutta fourth-order forward and backward methods, as well as MATLAB software. Simulation in Figure1 shows that forward bifurcation occurs at R0=1 and the endemic equilibrium exists only if R0>1. Thus, the disease persists in the population if and only if the average number of infected individuals produced by an infectious individual in the entire period of infection is greater than one. Also, the endemic equilibrium is globally asymptotically stable if R0>1 and unstable if R0<1. Figure 2 shows the simulation of the infected individuals equilibrium versus reproduction number; the size of all HIV-infected individuals becomes zero at R0=1, and HIV-infected individuals survive if R0>1. A reproduction number’s sensitivity to highly sensitive parameters is investigated in Figure 3. The value of the reproduction number drops as we raise the parameter with a negative sensitive index and increases as we increase the parameter with a positive sensitivity index (see, Table 1). The inclusion of preventative and treatment control measures spared more susceptible persons from becoming infected with the human immunodeficiency virus, as shown in Figure 4. The presence of control measures lowered the number of window persons, as shown in Figure 5. The existence of preventative and treatment control measures lowered the number of HIV patients who were transitioned to acquired immunodeficiency syndrome (AIDS), as shown in Figure 6. Figure 7 illustrates that the availability of preventative and treatment control measures has no discernible effect on the number of people living with AIDS. Figure 8 indicates that the preventative control measure falls from the start to the end of the time, whereas the treatment control measure climbs for about 5 years before gradually reducing until it reaches zero at the end. The adjoint variables function satisfies the transversal criteria necessary in the optimal control problem, as shown in Figure 9. Furthermore, the window period stage of HIV is a narrow stage that might open a gate for HIV transmission in the population without effective preventative control efforts, but the other stages of HIV development can be controlled with good preventive and treatment control interventions.Figure 2
The simulations of HIV-infected populations versus the reproduction number.Figure 3
The simulation of the reproduction number versus parametersβ2 and η.Figure 4
Simulation of susceptible population size without and with controls.Figure 5
Simulation of window population size without and with controls.Figure 6
Simulation of HIV population size without and with controls.Figure 7
Simulation of AIDS population size without and with controls.Figure 8
Simulation of preventive and treatment controls.Figure 9
Simulation of susceptible population size without and with controls.


## 7. Conclusion

In this study, a new method is applied to determine a kind of bifurcation atR0=1, and the stability analysis is carried out if R0<1 and R0>1. The current study finding shows that the endemic equilibrium of model (4) exhibits forward bifurcation at R0=1. The transmission of HIV persists in the population if R0>1 and extincts if R0<1. According to the findings of the sensitivity analysis, an increase in the rate of HIV transmission from HIV persons contributes to the prevalence of HIV in the community. When both preventative and treatment control strategies are employed combined, they are more effective in limiting the transmission dynamics of the human immunodeficiency virus. Moreover, the cost-effectiveness analysis shows that the control measures are more economically consumed if both controls are applied together in controlling the human immunodeficiency virus.



---

*Source: 7471290-2021-12-20.xml*