Backward bifurcation on HIV/AIDS SEI₁I₂TAR model with multiple interactions between sub-populations

Abstract

The HIV/AIDS model was dynamically analyzed in this study. The model has seven compartments: the uneducated, the educated, the HIV-positive who take antiretroviral therapy (ART), the HIV-positive who do not take ART, people receiving ART treatment, people with AIDS who do not receive any treatment (full-blown AIDS), and the recovered. This model takes into account the analysis of the multiple interactions between the uneducated and the educated subpopulations, the HIV-positive who take and who do not take ART. The free-disease and endemic equilibrium points, as well as the basic reproduction number $(R_0)$ as a limit condition for infection-free and endemic occurrence, were produced by a mathematical analysis. The center-manifold hypothesis was used to prove that a backward bifurcation exists. The free-disease and endemic equilibrium points coexist when $R_0<1.$ This means that HIV/AIDS is still spreading. A basic reproduction number below one is insufficient to constitute a free-disease condition. In order to determine essential parameters that significantly contribute to HIV/AIDS transmission, we computed sensitivity index values using a sensitivity analysis. The HIV/AIDS model and bifurcation parameter both identified the rate of HIV transmission from uneducated individuals to HIV-positive individuals who do not receive ART as the most crucial parameter. A numerical simulation supports the dynamical analysis.

Keywords:

• HIV/AIDS model
• dynamical system
• multiple interactions
• backward bifurcation

1. Introduction

A disease caused by the human immunodeficiency virus (HIV) is called HIV/AIDS (HIV/Acquired Immunodeficiency Syndrome). A person with an early HIV infection may not notice symptoms like a healthy person. As the infection progresses, it further damages the immune system. This infection’s ultimate symptom is acquired immunodeficiency syndrome (AIDS). HIV is mainly spread through unprotected sex, sullied blood transfusions, needles, pregnancy, mother-to-child transmission, childbirth, and breastfeeding.

Researchers proposed mathematical formulas to understand the spread of HIV/AIDS. The first model was constructed by Greenhalgh and Hay (1997) with developed Kaplan’s model (1989). They used a dynamical analysis to simulate the spread of HIV/AIDS in a community of injectable drug users. Cai, Li, Ghosh, and Guo (2009) investigated HIV/AIDS with a treatment in the SIAJ (Susceptible, Infected, AIDS, and Symptomatic) model. Huo, Chen, and Wang (2016) proposed the SIATR (susceptible, infected, AIDS, treatment, and recovered) model and studied its dynamical analysis. As the traditional model, Naik, Zu, and Owolabi (2021) created and examined a model of viral kinetics for HIV-1 primary infection in the presence of immune control and treatment. The literature has examined mathematical models in epidemiology to gain a better understanding of the processes underlying cancers connected to AIDS, giving us insight into viral oncogenesis and cancer immunity (Naik, Owolabi, Yavuz, & Zu, 2020).

Antiretroviral medication (ART) should be provided and started for those who have been diagnosed with HIV as soon as feasible following the diagnosis. Patients should also have routine clinical and laboratory evaluations, including tests to check for the presence of the virus in the blood (viral load). Regular use of ART prevents the spread of HIV to others (World Health Organization, 2023). By enhancing human abilities to process and comprehend the risks connected with the HIV/AIDS epidemic, education plays a significant role in enhancing population health. Insufficient knowledge makes it difficult for people to consider the potential health risks of their activities. A risk factor for the transmission of sexually transmitted illnesses is education. On the other hand, the incorrect portrayal of HIV/AIDS as a homosexual disease has increased the number of heterosexuals engaging in unsafe sexual conduct without taking transmission risk into account. Therefore, education becomes an essential factor to consider in the model. The HIV/AIDS model was developed by including an educated compartment (Ulfa, Trisilowati, & Kusumawinahyu, 2018; Wibowo, Hidayat, & Marsudi, 2016).

The number of new HIV infections is projected to decrease from 1.5 million in 2020 to 335,000 by 2030. The number of HIV-related deaths is also expected to decline, from 680,000 in 2020 to less than 240,000 by 2030. These projections are based on the WHO’s 2022-2030 global health sector strategy on HIV (World Health Organization, 2023). In accordance with the WHO mission, the Ministry of Health of the Republic of Indonesia planned to eradicate the HIV endemic by 2030 (Ministry of Health Republic of Indonesia, 2022). As a manifestation of this commitment, the Ministry of Health is making efforts to tackle HIV/AIDS by taking the 95-95-95 fast track, which refers to meeting the target indicators of 95% estimated people living with HIV (ODHIV) whose HIV status is known, 95% people living with HIV and treated, and 95% people living with HIV who received a treatment and experienced viral suppression. However, according to the 2018-2022 data, the achievement of these targets, especially for women, children, and adolescents, is still not optimal. This is because only 79% of people living with HIV (ODHIV) knew their HIV status, only 41% of people living with HIV were treated, and 16% of people living with HIV who received a treatment experienced viral suppression. In Habibah, Trisilowati, Pradana, and Villadystian (2021), there are 640,633 people living with HIV in June 2018. By 2021, it decreased and was estimated that there would be 526,841 people living with HIV, with an estimated 27 thousand new cases.

According to the above data, we derived a mathematical model of HIV/AIDS with compartments $I_2$ representing people living with HIV (ODHIV) whose HIV status is known, $I_1$ representing people living with HIV who received a treatment and experienced viral suppression, T representing people living with HIV who are treated, and E representing educated individuals.

In the last couple of years, different stages of infection have become interesting topics in the HIV/AIDS model. Huo and Chen (2015), Omondi, Mbogo, and Luboobi (2018), Ulfa et al. (2018), Mushayabasa and Bhunu (2011), Naik, Owolabi, et al. (2020) and Habibah, Trisilowati, Pradana, et al. (2021) studied mathematical analyses in the development of a model of HIV/AIDS infection with distinct stages. The model includes infections from the mildest stage to the most severe. Now, individuals with full-blown AIDS, which is the most serious stage of the infection, who are not receiving any treatment make up A. Habibah, Trisilowati, Muzaqi, Tania, and AlFaruq (2021a), Habibah, Trisilowati, Muzaqi, Tania, and AlFaruq (2021b), and Habibah, Owati, and Muzaqi (2021) constructed an HIV/AID model by taking into account multiple-interactions between sub-populations. We call it the $SE{I}_1{I}_{2}TAR$ model. Here, people who alter their sexual behaviors and stick with this alteration for the remainder of their lives are known as R. The dynamics of the model were analyzed in several cases, but the analyses were based on the assumption that not all interactions were considered. The stability analysis of the model can be determined from the basic reproduction number ($R_0$). According to the findings, when $R_0>1,$ asymptotically, the endemic equilibrium point is stable. Conversely, it is unstable. A dynamic analysis involving the relationship between all compartments remains to be conducted.

In this paper, we analyzed the dynamics of the $SE{I}_1{I}_{2}TAR$ model comprehensively by considering multiple-interactions between all sub-populations, especially the multiple interactions between S, E, $I_1,$ and $I_2$ sub-populations. The resultant model underwent a dynamical analysis. We demonstrate the positivity and bounds of the system’s solutions (Brauer & Chavez, 2016; Murray, 2002). Additionally, we determined the model’s solution, the basic reproduction number $R_0,$ and the disease-free and endemic equilibrium points (Abidemi, Owolabi, & Pindza, 2022; Heffernan, Smith, & Wahl, 2005; Karaagac & Owolabi, 2021; Mishra, Purohit, Owolabi, & Sharma, 2020). A mathematical analysis of the model was derived to determine the free-disease and endemic equilibrium points, existence, and local stability analysis. In this analysis, the reproduction number ($R_0<1$) is considered insufficient to meet the free-disease condition. It means that the infection can continue to spread.

Recently, the study of the backward bifurcation phenomenon has become an interesting problem. A backward bifurcation occurs when stable free-disease and endemic equilibria coexist. It happens when the basic reproduction number is below one. The center-manifold theory is used to analyze backward bifurcations. Castillo-Chavez and Song (2004) were the first to investigate this method. The center-manifold theory was applied to the tuberculosis model. Gumel (2012) studied the causes of backward bifurcations in some epidemiological models, such as the TB model, the dengue-fever model, and the West Nile virus model. There are several papers presenting backward bifurcation analyses, such as those by Buonomo & Vargas-De-León (2013), Gerberrya and Philipa (2016), Martcheva (2015), Liu and Liu (2018), Bi, Chen, Wu, and Ben-Arieh (2020), and Fatimah, Aldila, and Handari (2021). In this paper, we will show and explain a backward bifurcation in the $SE{I}_1{I}_{2}TAR$ model. The bifurcation parameter now plays a crucial function as a limit condition of a backward bifurcation. We calculated the sensitivity index values using a sensitivity analysis to obtain parameters that significantly contribute to HIV/AIDS transmission (Chitnis, Cushing, & Hyman, 2006; Rois, Trisilowati, & Habibah, 2021).

Finally, using the Runge-Kutta method with $O(h^4),$ we numerically solved the model and simulated it with the given parameters to back up the findings of the mathematical analysis.

2. Model analysis

2.1. Model formulation

A mathematical model of HIV/AIDS is proposed, with multiple interactions between uneducated people, educated people, HIV-positive people who take ARV, and HIV-positive people who do not take ARV. The model includes seven sub-populations: susceptible/uneducated people (S), educated people (E), people with HIV who are taking ARV ($I_1$), people with HIV who are not taking ARV ($I_2$), people with AIDS who are not receiving any treatment (A), people who are receiving an ARV treatment (T), and people who have recovered (R). The compartment diagram of the HIV/AIDS $SE{I}_1{I}_{2}TAR$ epidemic model, which includes multiple interactions, is shown in Figure 1.

Figure 1. Diagram of $SE{I}_1{I}_{2}TAR$ formula (Habibah et al., 2021a).

The compartment diagram pertaining to HIV/AIDS is developed and can be found in Habibah et al. (2021a). However, the model constructed did not account for the interaction between the infected people $I_1$ and $I_2$ and the educated subpopulation (E). Since $\left\{{\beta }_3,{\beta }_4\right\}=0$ was assumed, the model is rather straightforward. The dynamic of each sub-population follows the transmission diagram given in Figure 1 and is described as follows. The number of uneducated or susceptible individuals S increases with the rate of recruitment $\Lambda$ and decreases due to interactions with $I_1$ and $I_2.$ This sub-population decreases due to the natural death rate d and due to its receiving HIV/AIDS education at the rate $\eta .$ Therefore, the dynamic of uneducated or susceptible individuals is given by $$\frac{dS}{dt}=\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS.$$

The number of educated individuals (E) grows as susceptible individuals receive education at a rate $\eta .$ The interactions with $I_1$ and $I_2,$ as well as the natural death rate d, cause a decrease in the number of educated individuals (E). The dynamic of educated individuals is as follows. $$\frac{dE}{dt}=\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE.$$

The number of people with HIV who are taking ARV ($I_1$) grows as a result of interactions with uneducated people (S) and educated people (E); S and E now become infected people $I_1.$ A successful treatment contributes to an increase in $I_1$ at a rate ${\alpha }_1.$ On the other hand, the number of people with HIV who are taking ARV ($I_1$) decreases due to the natural death rate d. Individuals with HIV who are taking ARV ($I_1$) will have the worst health, so they require medical treatment at a rate $k_1.$ The dynamic of individuals with HIV who are taking ARV ($I_1$) is given by $$\frac{d{I}_1}{dt}={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1.$$

Following that, the number of people with HIV who are not taking ARV ($I_2$) grows by interactions with uneducated people (S) and educated people (E); S and E now become infected people ($I_2$). On the other hand, the number of people with HIV who are not taking ARV ($I_2$) decreases due to the natural death rate d. Individuals with HIV who are not taking ARV ($I_2$) will have the worst health and develop full-blown AIDS (A) at a rate $k_2,$ so they will require medical treatment at a rate $k_3.$ The dynamic of individuals with HIV who are not taking ARV ($I_2$) is given by $$\frac{d{I}_2}{dt}={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2.$$

The number of people receiving an ARV treatment (T) increases with the progression rates of people with HIV who are taking ARV ($I_1$) and people with HIV who are not taking ARV ($I_2$) who then receive a treatment. A Successful treatment, in turn, contributes to the decrease in T at a rate ${\alpha }_1,$ with the individuals (T) entering the $I_1$ compartment. A treatment failure contributes to a decrease in T at a rate ${\alpha }_2,$ with the individuals entering the A compartment. It means that the individuals receiving treatment develop full-blown AIDS (A). The natural mortality rate d and the disease-related mortality rate for treatment individuals ${\delta }_2$ both contribute to the reduction of treatment individuals. The dynamic of individuals who are receiving ARV treatment (T) is given by $$\frac{dT}{dt}=k_1{I}_1+k_3{I}_2-eT.$$

The number of people with AIDS who are not receiving any treatment (full-blown AIDS) (A) rises in tandem with the progression rate of individuals with HIV who are not taking ARV ($I_2$) and treatment failure at a rate ${\alpha }_2.$ On the other hand, the disease-related death rate ${\alpha }_1$ for AIDS (A) patients and the natural death rate d contribute to the decreasing number of AIDS patients (A). The dynamic of individuals who are not receiving any ARV treatment (A) is given by $$\frac{dA}{dt}=k_2{I}_2+{\alpha }_{2}T-fA.$$

The final compartment is occupied by recovered individuals (R). The R class consists of susceptible individuals (S) who adopt safe sexual habits and stick to them for the rest of their lives at a rate of $\mu .$ This class decreases at the natural death rate d. The dynamic of recovered individuals R is given by $$\frac{dR}{dt}=\mu S-dR.$$

Based on the above description, the mathematical model of the following system of nonlinear differential equations describes HIV/AIDS. (1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) with $a=\mu +\eta +d,$ $b=k_1+d,$ $c=k_2+k_3+d,$ $e={\alpha }_1+{\alpha }_2+{\delta }_2+d,$ and $f={\delta }_1+d.$ The initial values of the variables are $S(0)\ge 0,$ $E(0)\ge 0,$ $I_1(0)\ge 0,$ $I_2(0)\ge 0,$ $T(0)\ge 0,$ $A(0)\ge 0,$ and $R(0)\ge 0.$ All parameters in the system (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) are positive values as given in Table 1.

Table 1. Parameter description of the $SE{I}_1{I}_{2}TAR$ model.

Symbol	Description	Value	Data source
$\Lambda$	The rate of recruitment	0.55	Huo et al., 2016
$\mu$	The rate of S to R	0.03	Huo et al., 2016
${\beta }_1$	Transmission rate from S to $I_1$	0.0023	Habibah, Trisilowati, Pradana, et al., 2021
${\beta }_2$	Rate of transmission from S to $I_2$	0.0033	Habibah, Trisilowati, Pradana, et al., 2021
${\beta }_3$	Rate of transmission from E to $I_1$	0.0013	Assumption
${\beta }_4$	Rate of transmission from E to $I_2$	0.0019	Assumption
d	The human natural death rate	0.0196	Huo et al., 2016
$\eta$	Education rate	0.01	Assumption
$k_1$	Rate of progression from $I_1$ to T	0.0498	Habibah, Trisilowati, Pradana, et al., 2021
$k_2$	Rate of progression from $I_2$ to A	0.008	Habibah, Trisilowati, Pradana, et al., 2021
$k_3$	Rate of progression from $I_2$ to T	0.05	Habibah, Trisilowati, Pradana, et al., 2021
${\alpha }_1$	Fraction of successful treatment	0.02	Habibah, Trisilowati, Pradana, et al., 2021
${\alpha }_2$	Fraction of failure of treatment	0.05	Habibah, Trisilowati, Pradana, et al., 2021
${\delta }_1$	Death rate related to disease for AIDS individuals	0.0909	Huo et al., 2016
${\delta }_2$	Death rate related to disease for Treatment individuals	0.0667	Huo et al., 2016

In the previous research (Habibah et al., 2021a), the interaction between the educated subpopulation (E) and infected individuals $I_1$ and $I_2$ did not consider in the mathematical model and its mathematical stability, so that $\left\{{\beta }_3,{\beta }_4\right\}=0.$ This model presents multiple interactions between uneducated people (S), educated people (E), HIV-positive people who take ARV ($I_1$), and HIV-positive people who do not take ARV ($I_2$). Because education plays a significant role in preventing the spread of HIV to others, the interaction between educated individuals (E) and infected individuals $I_1$ and $I_2$ is less than the interaction between uneducated or susceptible individuals (S) and infected individuals $I_1$ and $I_2.$ Furthermore, the infection rate can be expressed mathematically as $\left\{{\beta }_1,{\beta }_2\right\}>\left\{{\beta }_3,{\beta }_4\right\}.$ On the other hand, people who do not take ARV ($I_2$) have the potential to transmit the HIV virus to a greater extent than people who take ARV ($I_1$). Therefore, we can also say that the infection rate (${\beta }_2$) among susceptible people $I_2$ is greater than the infection rate among susceptible people $I_1,$ or we can put it as ${\beta }_1<{\beta }_2.$ It will be further discussed in the sensitivity analysis in sub-section 2.4.

2.2. Boundedness and positivity of model solution

The boundedness and positivity of model solutions will be proved in the following theorems.

Theorem 2.1

(Boundedness). All feasible solutions $S(t),$ $E(t),$ $I_1(t),$ $I_2(t),$ $A(t),$ $T(t)$ and $R(t)$ of system (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) are bounded by a region $D=\{(S,E,I_1,I_2,A,T,R)\in {\mathfrak{R}}^7:S+E+I_1+I_2+A+T+R\le \Lambda /d\}.$

Proof.

Total population of system equation (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) can be written as (2) $$\begin{array}{cc}\frac{dN}{dt}=\frac{dS}{dt}+\frac{dE}{dt}+\frac{d{I}_1}{dt}+\frac{d{I}_2}{dt}+\frac{dA}{dt}+\frac{dT}{dt}+\frac{dR}{dt},\hfill & \hfill \\ =\Lambda -dN(t)-{\delta }_{1}A-{\delta }_{2}T.\hfill & \hfill \end{array}$$(2)

It implies that (3) $$\frac{dN}{dt}\le \Lambda -dN(t),$$(3) and by using the method to solve an ordinary differential equation, we obtain (4) $$N\le \Lambda /d+N(0)e^{-dt},$$(4) with the initial value $N(0).$ Then, by taking limit $t\to \infty$ we get (5) $$\underset{t\to \infty }{\lim }\sup N(t)\le \Lambda /d.$$(5)

We can write $N=S+E+I_1+I_2+A+T+R\le \Lambda /d.$ The system(1) is well defined in the region $D=\{(S,E,I_1,I_2,A,T,R)\in {\mathfrak{R}}^7:S+E+I_1+I_2+A+T+R\le \Lambda /d\}.$ Furthermore, we prove the system (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) in region D has solutions with a non-negative sign. □

Theorem 2.2

(Positivity). If $S(0)\ge 0,$ $I_1(0)\ge 0,$ $I_2(0)\ge 0,$ $A(0)\ge 0,$ $T(0)\ge 0$, and $R(0)\ge 0$, then the solution of system (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) $S(t),$ $I_1(t),$ $I_2(t),$ $A(t),$ $T(t)$ and $R(t)$ are positive for all $t>0.$

Proof.

We have from the first equation of (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) (6) $$\begin{array}{cc}\frac{dS}{dt}=\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill & \hfill \\ =\Lambda -S[{\beta }_1{I}_1-{\beta }_2{I}_2-a],\hfill & \hfill \\ =\Lambda -Q_1(t)S,\hfill & \hfill \end{array}$$(6) where $Q_1(t)={\beta }_1{I}_1-{\beta }_2{I}_2-a.$ Equation (2.2) is multiplied by $e^{{\int }_0^t{Q}_1(r)dr}$ to result (7) $$\frac{dS}{dt}{e}^{{\int }_0^t{Q}_1(r)dr}=\{\Lambda -Q_1(t)S(t)\}{e}^{{\int }_0^t{Q}_1(r)dr},$$(7) which implies (8) $$\frac{dS}{dt}{e}^{{\int }_0^t{Q}_1(r)dr}+Q_1(t)S(t)e^{{\int }_0^t{Q}_1(r)dr}=\Lambda e^{{\int }_0^t{Q}_1(r)dr}.$$(8)

Additionally, the derivative of $S(t)e^{{\int }_0^t{Q}_1(r)dr}$ with respect to t can be used to express the left side of Equation (8)(8) $$\frac{dS}{dt}{e}^{{\int }_0^t{Q}_1(r)dr}+Q_1(t)S(t)e^{{\int }_0^t{Q}_1(r)dr}=\Lambda e^{{\int }_0^t{Q}_1(r)dr}.$$(8) , to yield (9) $$\frac{d}{dt}\{S(t)e^{{\int }_0^t{Q}_1(r)dr}\}=\Lambda e^{{\int }_0^t{Q}_1(r)dr},$$(9) calculating the integral with respect to q from 0 to t, then (10) $$S(t)e^{{\int }_0^t{Q}_1(r)dr}-S(0)=\Lambda \left\{{\int }_0^t{e}^{{\int }_0^q{Q}_1(r)dr}dq\right\}.$$(10)

Equation (10)(10) $$S(t)e^{{\int }_0^t{Q}_1(r)dr}-S(0)=\Lambda \left\{{\int }_0^t{e}^{{\int }_0^q{Q}_1(r)dr}dq\right\}.$$(10) is multiplied by $e^{-{\int }_0^t{Q}_1(r)dr}$ to get (11) $$S(t)-S(0)e^{-{\int }_0^t{Q}_1(r)dr}=\Lambda e^{-{\int }_0^t{Q}_1(r)dr}\left\{{\int }_0^t{e}^{{\int }_0^q{Q}_1(r)dr}dq\right\}.$$(11)

Eventually, we achieve (12) $$S(t)=S(0)e^{-{\int }_0^t{Q}_1(r)dr}+\Lambda e^{-{\int }_0^t{Q}_1(r)dr}\left\{{\int }_0^t{e}^{{\int }_0^q{Q}_1(r)dr}dq\right\}\ge 0,$$(12) signifies that the system (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) solution for $S(t)$ is positive The solution for $E(t),$ $I_1,$ $I_2,$ T, A, and R can be obtained with the same procedures as follows (13) $$E(t)=E(0)e^{-{\int }_0^t{Q}_2(r)dr}+\eta e^{-{\int }_0^t{Q}_2(r)dr}\{{\int }_0^{t}E(t)e^{{\int }_0^q{Q}_2(r)dr}dq\}\ge 0,$$(13) (14) $$I_1(t)=I_1(0)e^{-{\int }_0^t{Q}_3(r)dr}+{\alpha }_1{e}^{-{\int }_0^t{Q}_3(r)dr}\{{\int }_0^{t}T(t)e^{{\int }_0^q{Q}_3(r)dr}dq\}\ge 0,$$(14) (15) $$I_2(t)=I_2(0)e^{{\int }_0^t{Q}_4(r)dr}\ge 0,$$(15) (16) $$T(t)=T(0)e^{et}+e^{et}{\int }_0^t{e}^{et}{Q}_5(r)dr\ge 0,$$(16) (17) $$A(t)=T(0)e^{ft}+e^{ft}{\int }_0^t{e}^{ft}{Q}_6(r)dr\ge 0,$$(17) (18) $$R(t)=R(0)e^{dt}+e^{dt}{\int }_0^t{e}^{dt}{Q}_7(r)dr\ge 0,$$(18) where $Q_2(t)={\beta }_3{I}_1(t)+{\beta }_4{I}_2(t)+d,$ $Q_3(t)=-{\beta }_{1}S(t)-{\beta }_{3}E(t)+b,$ $Q_4(t)={\beta }_{2}S(t)+{\beta }_{4}E(t)-c,$ $Q_5(t)=k_1{I}_1(t)+k_3{I}_2(t),$ $Q_6(t)=k_2{I}_2(t)+{\alpha }_{2}T(t)$ and $Q_7(t)=\mu S(t).$ Consequently, we can state that $S(t)\ge 0,$ $I_1(t)\ge 0,$ $I_2(t)\ge 0,$ $A(t)\ge 0,$ $T(t)\ge 0,$ and $R(t)\ge 0$ for all $t\ge 0,$ and this completes the proof. □

2.3. Existence of equilibrium points and basic reproduction number

The equilibrium points are obtained by equating Equation (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) to zero. For simplicity of analysis, we exclude the R variable because this variable does not affect directly the spread of HIV in the model. We get equilibrium points, called free of infection and endemic points. Free of infection equilibrium point $(K^0)$ is (19) $$\begin{array}{cc}{K}^0=(S^0,E^0,I_1{}^0,I_2{}^0,T^0,A^0)\hfill & =\left(\frac{\Lambda }{a},\frac{\eta \Lambda }{da},0,0,0,0\right).\hfill \end{array}$$(19)

The sign of each component of the point reveals a disease-free equilibrium point. The fact that the sign is positive indicates that the system’s problem has an existing solution because it displays the population. From the free-disease point obtained, we can say that there is no infection transmission from the infected subpopulation to all the subpopulations in the model. Next, we get the endemic equilibrium point as follows $K^{\ast }=(S^{\ast },E^{\ast },I_1{}^{\ast },I_2{}^{\ast },T^{\ast },A^{\ast })$ (20) $$\begin{array}{c}{S}^{\ast }=\frac{\Lambda }{{\beta }_1{I}_1{}^{\ast }+{\beta }_2{I}_2{}^{\ast }+a},\\ E^{\ast }=\frac{\eta \Lambda }{\left({\beta }_1{I}_1{}^{\ast }+{\beta }_2{I}_2{}^{\ast }+a\right)\left({\beta }_3{I}_1{}^{\ast }+{\beta }_4{I}_2{}^{\ast }+d\right)},\\ T^{\ast }=\frac{k_1{I}_1{}^{\ast }+k_3{I}_2{}^{\ast }}{e},\\ A^{\ast }=\frac{(e{k}_2+{\alpha }_2{k}_3)I_2{}^{\ast }+{\alpha }_2{k}_1{I}_1{}^{\ast }}{fe},\\ {\left(I_2{}^{\ast }\right)}_{1,2}=\frac{-b_2\pm \sqrt{b_2{}^2-4b_1{b}_3}}{2b_1}.\end{array}$$(20) with $$\begin{array}{cc}{b}_1\hfill & =c{\beta }_2{\beta }_4,\hfill \\ b_2\hfill & =cd{\beta }_2+c\left({\beta }_1{\beta }_4+{\beta }_2{\beta }_3\right)I_1{}^{\ast }+ca{\beta }_4-{\beta }_2\Lambda {\beta }_4,\hfill \\ b_3\hfill & =c{\beta }_1{\beta }_3{I}_1{{}^{\ast }}^2+cd{\beta }_1{I}_1{}^{\ast }+ca{\beta }_3{I}_1{}^{\ast }+\mathit{acd}-{\beta }_2\Lambda {\beta }_3{I}_1{}^{\ast }-\mathit{acd}-{\beta }_4\eta \Lambda .\hfill \end{array}$$

Furthermore, $I_1{}^{\ast }$ can be obtained by solving the following equation: (21) $$e{\beta }_1\Lambda I_1{}^{\ast }\varphi +e{\beta }_3\eta \Lambda I_1{}^{\ast }+{\alpha }_1\theta \varphi \left(k_1{I}_1{}^{\ast }+k_3{I}_2{}^{\ast }\right)-\mathit{\text{eθϕb}}{I}_1{}^{\ast }=0,$$(21) with $$\begin{array}{cc}\theta \hfill & ={\beta }_1{I}_1{}^{\ast }+{\beta }_2{I}_2{}^{\ast }+a\hfill \\ \varphi \hfill & ={\beta }_3{I}_1{}^{\ast }+{\beta }_4{I}_2{}^{\ast }+d.\hfill \end{array}$$

Additionally, an endemic equilibrium point must show a positive sign. All points in the endemic point depend on the $I_1^{\ast }$ and $I_2^{\ast },$ thus we need to find a solution that works for both. Since we established the viability of the solutions in the preceding section, we can affirm that the endemic equilibrium exists.

In dynamical analysis, the basic reproduction number $(R_0)$ is a unique number determining whether infectious disease occurs or not (Murray, 2002). The basic reproduction number $(R_0)$ is calculated using the Next Generation Matrix (NGM) (Heffernan et al., 2005). From Equation (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) , we consider only the infected individuals, so we have $$x_i^{{}^{\prime }}=F_i-V_i,>i=1,2$$ where $F=({\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1,{\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2)$ and $V=(b{I}_1,c{I}_2).$ By taking partial derivative with respect to $I_1$ and $I_2$ at point $K_0,$ we get the Jacobian matrix F and V as follows $$F(K^0)=\left[\begin{array}{cc}{\beta }_1{S}^0+{\beta }_3{E}^0& 0\\ 0& {\beta }_2{S}^0+{\beta }_4{E}^0\end{array}\right],$$ and $$V(K^0)=\left[\begin{array}{cc}b& 0\\ 0& c\end{array}\right].$$

Thus, by taking invers $V(K^0)$ and multiply with $F(K^0),$ we end up the NGM as follows (22) $$\begin{array}{cc}K=F(K_0)V^{-1}(K_0),\hfill & \hfill \\ =\left[\begin{array}{cc}\frac{{\beta }_1{S}^0+{\beta }_3{E}^0}{b}& 0\\ 0& \frac{{\beta }_2{S}^0+{\beta }_4{E}^0}{c}\end{array}\right].\hfill & \hfill \end{array}$$(22)

Thus, by calculating $|K-\lambda I|=0,$ we get eigen values as follows (23) $${\lambda }_1=\frac{{\beta }_1{S}^0+{\beta }_3{E}^0}{b},$$(23) and (24) $${\lambda }_2=\frac{{\beta }_2{S}^0+{\beta }_4{E}^0}{c}.$$(24)

The basic reproduction number $R_0$ is equal to the maximum modulus of matrix eigenvalues K, i.e. $R_0=$ max $({\lambda }_1,{\lambda }_2).$ Substituting $S^0=\frac{\Lambda }{a}$ and $E^0=\frac{\eta \Lambda }{da}$ from free-disease equilibrium point (19)(19) $$\begin{array}{cc}{K}^0=(S^0,E^0,I_1{}^0,I_2{}^0,T^0,A^0)\hfill & =\left(\frac{\Lambda }{a},\frac{\eta \Lambda }{da},0,0,0,0\right).\hfill \end{array}$$(19) into Equations (23)(23) $${\lambda }_1=\frac{{\beta }_1{S}^0+{\beta }_3{E}^0}{b},$$(23) and (24)(24) $${\lambda }_2=\frac{{\beta }_2{S}^0+{\beta }_4{E}^0}{c}.$$(24) , as well as using the information described in subsection 2.1, ${\beta }_1,{\beta }_2>{\beta }_3,{\beta }_4$ and ${\beta }_1<{\beta }_2$ as a condition for the transmission of HIV, then we get (25) $$R_0=\frac{{\beta }_{2}d\Lambda +{\beta }_4\eta \Lambda }{\mathit{acd}}.$$(25)

The particular cases of multi-interaction between subpopulations in the HIV/AIDS model and their analysis have been written in Habibah et al. (2021a), Habibah, Trisilowati, Pradana, et al. (2021), and Habibah, Owati, and Muzaqi (2021).

2.4. Sensitivity analysis

A sensitivity analysis was done to examine the effect of each parameter on infection spread depending on the $R_0,$ which can be helpful for choosing the best control methods to stop an epidemic from spreading (Rois et al., 2021). The following formula can be used to determine the sensitivity index value: $$C_p^A=\frac{\partial A}{\partial p}\frac{p}{A},$$ where A is a differentiated function of a parameter p.

It is crucial to identify the virus’s primary route of transmission in order to decrease the number of HIV-related illnesses and fatalities. This sensitivity index can be used to determine how important certain factors are for the spread of the disease. As a result, this method was used to identify the parameters that need to be regulated and those that have a substantial impact on the $R_0.$ The basic reproduction number $R_0$ has the following parameters: ${\beta }_2,$ ${\beta }_4,$ d, $\Lambda ,$ $\eta ,$ $\mu ,$ $k_2,$ and $k_3.$ Using the technique described by Erfanian and Noori Skandari (2011) and Rois et al. (2021), the parameter sensitivity analysis can be calculated as follows: (26) $$\begin{array}{cc}{C}_{{\beta }_2}^{R_0}\hfill & =\frac{\partial R_0}{\partial {\beta }_2}\times \frac{{\beta }_2}{R_0}=\frac{{\beta }_2\Lambda d}{{\beta }_{2}d\Lambda +{\beta }_4\eta \Lambda },\hfill \\ C_{{\beta }_4}^{R_0}\hfill & =\frac{\partial R_0}{\partial {\beta }_4}\times \frac{{\beta }_2}{R_0}=\frac{{\beta }_4\Lambda \eta }{{\beta }_{2}d\Lambda +{\beta }_4\eta \Lambda },\hfill \\ C_{\Lambda }^{R_0}\hfill & =\frac{\partial R_0}{\partial \Lambda }\times \frac{\Lambda }{R_0}=\frac{\Lambda ({\beta }_{2}d+{\beta }_4\eta )}{{\beta }_{2}d\Lambda +{\beta }_4\eta \Lambda },\hfill \\ C_{\mu }^{R_0}\hfill & =\frac{\partial R_0}{\partial \mu }\times \frac{\mu }{R_0}=-\frac{\mu }{\mu +d+\eta },\hfill \\ C_{k_2}^{R_0}\hfill & =\frac{\partial R_0}{\partial k_2}\times \frac{k_2}{R_0}=-\frac{k_2}{k_2+k_3+d},\hfill \\ C_{k_3}^{R_0}\hfill & =\frac{\partial R_0}{\partial k_3}\times \frac{k_3}{R_0}=-\frac{k_3}{k_2+k_3+d},\hfill \\ C_d^{R_0}\hfill & =\frac{\partial R_0}{\partial d}\times \frac{d}{R_0}=DD,\hfill \\ C_{\eta }^{R_0}\hfill & =\frac{\partial R_0}{\partial \eta }\times \frac{d}{R_0}=EE,\hfill \end{array}$$(26) where (27) $$\begin{array}{cc}DD=\frac{d^2(\mu +\eta +d)(k_2+k_3+d)}{{\beta }_{2}d\Lambda +{\beta }_4\eta \Lambda }(\frac{{\beta }_2\Lambda }{(\mu +\eta +d)d(k_2+k_3+d)}-\frac{{\beta }_{2}d\Lambda +{\beta }_4\eta \Lambda }{{(\mu +\eta +d)}^{2}d(k_2+k_3+d)}\hfill & \hfill \\ -\frac{{\beta }_{2}d\Lambda +{\beta }_4\eta \Lambda }{(\mu +\eta +d)d^2(k_2+k_3+d)}-\frac{{\beta }_{2}d\Lambda +{\beta }_4\eta l\Lambda }{(\mu +\eta +d)d{(k_2+k_3+d)}^2}),\hfill & \hfill \end{array}$$(27) and (28) $$\begin{array}{cc}EE\hfill & =\frac{d\eta (\mu +\eta +d)(k_2+k_3+d)}{{\beta }_{2}d\Lambda +{\beta }_4\eta \Lambda }\left(\frac{{\beta }_4\Lambda }{(\mu +\eta +d)d(k_2+k_3+d)}-\frac{{\beta }_{2}d\Lambda +{\beta }_4\eta \Lambda }{{(\mu +\eta +d)}^{2}d(k_2+k_3+d)}\right).\hfill \end{array}$$(28)

Considering the values of the parameters in Table 1, we obtained $R_0=0.5077$ and the sensitivity index values of all parameters in the reproduction number $R_0$ can be seen in Table 2. A visual representation of the data in Table 2 can be seen in Figure 2.

Figure 2. Sensitivity index of parameters in the reproduction number $R_0$ in the form of a diagram.

Table 2. Sensitivity index of parameters in the reproduction number $R_0.$

Parameter	Sensitivity index value
β2	0.7729
${\beta }_4$	0.2270
d	−0.8085
$\Lambda$	1
$\eta$	0.0593
$\mu$	−0.5034
$k_2$	−0.1031
$k_3$	−0.6443

Based on Table 2 or Figure 2, four parameters have positive sensitivity index values (${\beta }_2,$ ${\beta }_4,$ $\Lambda ,$ and $\eta$), while the others are negative. A positive mark of sensitivity index value means that if the parameter value is increased (decreased), then the $R_0$ will increase (decrease). This causes the spread of HIV/AIDS to increase (decrease). As with positive values, a negative mark of sensitivity index value (${\mu }_1,$ d, $k_2,$ and $k_3$) means that if the parameter value is increased (decreased), then the $R_0$ will decrease (increase). As a result, the spread of HIV/AIDS decreases (increases). Positive and negative sensitivity index values reflect opposite parameter behavior. Now, we focus on the positive sensitivity index values. $\Lambda$ has a higher value of sensitivity index, but it does not contribute to the transmission of infection. Therefore, we avoid this parameter. On the other hand, ${\beta }_2$ is the parameter that has contributed to the increase in infection the most. The discussion of the parameters that affect how this disease spreads were discovered using parameter sensitivity analysis can be seen in Sangsawang, Humphries, Khan, and Pongsumpun (2023).

2.5. Stability analysis

The free-of-infection equilibrium point can be analyzed for its stability by calculating those of the Jacobian matrix’s eigenvalues of the system (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) . If all the eigenvalues have negative values, then the equilibrium point $K^0$ is asymptotically stable. To show this, first, the system (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) should be linearized to have the Jacobian matrix at $K^0$ as follows (29) $$J\left(K^0\right)=\left(\begin{array}{cccccc}-a& 0& -{\beta }_1\frac{\Lambda }{a}& -{\beta }_2\frac{\Lambda }{a}& 0& 0\\ \eta & -d& -{\beta }_3\frac{\eta \Lambda }{da}& -{\beta }_4\frac{\eta \Lambda }{da}& 0& 0\\ 0& 0& {\beta }_1\frac{\Lambda }{a}+{\beta }_3\frac{\eta \Lambda }{da}-b& 0& {\alpha }_1& 0\\ 0& 0& 0& {\beta }_2\frac{\Lambda }{a}+{\beta }_4\frac{\eta \Lambda }{da}-c& 0& 0\\ 0& 0& k_1& k_3& -e& 0\\ 0& 0& 0& k_2& {\alpha }_2& -f\end{array}\right).$$(29)

By solving $\left|J(K^0)-rI\right|=0,$ we get the characteristic equation of (29)(29) $$J\left(K^0\right)=\left(\begin{array}{cccccc}-a& 0& -{\beta }_1\frac{\Lambda }{a}& -{\beta }_2\frac{\Lambda }{a}& 0& 0\\ \eta & -d& -{\beta }_3\frac{\eta \Lambda }{da}& -{\beta }_4\frac{\eta \Lambda }{da}& 0& 0\\ 0& 0& {\beta }_1\frac{\Lambda }{a}+{\beta }_3\frac{\eta \Lambda }{da}-b& 0& {\alpha }_1& 0\\ 0& 0& 0& {\beta }_2\frac{\Lambda }{a}+{\beta }_4\frac{\eta \Lambda }{da}-c& 0& 0\\ 0& 0& k_1& k_3& -e& 0\\ 0& 0& 0& k_2& {\alpha }_2& -f\end{array}\right).$$(29) . Furthermore, we obtain eigenvalues as follows $r_1=-d,$ $r_2=-f,$ $r_3=-a,$ $r_4={\beta }_2\frac{\Lambda }{a}+{\beta }_4\frac{\eta \Lambda }{da}-c,$ where value of $r_4$ can be written as $r_4={\beta }_2\frac{\Lambda }{a}+{\beta }_4\frac{\eta \Lambda }{da}-c=c(R_0-1).$ Value of $r_4$ is negative when $R_0<1.$ The eigenvalue $r_5$ and $r_6$ can be obtained by calculating the following equation (30) $$\left|\begin{array}{cc}{\beta }_1\frac{\Lambda }{a}+{\beta }_3\frac{\eta \Lambda }{da}-b-r& {\alpha }_1\\ k_1& -e-r\end{array}\right|=0,$$(30) we get characteristic equation (31) $$r^2+(e-g)r-ge-k_1{\alpha }_1=0,$$(31) where $g={\beta }_1\frac{\Lambda }{a}+{\beta }_3\frac{\eta \Lambda }{da}-b.$ The Equation (31)(31) $$r^2+(e-g)r-ge-k_1{\alpha }_1=0,$$(31) has negative roots ($r_5<0$ and $r_6<0$) when satisfies the discriminant of Equation (31)(31) $$r^2+(e-g)r-ge-k_1{\alpha }_1=0,$$(31) is positive, $r_5\cdot r_6=-ge-{\alpha }_1{k}_1>0$ and $r_5+r_6=-(e-g)<0.$ Finally, we obtain all of the eigenvalues that are negative which can conclude that the disease-free equilibrium point $(K^0)$ is locally asymptotically stable when $R_0<1.$

The stability of the endemic equilibrium point $(K^{\ast })$ is obtained using the Routh-Hurwitz criteria (Martcheva, 2015; Tunc, 2007, 2002; 2009). Firstly, we find the eigenvalues of the Jacobian matrix of the system (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) at $K^{\ast }.$ The following is the Jacobian matrix at $K^{\ast }$ (32) $$J\left(K^{\ast }\right)=\left(\begin{array}{cccccc}-H_1& 0& -{\beta }_1{S}^{*}& -{\beta }_2{S}^{*}& 0& 0\\ \eta & -H_2& -{\beta }_3{E}^{*}& -{\beta }_4{E}^{*}& 0& 0\\ {\beta }_1{I}_1^{*}& {\beta }_3{I}_1^{*}& \Psi & 0& {\alpha }_1& 0\\ {\beta }_2{I}_2^{*}& {\beta }_4{I}_2^{*}& 0& \Omega & 0& 0\\ 0& 0& k_1& k_3& -e& 0\\ 0& 0& 0& k_2& {\alpha }_2& -f\end{array}\right),$$(32) with (33) $$\begin{array}{cc}{H}_1\hfill & ={\beta }_1{I}_1{}^{\ast }+{\beta }_2{I}_2{}^{\ast }+a,\hfill \end{array}$$(33) (34) $$\begin{array}{cc}\Omega \hfill & ={\beta }_2{S}^{\ast }+{\beta }_4{E}^{\ast }-c,\hfill \end{array}$$(34) (35) $$\begin{array}{cc}\Psi \hfill & ={\beta }_1{S}^{\ast }+{\beta }_3{E}^{\ast }-b,\hfill \end{array}$$(35) (36) $$\begin{array}{cc}{H}_2\hfill & ={\beta }_3{I}_1{}^{\ast }+{\beta }_4{I}_2{}^{\ast }+d.\hfill \end{array}$$(36)

The equation of characteristic of matrix $J(K^{\ast })$ is yielded by calculating $\left|J(K^{\ast })-rI\right|=0$ (37) $$\left|\begin{array}{cccccc}-H_1-r& 0& -{\beta }_1{S}^{\ast }& -{\beta }_2{S}^{\ast }& 0& 0\\ \eta & -H_2-r& -{\beta }_3{E}^{\ast }& -{\beta }_4{E}^{\ast }& 0& 0\\ {\beta }_1{I}_1{}^{\ast }& {\beta }_3{I}_1{}^{\ast }& \Psi -r& 0& {\alpha }_1& 0\\ {\beta }_2{I}_2{}^{\ast }& {\beta }_4{I}_2{}^{\ast }& 0& \Omega -r& 0& 0\\ 0& 0& k_1& k_3& -e-r& 0\\ 0& 0& 0& k_2& {\alpha }_2& -f-r\end{array}\right|=0,$$(37) hence we obtain the following eigenvalues $r_1=-f,$ and $r_2,$ $r_3,$ $r_4,$ $r_5,$ $r_6$ that meets (38) $$r^5+b_1{r}^4+b_2{r}^3+b_3{r}^2+b_{4}r+b_5=0,$$(38) with $$\begin{array}{c}{b}_1=-\left(\Psi +\Omega -H_2-H_1-e\right),\\ b_2=\Psi \Omega +H_2\Psi +H_2\Omega +{\beta }_4{E}^{\ast }{\beta }_4{I}_2{}^{\ast }+{\beta }_3{E}^{\ast }{\beta }_3{I}_1{}^{\ast }\\ -H_1\Psi -H_1\Omega +H_1{H}_2+{\beta }_1{S}^{\ast }{\beta }_1{I}_1{}^{\ast }\\ +{\beta }_2{S}^{\ast }{\beta }_2{I}_2{}^{\ast }-\Psi e-\Omega e+H_{2}e+e{H}_1-{\alpha }_1{k}_1,\\ b_3=-H_2\Psi \Omega -{\beta }_4{E}^{\ast }{\beta }_4{I}_2{}^{\ast }\Psi -{\beta }_3{E}^{\ast }{\beta }_3{I}_1{}^{\ast }\Omega +H_1\Psi \Omega \\ +H_1{H}_2\Psi +H_1{H}_2\Omega +H_1{\beta }_4{E}^{\ast }{\beta }_4{I}_2{}^{\ast }+H_1{\beta }_3{E}^{\ast }{\beta }_3{I}_1{}^{\ast }\\ -{\beta }_1{S}^{\ast }\left({\beta }_1{I}_1{}^{\ast }\Omega -{\beta }_1{I}_1{}^{\ast }{H}_2-\eta {\beta }_3{I}_1{}^{\ast }\right)+{\beta }_2{S}^{\ast }{\beta }_2{I}_2{}^{\ast }{H}_2\\ -{\beta }_2{S}^{\ast }{\beta }_2{I}_2{}^{\ast }\Psi +{\beta }_2{S}^{\ast }\eta {\beta }_4{I}_2{}^{\ast }+e\Psi \Omega +e{H}_2\Psi \\ +e{H}_2\Omega +e{\beta }_4{E}^{\ast }{\beta }_4{I}_2{}^{\ast }+e{\beta }_3{E}^{\ast }{\beta }_3{I}_1{}^{\ast }-e{H}_1\Psi \\ -e{H}_1\Omega +e{H}_1{H}_2+e{\beta }_1{S}^{\ast }{\beta }_1{I}_1{}^{\ast }+e{\beta }_2{S}^{\ast }{\beta }_2{I}_2{}^{\ast }\\ -{\alpha }_1{H}_1{k}_1-{\alpha }_1{H}_2{k}_1+{\alpha }_1\Omega k_1,\\ b_4=-H_1{H}_2\Psi \Omega -H_1{\beta }_4{E}^{\ast }{\beta }_4{I}_2{}^{\ast }\Psi -H_1{\beta }_3{E}^{\ast }{\beta }_3{I}_1{}^{\ast }\Omega \\ -{\beta }_1{}^2{I}_1{}^{\ast }{H}_2\Omega S^{\ast }-\eta {\beta }_3{I}_1{}^{\ast }\Omega S^{\ast }+{\beta }_4{E}^{\ast }{\beta }_1{}^2{I}_1{}^{\ast }{\beta }_4{I}_2{}^{\ast }{S}^{\ast }\\ -{\beta }_1{S}^{\ast }{\beta }_4{E}^{\ast }{\beta }_3{I}_1{}^{\ast }{\beta }_2{I}_2{}^{\ast }+{\beta }_3{E}^{\ast }{\beta }_3{I}_1{}^{\ast }{\beta }_2{}^2{I}_2{}^{\ast }{S}^{\ast }\\ -{\beta }_2{}^2{I}_2{}^{\ast }{H}_2\Psi S^{\ast }-{\beta }_2{S}^{\ast }{\beta }_3{E}^{\ast }{\beta }_1{I}_1{}^{\ast }{\beta }_4{I}_2{}^{\ast }\\ -{\beta }_2{S}^{\ast }\eta {\beta }_4{I}_2{}^{\ast }\Psi -e{H}_2\Psi \Omega -e{\beta }_4{E}^{\ast }{\beta }_4{I}_2{}^{\ast }\Psi \\ -e{\beta }_3{E}^{\ast }{\beta }_3{I}_1{}^{\ast }\Omega +e{H}_1\Psi \Omega +e{H}_1{H}_2\Psi +e{H}_1{H}_2\Omega \\ +e{H}_1{\beta }_4{E}^{\ast }{\beta }_4{I}_2{}^{\ast }+e{H}_1{\beta }_3{E}^{\ast }{\beta }_3{I}_1{}^{\ast }-{\beta }_1{}^2{I}_1{}^{\ast }\Omega e{S}^{\ast }\\ +{\beta }_1{}^2{I}_1{}^{\ast }{H}_{2}e{S}^{\ast }+e{\beta }_1{S}^{\ast }\eta {\beta }_3{I}_1{}^{\ast }+e{\beta }_2{S}^{\ast }{\beta }_2{I}_2{}^{\ast }{H}_2\\ -e{\beta }_2{S}^{\ast }{\beta }_2{I}_2{}^{\ast }\Psi +e{\beta }_2{S}^{\ast }\eta {\beta }_4{I}_2{}^{\ast }-{\alpha }_1{H}_1{H}_2{k}_1\\ +{\alpha }_1{H}_1\Omega k_1+{\alpha }_1{H}_2\Omega k_1-{\alpha }_1{\beta }_4{E}^{\ast }{\beta }_4{I}_2{}^{\ast }{k}_1\\ -{\alpha }_1{\beta }_2{S}^{\ast }{\beta }_2{I}_2{}^{\ast }{k}_1+{\alpha }_1{k}_3{\beta }_1{S}^{\ast }{\beta }_2{I}_2{}^{\ast }+{\alpha }_1{k}_3{\beta }_3{E}^{\ast }{\beta }_4{I}_2{}^{\ast },\\ b_5=-e{H}_1{H}_2\Psi \Omega -e{H}_1{\beta }_4{E}^{\ast }{\beta }_4{I}_2{}^{\ast }\Psi -e{H}_1{\beta }_3{E}^{\ast }{\beta }_3{I}_1{}^{\ast }\Omega \\ -{\beta }_1{}^2{I}_1{}^{\ast }{H}_2\Omega e{S}^{\ast }-e{\beta }_1{S}^{\ast }\eta {\beta }_3{I}_1{}^{\ast }\Omega \\ +e{S}^{\ast }{\beta }_4{E}^{\ast }{\beta }_1{}^2{I}_1{}^{\ast }{\beta }_4{I}_2{}^{\ast }-e{\beta }_1{S}^{\ast }{\beta }_4{E}^{\ast }{\beta }_3{I}_1{}^{\ast }{\beta }_2{I}_2{}^{\ast }\\ +E^{\ast }{\beta }_3{}^3{I}_1{}^{\ast }{\beta }_2{}^2{I}_2{}^{\ast }e{S}^{\ast }-e{S}^{\ast }{\beta }_2{}^2{I}_2{}^{\ast }{H}_2\Psi \\ -e{\beta }_2{S}^{\ast }{\beta }_3{E}^{\ast }{\beta }_1{I}_1{}^{\ast }{\beta }_4{I}_2{}^{\ast }-e{\beta }_2{S}^{\ast }\eta {\beta }_4{I}_2{}^{\ast }\Psi \\ -{\alpha }_1{H}_1{H}_2\Omega k_1+{\alpha }_1\eta {\beta }_2{S}^{\ast }{\beta }_4{I}_2{}^{\ast }{k}_1\\ +{\alpha }_1{\beta }_2{S}^{\ast }{\beta }_2{I}_2{}^{\ast }{H}_2{k}_1+{\alpha }_1{\beta }_4{E}^{\ast }{\beta }_4{I}_2{}^{\ast }{H}_1{k}_1\\ -{\alpha }_1\eta {\beta }_1{S}^{\ast }{\beta }_4{I}_2{}^{\ast }{k}_3-{\alpha }_1{\beta }_1{S}^{\ast }{\beta }_2{I}_2{}^{\ast }{H}_2{k}_3\\ -{\alpha }_1{k}_3{\beta }_3{E}^{\ast }{\beta }_4{I}_2{}^{\ast }{H}_1.\end{array}$$

Since the characteristic equation has a complicated formula, it is rather difficult to determine the eigenvalues. Hence, we apply the Routh-Hurwitz criteria (Martcheva, 2015; Tunc, 2009, 2007, 2002). The characteristics equation has negative roots if and only if $b_1>0$ and

1. $b_1{b}_2-b_3>0,$

2. $b_1{b}_2{b}_3-b_3{}^2-b_1{}^2{b}_4>0,$

3. $b_1{b}_2{b}_3{b}_4+2b_1{b}_4{b}_5+b_2{b}_3{b}_5-b_1{}^2{b}_4{}^2-b_1{b}_2{}^2{b}_5-b_3{}^2{b}_4-b_5{}^2>0,$

4. $b_5>0.$

In the next subsection, we present a bifurcation analysis of the HIV/AIDS model.

2.6. Backward bifurcation analysis

To analyze the backward bifurcation, we apply the Center-Manifold Theorem which is introduced in Castillo-Chavez and Song (2004). In the Center-Manifold Theorem, we should calculate two important values, P and Q. These values as a key to determining the type of bifurcation when $R_0=1.$ A forward bifurcation occurs when $P<0$ and $Q>0.$ On the other hand, a backward bifurcation occurs when $P>0$ and $Q>0.$

Theorem 2.3

(Backward Bifurcation). If $R_0<1$ and (39) $$\begin{array}{cc}P\hfill & =\frac{2}{\left(k_3{\alpha }_1{a}^{2}d\right)\left(k_3{\alpha }_{1}ad\right)\left(d\Lambda \right)}\hfill \\ (\left(-k_3{\alpha }_{1}acd+k_3{\alpha }_1{\beta }_4\eta \Lambda \right)(k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd})-{\beta }_1\Lambda k_3{\alpha }_{1}d)\\ \left(-k_1{\alpha }_{1}ad-e{\beta }_{1}d\Lambda -e{\beta }_3\eta \Lambda +\mathit{eabd}\right)\left(\mathit{acd}-{\beta }_4\eta \Lambda \right)\\ +\frac{2{\beta }_4}{\left(k_3{\alpha }_1{a}^2{d}^2\right)\left(k_3{\alpha }_{1}ad\right)}\\ (\eta d\left(\left(-k_3{\alpha }_{1}acd+k_3{\alpha }_1{\beta }_4\eta \Lambda \right)(k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd})-{\beta }_1\Lambda k_3{\alpha }_{1}d\right)\\ -{\beta }_3\eta \Lambda k_3{\alpha }_{1}ad+\left({\beta }_4\eta \Lambda \right)(k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd}))\\ \left(-k_1{\alpha }_{1}ad-e{\beta }_{1}d\Lambda -e{\beta }_3\eta \Lambda +\mathit{eabd}\right)>0,\end{array}$$(39) then $SE{I}_1{I}_{2}TAR$ model exhibits a backward bifurcation at $R_0=1$. Oppositely, if $P<0$ then $SE{I}_1{I}_{2}TAR$ model exhibits a forward bifurcation at $R_0=1.$

Proof.

First, let the system (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) excluding the last equation, be redefined as follows (40) $$\begin{array}{cc}{f}_1\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \end{array}$$(40) (41) $$\begin{array}{cc}{f}_2\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \end{array}$$(41) (42) $$\begin{array}{cc}{f}_3\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \end{array}$$(42) (43) $$\begin{array}{cc}{f}_4\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \end{array}$$(43) (44) $$\begin{array}{cc}{f}_5\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \end{array}$$(44) (45) $$\begin{array}{cc}{f}_6\hfill & =k_2{I}_2+{\alpha }_{2}T-fA.\hfill \end{array}$$(45)

Let ${\beta }_2$ be the bifurcation parameter. We choose ${\beta }_2$ as a bifurcation parameter since it shows the interaction between uneducated and HIV-positive without consuming ARV individuals that give the most contribution to the spread of HIV. It is consistent with the real phenomena in the spread of HIV/AIDS. The value of ${\beta }_2$ is obtained from the basic reproduction number and evaluated when $R_0=1,$ such that we have (46) $${\beta }_2^{\ast }=\frac{\mathit{acd}-{\beta }_4\eta \Lambda }{d\Lambda }.$$(46)

By substituting ${\beta }_2^{\ast }$ and $K^0$ into the Jacobian matrix (29)(29) $$J\left(K^0\right)=\left(\begin{array}{cccccc}-a& 0& -{\beta }_1\frac{\Lambda }{a}& -{\beta }_2\frac{\Lambda }{a}& 0& 0\\ \eta & -d& -{\beta }_3\frac{\eta \Lambda }{da}& -{\beta }_4\frac{\eta \Lambda }{da}& 0& 0\\ 0& 0& {\beta }_1\frac{\Lambda }{a}+{\beta }_3\frac{\eta \Lambda }{da}-b& 0& {\alpha }_1& 0\\ 0& 0& 0& {\beta }_2\frac{\Lambda }{a}+{\beta }_4\frac{\eta \Lambda }{da}-c& 0& 0\\ 0& 0& k_1& k_3& -e& 0\\ 0& 0& 0& k_2& {\alpha }_2& -f\end{array}\right).$$(29) , we have (47) $$J\left(K^0,{\beta }_2^{*}\right)=\left(\begin{array}{cccccc}-a& 0& -{\beta }_1\frac{\Lambda }{a}& \frac{-\mathit{acd}+{\beta }_4\eta \Lambda }{ad}& 0& 0\\ \eta & -d& -{\beta }_3\frac{\eta \Lambda }{da}& -{\beta }_4\frac{\eta \Lambda }{da}& 0& 0\\ 0& 0& {\beta }_1\frac{\Lambda }{a}+{\beta }_3\frac{\eta \Lambda }{da}-b& 0& {\alpha }_1& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& k_1& k_3& -e& 0\\ 0& 0& 0& k_2& {\alpha }_2& -f\end{array}\right).$$(47)

Then, we calculate the characteristic equation of $\left|J(K^0,{\beta }_2^{*})-rI\right|=0$ to yield $r_1=a,$ $r_2=d,$ $r_3=0,$ $r_4=f,$ and $r_{5,6}$ that satisfies (48) $$\left|\begin{array}{cc}{\beta }_1\frac{\Lambda }{a}+{\beta }_3\frac{\eta \Lambda }{da}-b-r& {\alpha }_1\\ k_1& -e-r\end{array}\right|=0,$$(48) which has the same form as the eigenvalues in the free-of-infection stability. From the above calculation, we obtain zero eigenvalues ($r_3=0$) and the other eigenvalues are negative real parts. It means this meets the condition of the Center-Manifold Theorem.

Furthermore, we need to calculate the left and right eigenvectors related to the zero eigenvalues. Then, the right eigenvector is ${\left[w_1{w}_2{w}_3{w}_4{w}_5{w}_6\right]}^T,$ where $$\begin{array}{c}{w}_1=\frac{\left(-k_3{\alpha }_{1}acd+k_3{\alpha }_1{\beta }_4\eta \Lambda \right)\left(k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd}\right)-{\beta }_1\Lambda k_3{\alpha }_{1}d}{k_3{\alpha }_1{a}^{2}d},\\ w_2=\eta \left(\frac{\left(-k_3{\alpha }_{1}acd+k_3{\alpha }_1{\beta }_4\eta \Lambda \right)\left(k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd}\right)-{\beta }_1\Lambda k_3{\alpha }_{1}d}{k_3{\alpha }_1{a}^{2}d}\right)\\ -{\beta }_3\frac{\eta \Lambda }{ad}+\left(\frac{{\beta }_4\eta \Lambda }{ad}\right)\left(\frac{k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd}}{k_3{\alpha }_{1}ad}\right),\\ w_3=w_3,\\ w_4=-\frac{k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd}}{k_3{\alpha }_{1}ad},\\ w_5=-\frac{{\beta }_{1}d\Lambda +{\beta }_3\eta \Lambda -\mathit{abd}}{{\alpha }_{1}ad},\\ w_6=\frac{-k_1{k}_2{\alpha }_{1}ad-k_{2}e{\beta }_{1}d\Lambda -k_{2}e{\beta }_3\eta \Lambda +k_{2}e\mathit{abd}-k_3{\alpha }_2{\beta }_{1}d\Lambda -k_3{\alpha }_2{\beta }_3\eta \Lambda +k_3{\alpha }_{2}abd}{f{k}_3{\alpha }_{1}ad}.\end{array}$$

On the other hand, the left eigenvector is $\left[v_1{v}_2{v}_3{v}_4{v}_5{v}_6\right],$ where $v_1=v_2=v_3=v_5=v_6=0,v_4=v_4.$ Since $v_1,v_2,v_3,v_5,v_6=0,$ we only need to calculate the second partial derivative of $f_4.$ Hence, we have (49) $$\begin{array}{cc}P\hfill & =v_4\left(w_1{w}_4\frac{{\partial }^2{f}_4}{\partial S\partial I_2}\left(K^0,{\beta }_2{}^{\ast }\right)+w_2{w}_4\frac{{\partial }^2{f}_4}{\partial E\partial I_2}\left(K^0,{\beta }_2{}^{\ast }\right)+w_4{w}_1\frac{{\partial }^2{f}_4}{\partial I_2\partial S}\left(K^0,{\beta }_2{}^{\ast }\right)\right)\hfill \\ +v_4\left(w_4{w}_2\frac{{\partial }^2{f}_4}{\partial I_2\partial E}\left(K^0,{\beta }_2{}^{\ast }\right)\right),\\ =\frac{2}{\left(k_3{\alpha }_1{a}^{2}d\right)\left(k_3{\alpha }_{1}ad\right)\left(d\Lambda \right)}\\ \left(\left(-k_3{\alpha }_{1}acd+k_3{\alpha }_1{\beta }_4\eta \Lambda \right)\left(k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd}\right)-{\beta }_1\Lambda k_3{\alpha }_{1}d\right)\\ \left(-k_1{\alpha }_{1}ad-e{\beta }_{1}d\Lambda -e{\beta }_3\eta \Lambda +\mathit{eabd}\right)\left(\mathit{acd}-{\beta }_4\eta \Lambda \right)\\ +\frac{2{\beta }_4}{\left(k_3{\alpha }_1{a}^2{d}^2\right)\left(k_3{\alpha }_{1}ad\right)}\\ (\eta d\left(\left(-k_3{\alpha }_{1}acd+k_3{\alpha }_1{\beta }_4\eta \Lambda \right)\left(k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd}\right)-{\beta }_1\Lambda k_3{\alpha }_{1}d\right)\\ -{\beta }_3\eta \Lambda k_3{\alpha }_{1}ad+\left({\beta }_4\eta \Lambda \right)\left(k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd}\right))\\ \left(-k_1{\alpha }_{1}ad-e{\beta }_{1}d\Lambda -e{\beta }_3\eta \Lambda +\mathit{eabd}\right),\end{array}$$(49) (50) $$\begin{array}{c}Q=v_4\left(w_1\frac{{\partial }^2{f}_4}{\partial S{\partial }_{{\beta }_2{}^{\ast }}}\left(K^0,{\beta }_2{}^{\ast }\right)+w_2\frac{{\partial }^2{f}_4}{\partial E{\partial }_{{\beta }_2{}^{\ast }}}\left(K^0,{\beta }_2{}^{\ast }\right)+w_4\frac{{\partial }^2{f}_4}{\partial I_2{\partial }_{{\beta }_2{}^{\ast }}}\left(K^0,{\beta }_2{}^{\ast }\right)\right),\\ =\left(-\frac{k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd}}{k_3{\alpha }_{1}ad}\right)\left(\frac{\Lambda }{a}\right).\end{array}$$(50)

It is difficult to determine a backward bifurcation occurs when $P>0$ and $Q>0$ directly from the Equations (49)(49) $$\begin{array}{cc}P\hfill & =v_4\left(w_1{w}_4\frac{{\partial }^2{f}_4}{\partial S\partial I_2}\left(K^0,{\beta }_2{}^{\ast }\right)+w_2{w}_4\frac{{\partial }^2{f}_4}{\partial E\partial I_2}\left(K^0,{\beta }_2{}^{\ast }\right)+w_4{w}_1\frac{{\partial }^2{f}_4}{\partial I_2\partial S}\left(K^0,{\beta }_2{}^{\ast }\right)\right)\hfill \\ +v_4\left(w_4{w}_2\frac{{\partial }^2{f}_4}{\partial I_2\partial E}\left(K^0,{\beta }_2{}^{\ast }\right)\right),\\ =\frac{2}{\left(k_3{\alpha }_1{a}^{2}d\right)\left(k_3{\alpha }_{1}ad\right)\left(d\Lambda \right)}\\ \left(\left(-k_3{\alpha }_{1}acd+k_3{\alpha }_1{\beta }_4\eta \Lambda \right)\left(k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd}\right)-{\beta }_1\Lambda k_3{\alpha }_{1}d\right)\\ \left(-k_1{\alpha }_{1}ad-e{\beta }_{1}d\Lambda -e{\beta }_3\eta \Lambda +\mathit{eabd}\right)\left(\mathit{acd}-{\beta }_4\eta \Lambda \right)\\ +\frac{2{\beta }_4}{\left(k_3{\alpha }_1{a}^2{d}^2\right)\left(k_3{\alpha }_{1}ad\right)}\\ (\eta d\left(\left(-k_3{\alpha }_{1}acd+k_3{\alpha }_1{\beta }_4\eta \Lambda \right)\left(k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd}\right)-{\beta }_1\Lambda k_3{\alpha }_{1}d\right)\\ -{\beta }_3\eta \Lambda k_3{\alpha }_{1}ad+\left({\beta }_4\eta \Lambda \right)\left(k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd}\right))\\ \left(-k_1{\alpha }_{1}ad-e{\beta }_{1}d\Lambda -e{\beta }_3\eta \Lambda +\mathit{eabd}\right),\end{array}$$(49) and (50)(50) $$\begin{array}{c}Q=v_4\left(w_1\frac{{\partial }^2{f}_4}{\partial S{\partial }_{{\beta }_2{}^{\ast }}}\left(K^0,{\beta }_2{}^{\ast }\right)+w_2\frac{{\partial }^2{f}_4}{\partial E{\partial }_{{\beta }_2{}^{\ast }}}\left(K^0,{\beta }_2{}^{\ast }\right)+w_4\frac{{\partial }^2{f}_4}{\partial I_2{\partial }_{{\beta }_2{}^{\ast }}}\left(K^0,{\beta }_2{}^{\ast }\right)\right),\\ =\left(-\frac{k_1{\alpha }_{1}ad+e{\beta }_{1}d\Lambda +e{\beta }_3\eta \Lambda -\mathit{eabd}}{k_3{\alpha }_{1}ad}\right)\left(\frac{\Lambda }{a}\right).\end{array}$$(50) , so we will show them in the numerical simulation by using MAPLE software.

3. Numerical simulation

In this part, we give some simulations to illustrate the analytical results. Numerical simulations are carried out using the Runge-Kutta 4th order since this method is suitable for solving the system of ordinary differential equations. Furthermore, we show the dynamics of the HIV/AIDS mathematical formula with the parameter values as in Table 1. The local stability of the disease-free equilibrium point is presented in Figure 3. The basic reproduction number is calculated and we obtain $R_0=0.0867<1.$ Figure 3 shows that with $NA=(100,5,25,35,20,16)$ as initial values, the model solutions lead to the no-infection equilibrium point $(K^0).$ We can see that for a long time, there is no HIV infection in the population. This numerical solution is consistent with the analytical results. If $R_0<1$ then the free-of-infection equilibrium point $(K^0)$ is locally asymptotically stable. As demonstrated by mathematical analysis, the solutions are bounded and positive, as is evident. However, we found that when $R_0=0.5077<1,$ the model was unstable disease-free equilibrium point $(K^0).$ In Figure 4, we see the model converges to the endemic equilibrium point. In this case, bi-stability occurs.

Figure 3. The solution of the model converges to a disease-free equilibrium point when $R_0=0.0867<1.$

Figure 4. The solution of the model converges to endemic equilibrium point when $R_0=0.5077<1.$

In addition, we simulate multiple parameters ${\beta }_2$ in order to provide effect of varying ${\beta }_2$ to the basic reproduction number $R_0$ that complies with the sensitivity analysis outcome. The rate of transmission from S to $I_2,$ or ${\beta }_2,$ was discovered to be the most crucial component of both the bifurcation parameter and the HIV/AIDS model. The dynamics of the $I_1$ and $I_2$ subpopulations when we employ a different value of parameter ${\beta }_2$ are displayed in Figure 5. The graph illustrates how raising ${\beta }_2$ results in rising $R_0.$ This indicates that ${\beta }_2$ is involved in the HIV/AIDS epidemic. As can be seen in Figure 5b, $I_2$ is greater than $I_1$ in Figure 5a. Consequently, reducing the rate of transmission from S to $I_2,$ for example, by refraining from HIV/AIDS-related sexual contact (Figure 6).

Figure 5. The dynamics of $I_1$ and $I_2$ subpopulations with different values of ${\beta }_2$ and $R_0.$

Figure 6. The dynamics of T and a subpopulations with different values of ${\beta }_2$ and $R_0.$

Figure 7 shows the existence of bi-stability phenomena. The condition in this work allows for the existence of a stable free-disease and endemic equilibrium points, one of which is stable (blue line) and the other unstable (red-dashed line), despite the fact that $R_0<1.$ In this instance, reverse bifurcation causes a bi-stability event. The initial state of the population may affect the final state of the population when backward bifurcation occurs. A backward bifurcation parameter now plays a crucial function as a limit condition of a backward bifurcation. In this paper, ${\beta }_2$ is a bifurcation parameter. Using the parameter values presented in Table 1, we calculate the values $P=5\times {10}^5>0$ and $Q=3\times {10}^1>0$ where it satisfies the requirement for backward bifurcation as in the analysis. We calculate P and Q by using MAPLE software. As can be observed, system (1)(1) $$\begin{array}{cc}\frac{dS}{dt}\hfill & =\Lambda -{\beta }_{1}S{I}_1-{\beta }_{2}S{I}_2-aS,\hfill \\ \frac{dE}{dt}\hfill & =\eta S-{\beta }_{3}E{I}_1-{\beta }_{4}E{I}_2-dE,\hfill \\ \frac{d{I}_1}{dt}\hfill & ={\beta }_{1}S{I}_1+{\beta }_{3}E{I}_1+{\alpha }_{1}T-b{I}_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\beta }_{2}S{I}_2+{\beta }_{4}E{I}_2-c{I}_2,\hfill \\ \frac{dT}{dt}\hfill & =k_1{I}_1+k_3{I}_2-eT,\hfill \\ \frac{dA}{dt}\hfill & =k_2{I}_2+{\alpha }_{2}T-fA,\hfill \\ \frac{dR}{dt}\hfill & =\mu S-dR,\hfill \end{array}$$(1) tends to two different final states, two stable equilibrium points, namely free-disease (FDE) and endemic (EE) equilibrium points, when ${\beta }_2<{\beta }_2^{\ast },$ for the same set of parameter values. Related phenomena of a backward bifurcation can be seen in Fatimah et al. (2021).

Figure 7. Backward bifurcation of HIV/AIDS model.

4. Conclusion

The HIV/AIDS model has seven compartments: the uneducated (S), the educated (E), the HIV-positive who are taking ARV ($I_1$), the HIV-positive who are not taking ARV ($I_2$), the people receiving an ARV treatment (T), the people with AIDS who are not receiving any treatment (full-blown AIDS (A), and the recovered (R). The interactions between S, E, $I_1,$ and $I_2$ are taken into account in this model. There are both disease-free points and endemic equilibrium points. Whether the equilibrium points are locally stable was investigated using the Routh-Hurwitz criteria. Bi-stability exists when the basic reproduction number is $(R_0<1).$ The basic reproduction number below one is insufficient to achieve a disease-free condition because HIV/AIDS is continuously spreading.

The center manifold theory can be used to identify the model’s backward bifurcation. The bifurcation parameter now plays a crucial function as a limit condition of a backward bifurcation. The basic reproduction number was equalized to one before being used to determine the bifurcation parameter value. The result is the rate of HIV transmission from uneducated people to HIV-positive people who do not take antiretroviral therapy (ART) was found as the most important parameter by the HIV/AIDS model. The bifurcation parameter is ${\beta }_2.$

We calculated the sensitivity index values using a sensitivity analysis to obtain parameters that significantly contribute to HIV/AIDS transmission. Again, ${\beta }_2$ (the rate of transmission from S to $I_2$) was found to be the most important factor in the HIV/AIDS model and bifurcation parameter as well. Therefore, lowering the rate of transmission from S to $I_2,$ for instance by avoiding sexual contact between HIV/AIDS infection subpopulations could be a successful intervention to stop the spread of HIV virus. Backward bifurcation is a key component of many mathematical models and is essential for creating control schemes (Buonomo & Vargas-De-León, 2013; Fatimah et al., 2021). For further work, we add some interventions to the model to stop the spread of HIV, by using a control optimal approach.

The dynamical analysis and proposed model suggest that community social health professionals can raise public awareness of HIV/AIDS to prevent the spread of HIV/AIDS. This can be accomplished through various outreach and education initiatives, including campaigns, counseling, and informational media releases. Comprehensive education and outreach programs must include information on HIV/AIDS transmission, prevention, and treatment.

Acknowledgments

The authors would like to thank Brawijaya University for giving support through Hibah Doktor 2021, No: 1629/UN10.F09/PN/2021.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

Universitas Brawijaya.