Understanding HIV/AIDS dynamics: insights from CD4+T cells, antiretroviral treatment, and country-specific analysis

In this article, we present a mathematical model for human immunodeficiency virus (HIV)/Acquired immune deficiency syndrome (AIDS), taking into account the number of CD4+T cells and antiretroviral treatment. This model is developed based on the susceptible, infected, treated, AIDS (SITA) framework, wherein the infected and treated compartments are divided based on the number of CD4+T cells. Additionally, we consider the possibility of treatment failure, which can exacerbate the condition of the treated individual. Initially, we analyze a simplified HIV/AIDS model without differentiation between the infected and treated classes. Our findings reveal that the global stability of the HIV/AIDS-free equilibrium point is contingent upon the basic reproduction number being less than one. Furthermore, a bifurcation analysis demonstrates that our simplified model consistently exhibits a transcritical bifurcation at a reproduction number equal to one. In the complete model, we elucidate how the control reproduction number determines the stability of the HIV/AIDS-free equilibrium point. To align our model with the empirical data, we estimate its parameters using prevalence data from the top four countries affected by HIV/AIDS, namely, Eswatini, Lesotho, Botswana, and South Africa. We employ numerical simulations and conduct elasticity and sensitivity analyses to examine how our model parameters influence the control reproduction number and the dynamics of each model compartment. Our findings reveal that each country displays distinct sensitivities to the model parameters, implying the need for tailored strategies depending on the target country. Autonomous simulations highlight the potential of case detection and condom use in reducing HIV/AIDS prevalence. Furthermore, we identify that the quality of condoms plays a crucial role: with higher quality condoms, a smaller proportion of infected individuals need to use them for the potential eradication of HIV/AIDS from the population. In our optimal control simulations, we assess population behavior when control interventions are treated as time-dependent variables. Our analysis demonstrates that a combination of condom use and case detection, as time-dependent variables, can significantly curtail the spread of HIV while maintaining an optimal cost of intervention. Moreover, our cost-effectiveness analysis indicates that the condom use intervention alone emerges as the most cost-effective strategy, followed by a combination of case detection and condom use, and finally, case detection as a standalone strategy.

1 Introduction

Human immunodeficiency virus (HIV) is a virus that infects CD4+ T lymphocytes, leading to a weakened immune system in individuals. On the other hand, Acquired immune deficiency syndrome (AIDS) refers to the symptoms that occur as a result of a weakened immune system due to HIV infection (1). HIV can be transmitted through bodily fluids such as blood, semen, genital fluids, and breast milk. According to data from the WHO and the United Nations Programme on HIV/AIDS (UNAIDS), in 2016, there were 36.7 million people living with HIV (PLHIV)/AIDS worldwide. The majority of people living with HIV are in low and middle-income countries, such as Eswatini, where the prevalence is 25.9% among the adult population, Lesotho with a prevalence of 19.3%, South Africa with a prevalence of 17.8%, Botswana with a prevalence of 16.4%, and Mozambique with a prevalence of 11.6% (2). According to the same source, the top 10 countries with the highest HIV prevalence include 10 African nations, with South Africa having the highest number of people living with HIV, surpassing 7.6 million in 2020. Beyond Africa, the spread of HIV is also a significant concern, particularly in Indonesia, where the latest data report 519,158 individuals affected by HIV as of June 2022 (3). This places Indonesia as the third-highest country with people living with HIV in Asia, following India and Thailand.

The HIV illness may be divided into four phases. People living with HIV (PLHIV) and have entered stage one may experience mild symptoms such as flu-like symptoms, diarrhea, and fever. Those who have entered stage two may experience symptoms such as Tuberculosis (TB), swollen lymph nodes, and skin disorders. In stage three, individuals may experience symptoms in the mucous membranes, such as Tuberculosis (TB) in the lymph nodes. Finally, in stage four, individuals may experience systemic meningoencephalitis. Stage four is commonly referred to as AIDS. HIV attacks CD4+ T cells in the bodies of infected individuals. CD4+ T cells play a crucial role in the immune system and perform many functions in immune activation, coordination, modulation, and regulation (4). AIDS is defined as an outcome among those living with HIV if the CD4+ T cell count is < 200.

Antiretroviral therapy (ART) is available as a treatment to reduce the risk of HIV transmission and lower the amount of HIV in the blood. Highly active antiretroviral therapy (HAART) has been successful in reducing morbidity and mortality among individuals infected with HIV (5). According to the United Nations Programme on HIV/AIDS (UNAIDS),¹ there were 1.3 million people newly infected with HIV in 2022. Among 39 million people living with HIV in 2022, 630,000 people died of AIDS. According to the global HIV and AIDS epidemic,² only 86% of people with HIV globally knew their HIV status in 2022, the rest of them do not know that they were living with HIV. Therefore, antiretroviral therapy has the potential to reduce the number of individuals infected with HIV in sub-Saharan Africa and other countries.

The mathematical models have been used widely to model the spread of diseases, such as dengue (6), malaria (7, 8), tuberculosis (9–11), HIV (12), COVID-19 (13–16), pneumonia (17). In the context of the HIV/AIDS infectious disease spread, the mathematical model can help researchers understand the impact of interventions and can be used to predict the potential outcome of scenarios in the field. An early mathematical model for HIV/AIDS can be found in Rahman et al.'s study (18). In 2009, Mukandavire et al. (19) presents a mathematical model for HIV/AIDS transmission dynamics, incorporating an explicit incubation period. It demonstrates that effective public health educational campaigns, when targeted at both sexually immature and sexually mature individuals, can significantly slow down the epidemic. The study also identifies the presence of backward bifurcation in the mathematical model, highlighting the complexity of disease dynamics and the potential impact of comprehensive intervention strategies. The presence of the backward bifurcation phenomenon in their model suggests that the extinction of HIV may not solely depend on the condition of the reproduction number being less than one. This is because another endemic possibility may emerge even when the reproduction number is already less than one. In lay terms, backward bifurcation in the model means that even if the conditions initially seem favorable for controlling and reducing HIV (for example, when the reproduction number is less than one), there is still a risk of a resurgence or a persistent presence of the infection. This phenomenon introduces a level of complexity, suggesting that the effectiveness of interventions and the possibility of HIV extinction are not solely determined by one factor (such as a low reproduction number). Instead, additional factors or conditions may influence the dynamics of the infection, making it more challenging to predict and control. Nyabadza and Mukandavire in (20) employ ordinary differential equations to investigate the dynamics of HIV/AIDS, particularly in the context of HIV testing and screening campaigns. The key findings include that having a basic reproduction number below one is necessary but not sufficient for disease control due to backward bifurcation phenomena. Additionally, the study fits the model to real data on HIV prevalence in South Africa, adding empirical validation to the model's insights that HIV counseling and testing itself has very little impact in reducing the prevalence of HIV unless the efficacy of the campaigns exceeds an evaluated threshold in the absence of backward bifurcation. Recently, Zhai et al. (21) develop a stochastic HIV/AIDS model that considers individuals with protection awareness. Their research revealed that HIV can become extinct when the $R_0^s$ is <1. Furthermore, the study highlights that enhancing the protection efficiency of individuals with awareness and the implementation of continuous antiretroviral therapy (ART) both contribute to reducing the number of people living with HIV (PLHIV), offering potential strategies for disease control. In Jamil et al.'s study (22), a fractal fractional HIV/AIDS model is introduced, using fractional order differential equations. The study utilizes the first and second derivatives of the Lyapunov function to conduct a global analysis of the model. The research suggests that measures to reduce the effective contact rate between susceptible and infected individuals, coupled with improved treatment for those who are already infected, can enhance the effectiveness of interventions against HIV/AIDS. Recently, due to the COVID-19 pandemic, many mathematical models have been introduced to understand the impact of co-infection between HIV/AIDS with COVID-19. Xu et al. (23) focus on developing a mathematical model for HIV-TB co-infection dynamics and validate it using real incidence data from different regions. The article also delves into the comparison of numerical schemes to determine the most effective computational approach for simulating the model. Pinto et al. (24) explore models for HIV and TB coinfection dynamics, considering both fractional and integer order models. Their analysis encompasses treatment strategies for both diseases and includes considerations for the vertical transmission of HIV. Ringa et al. (25) conduct an analysis on sub-models (and co-infection model) related to HIV and COVID-19 co-infection. They apply an optimal control approach to control variables, finding that preventive measures can substantially reduce the burden of co-infections with COVID-19, and effective treatment of COVID-19 could, in turn, reduce co-infections with opportunistic infections such as HIV/AIDS. Please see Omami et al.'s study (26) for another model which incorporates a coinvection between HIV, dengue, and COVID-19. Another classical model was presented by Garba et al. (27), wherein they employed a mathematical framework to simulate the spread of HIV in Nigeria. Their model incorporates factors such as condom use and asymptomatic cases. The researchers utilized incidence data from Nigeria to calibrate the parameters of their model. Readers may refer to (28–34) for more HIV/AIDS related models.

This research is an extension of the study conducted by Rahman et al. (18) with modifications that include the addition of an AIDS compartment as a variable, taking into consideration the population engaged in unprotected sexual intercourse, and the intervention of antiretroviral therapy. This research aims to develop a comprehensive model for the spread of HIV, integrating the dynamics of CD4+T cells and considering key interventions such as condom use and treatment. By fitting model parameters based on data from four top countries with high HIV incidence rates, we seek to provide a nuanced understanding of how these factors influence the trajectory of the epidemic. Additionally, sensitivity analysis was conducted to assess the impact of various model parameters on the number of cases in each country, offering insights into the relative importance of different factors. Furthermore, optimal control techniques are employed to forecast potential optimal scenarios in the field, aiding in the design of effective strategies for HIV prevention and management.

In this study, a mathematical model for the spread of HIV/AIDS through unprotected sexual intercourse has been constructed based on the classification of the number of CD4+ T cells in the body, incorporating the intervention of antiretroviral therapy. The number of CD4+ T cells is crucial in constructing an HIV mathematical model because these cells play a central role in the immune system, and their depletion is a hallmark of HIV infection. CD4+ T cells are a type of white blood cell that helps coordinate the immune response to infections. HIV specifically targets and infects these cells, leading to a decline in their numbers over time. Incorporating the CD4+ T cell count into the mathematical model allows for a more realistic representation of the dynamics of HIV infection and disease progression. Applying some mathematical analysis tools, we conduct an analytical study on our model, including an analysis of existence, the stability analysis of equilibrium points, and the analysis of the basic reproduction number ℛ₀. With this study, we can understand the long-term behavior of HIV transmission as it changes over time. The analysis of the basic reproduction number can quantify which factors play a significant role in efforts to control the spread of HIV. Following that, we conducted numerical simulations, which involved analyzing the elasticity and sensitivity of ℛ₀, in addition to performing autonomous and optimal control simulations using the constructed model. The goal is to gain insights into the impact of antiretroviral therapy on the transmission of HIV/AIDS through unprotected sexual intercourse. This analysis is based on the classification of the number of CD4+ T cells in the body. By interpreting the outcomes from both the analytical study and numerical simulations, we aim to better understand the dynamics and effects of the treatment and condom use on the spread of the virus.

The article is structured as follows: Section 2 presents the construction and analysis of a simple case model where the number of CD4+T cells is not considered in the model. This section includes a global stability of equilibrium points and a sensitivity analysis of the basic reproduction number. Section 3 delves into the mathematical model analysis of the complete model, discussing equilibrium points, control reproduction number, data fitting, and model sensitivity analysis. Section 4 describes the modification of the complete model into a control optimal model. This section includes numerical experiments for different strategies, as well as a cost-effectiveness analysis. Finally, Section 5 provides conclusion.

2 A mathematical model of HIV/AIDS with antiretroviral treatment without considering the number of CD4+T cell number class

2.1 Model construction

In this section, we consider an antiretroviral treatment for an HIV-infected individual. First, let us consider the total human population (aged 15–49 years), denoted by N, to be be categorized into four different compartments based on their health statuses, namely, the susceptible individual (S), PLHIV with infection only (I), PLHIV receiving treatment (T), and PLHIV with AIDS illness (A). In this model, we assume that without any test, the HIV-infected individual cannot be detected. Hence, they cannot get a proper treatment. To construct the model, we use the transmission diagram as shown in Figure 1.

Figure 1

Figure 1. HIV transmission diagram with antiretroviral treatment.

The construction of the model is as follows. The recruitment rate of N is always considered as a susceptible person from a younger age class of <15 years, with a constant rate Λ. Infection only occurs due to successful contact between susceptible individuals with infected individuals I and T with a probability of β_u and β_t, respectively. The most effective way to prevent sexual transmission of HIV is through abstinence, that is, avoiding all vaginal, anal, and oral sex. In our model, we include the use of a male condom or a female condom that can avoid the transmission of HIV. The variables ϵ and κ are denoted as condom efficacy and the proportion of people who use condoms, then the transmission rate β_u and β_t can be reduced by a factor of 1 − ϵκ. With this assumption, we have κ ∈ [0, 1] and ϵ ∈ [0, 1]. A bigger value of κ represents a better quality of a condom, and a bigger value of ϵ represents a condition where more people use condoms during their sexual activity. With the mass action infection process, then the infection force of infection is given by (1 − ϵκ)(β_uI + β_tT). In this model, we have a transition from I to T due to a medical test for HIV detection, which is denoted by τ. Furthermore, the transition from the HIV stage to AIDS is given by γ_u and γ_t for I and T individuals, respectively. Due to antiretroviral treatment, we assume that γ_t < γ_u. Each compartment has a natural death rate denoted by μ, except for A compartment which has an additional death rate due to AIDS, which is denoted by η. Based on this assumption, the complete model of HIV/AIDS transmission with the usage of condoms and the intervention of antiretroviral treatment is given by the following system of ordinary differential equations.

$\begin{array}{ll}\begin{array}{l}\frac{dS}{dt}=\Lambda -((1-ϵ\kappa )({\beta }_{u}I+{\beta }_{t}T)+\mu )S,\\ \frac{dI}{dt}=(1-ϵ\kappa )({\beta }_{u}I+{\beta }_{t}T)S-(\tau +{\gamma }_u+\mu )I,\\ \frac{dT}{dt}=\tau I-({\gamma }_t+\mu )T,\\ \frac{dA}{dt}={\gamma }_{u}I+{\gamma }_{t}T-(\mu +\eta )A,\end{array}& (1)\end{array}$

completed with the following initial condition

$S(0)>0,I(0)\ge 0,T(0)\ge 0,A(0)\ge 0.$

2.2 Preliminary analysis

In this section, we describe two theorems to guarantee the positiveness of the solution and also the positive invariant region of the system (1).

Theorem 1. All solutions of the HIV/AIDS model in Equation (1) with a non-negative initial condition in ${ℝ}_{+}^4$ remain positive for all time t > 0.

Proof: From the first equation on the (1), we have

$\begin{array}{l}\frac{dS}{dt}=\Lambda -((1-ϵ\kappa )({\beta }_{u}I+{\beta }_{t}T)+\mu )S\\ >-((1-ϵ\kappa )({\beta }_{u}I+{\beta }_{t}T)+\mu )S.\end{array}$

The solution of S(t) is given by

$S(t)>S(0)exp\left[-\left({\displaystyle {\int }_0^t}(1-ϵ\kappa )({\beta }_{u}I(\varphi )+{\beta }_{t}T(\varphi ))+\mu t\right)\right].$

Since S(0) > 0, then S(t) is always positive for all t > 0. The remaining variables I(t), T(t), and A(t) can be shown in a similar way. Hence the solution set {S, I, T, A} is always non-negative for all time t > 0.□

Theorem 2. The region

$\mathcal{D}=\left\{(S,I,T,A)\in {ℝ}_{+}^4:0<S+I+T+A\le max\left(\frac{\Lambda }{\mu },N(0)\right)\right\}$

is positively invariant and attracting with respect to the system (1) with a non-negative initial condition in ${ℝ}_{+}^4$.

Proof: Adding all equations in the system (1), we have

$\frac{dN}{dt}=\Lambda -\mu N-\eta A\le \Lambda -\mu N.$

Solving this differential equation with respect to N(t), we have

$N(t)\le N(0)exp(-\mu t)+\frac{\Lambda }{\mu }(1-exp(-\mu t)).$

Therefore, we can see that $N(t)\to \frac{\Lambda }{\mu }$ for t → ∞. If $N(0)>\frac{\Lambda }{\mu }$, then N(t) will monotonically decrease and tends to $\frac{\Lambda }{\mu }$. If $N(0)<\frac{\Lambda }{\mu }$, then N(t) will monotonically increase and tends to $\frac{\Lambda }{\mu }$. On the other hand, if $N(0)=\frac{\Lambda }{\mu }$, then N(t) will remain constant for all time t, where $N(t)=\frac{\Lambda }{\mu }$. Hence, according to Theorem 1 and previous calculation, we have $0<S+I+T+A\le max\left\{\frac{\Lambda }{\mu },N(0)\right\}$.□

Based on Theorems 1 and 2, it is sufficient to consider the dynamics of the system (1) in the region $\mathcal{D}$ where the existence, uniqueness, and positiveness of the solution hold.

2.3 The HIV/AIDS-free equilibrium point and the basic reproduction number of the model (1)

In this section, we analyze the existence and stability of the HIV/AIDS-free equilibrium points, and how they relate to the basic reproduction number of the system (1). The HIV/AIDS-free equilibrium point of the system (1) is given by

$\begin{array}{l}{\mathcal{E}}_1=(S^0,I^0,T^0,A^0)=\left(\frac{\Lambda }{\mu },0,0,0\right).\end{array}$

Before we calculate the local and global stability criteria of the HIV/AIDS-free equilibrium points, we first calculate the basic reproduction number of the system (1). First, we decompose the Jacobian matrix of an infected subsystem of the system (1) which is evaluated in ${\mathcal{E}}_0$ in a transition Σ and transmission T matrix as follows:

$\begin{array}{r}\text{T}=\left[\begin{array}{cc}\frac{(1-ϵ\kappa ){\beta }_u\Lambda }{\mu }& \frac{(1-ϵ\kappa ){\beta }_t\Lambda }{\mu }\\ 0& 0\end{array}\right],\\ \Sigma =\left[\begin{array}{cc}-\tau -{\gamma }_u-\mu & 0\\ \tau & -{\gamma }_t-\mu \end{array}\right].\end{array}$

Since the elements of the second row of T are all zero, then defining E = [10]^t, the next-generation matrix is given by

$\begin{array}{l}\begin{array}{l}\text{K}=-{\text{E}}^t\text{T}{\Sigma }^{-1}\text{E}=\left[\frac{\Lambda ({\beta }_u({\gamma }_t+\mu )+{\beta }_t\tau )(1-ϵ\kappa )}{\mu (\tau +{\gamma }_u+\mu )({\gamma }_t+\mu )}\right].\end{array}\end{array}$

The basic reproduction number of the system (1) is taken by the spectral radius of K, which is given by

$\begin{array}{l}{\mathcal{R}}_0=\frac{\Lambda ({\beta }_u({\gamma }_t+\mu )+{\beta }_t\tau )(1-ϵ\kappa )}{\mu (\tau +{\gamma }_u+\mu )({\gamma }_t+\mu )}.\end{array}$

The basic reproduction number in many epidemiological models play an important role in determining the local and global stability of the equilibrium points of their model. The basic reproduction number in our model represents the total number of secondary cases of HIV caused by one primary case of HIV in a completely virgin population during its infection period. The following theorem gives the local stability criteria of ${\mathcal{E}}_0$.

Theorem 3. The HIV/AIDS-free equilibrium point ${\mathcal{E}}_1$ is locally asymptotically stable if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$.

Proof: The Jacobian matrix of the system (1) evaluated in ${\mathcal{E}}_0$ is given by

$\mathcal{J}({ϵ}_0)=\left[\begin{array}{cccc}-\mu & -\frac{(1-ϵ\kappa ){\beta }_u\Lambda }{\mu }& -\frac{(1-ϵ\kappa ){\beta }_t\Lambda }{\mu }& 0\\ 0& \frac{(1-ϵ\kappa ){\beta }_u\Lambda }{\mu }-\mu -\tau -{\gamma }_u& \frac{(1-ϵ\kappa ){\beta }_t\Lambda }{\mu }& 0\\ 0& \tau & -\mu -{\gamma }_t& 0\\ 0& {\gamma }_u& {\gamma }_t& -\mu -\eta \end{array}\right],$

which has two explicit eigenvalues, i.e, λ₁ = − μ and λ₂ = − (μ + η), while the other two eigenvalues are taken by the solution of the following polynomial

$\begin{array}{c}f(\lambda )=\mu {\lambda }^2+(2\mu +\tau +{\gamma }_u+{\gamma }_t)(1-{\mathcal{R}}_1)\lambda +\\ \mu (\tau +{\gamma }_u+\mu )({\gamma }_t+\mu )(1-{\mathcal{R}}_0)=0,\end{array}$

where ${\mathcal{R}}_1=\frac{\Lambda {\beta }_u(1-ϵ\kappa )}{\mu (2\mu +\tau +{\gamma }_u+{\gamma }_t)}$ and ${\mathcal{R}}_0=\frac{\Lambda ({\beta }_u({\gamma }_t+\mu )+{\beta }_t\tau )(1-ϵ\kappa )}{\mu (\tau +{\gamma }_u+\mu )({\gamma }_t+\mu )}$. Since ${\mathcal{R}}_1<{\mathcal{R}}_0$, then the other two eigenvalues will have a negative real part if ${\mathcal{R}}_0<1$. Therefore, the HIV/AIDS-free equilibrium ${\mathcal{E}}_0$ is locally asymptotically stable if ${\mathcal{R}}_0<1$. □

The following Theorem gives the global stability criteria of ${\mathcal{E}}_0$.

Theorem 4. The HIV/AIDS-free equilibrium ${\mathcal{E}}_1$ of the system (1) is globally stable if ${\mathcal{R}}_0<1$.

Proof: Using a same approach as authors in (35), let us define

$\begin{array}{l}\begin{array}{l}\frac{dX}{dt}=F(X,Z),\\ \frac{dZ}{dt}=G(X,Z),G(X,0)=0,\end{array}\end{array}$

where X = S ∈ ℝ₊ is the compartment of non-infected individuals, and $Z=(I,T,A)\in {ℝ}_{+}^3$ is the infected compartments. Let $X^{*}=(\frac{\Lambda }{\mu },0)$.

From direct calculation, we have

$\begin{array}{l}F(X,0)=\left[\Lambda -\mu S\right],\\ G(X,Z)=AZ-Ĝ(X,Z),\end{array}$

where

$\begin{array}{l}A=D_{Z}G(X^{*},0)=\\ \left[\begin{array}{ccc}\frac{(1-ϵ\kappa ){\beta }_u}{\mu }-\tau -{\gamma }_u-\mu & \frac{(1-ϵ\kappa ){\beta }_t}{\mu }& 0\\ \tau & -{\gamma }_t-\mu & 0\\ {\gamma }_u& {\gamma }_t& -\mu -\eta \end{array}\right],\\ Ĝ(X,Z)=\left[\begin{array}{c}{Ĝ}_1(X,Z)\\ {Ĝ}_2(X,Z)\\ {Ĝ}_3(X,Z)\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}(1-ϵ\kappa )({\beta }_{u}I+{\beta }_{t}T)(\frac{\Lambda }{\mu }-S)\\ 0\\ 0\end{array}\end{array}\right].\end{array}$

Since I and T are always positive, then it is clear that Ĝ(X, Z) ≥ 0. It is also clear that X^* = (Λ/μ) is globally stable for F(X, 0). Hence, ${\mathcal{E}}_0$ is globally asymptotically stable.□

2.4 HIV/AIDS endemic equilibrium point of the system (1)

The HIV/AIDS endemic equilibrium point of system (1) is given by

${\mathcal{E}}_2=(S^{*},I^{*},T^{*},A^{*}),$

where

$\begin{array}{l}{S}^{*}=\frac{\left({\gamma }_t+\mu \right)\left(\tau +{\gamma }_u+\mu \right)}{\left({\beta }_u\left({\gamma }_t+\mu \right)+{\beta }_t\tau \right)\left(1-ϵ\kappa \right)},\\ =\frac{\Lambda }{{ℛ}_0\mu },\\ I^{*}=\frac{\Lambda (1-ϵ\kappa )({\beta }_u(\mu +{\gamma }_t)+{\beta }_t\tau )-\mu (\mu +{\gamma }_t)(\mu +\tau +{\gamma }_u)}{((1-ϵ\kappa )({\beta }_u(\mu +{\gamma }_t)+{\beta }_t\tau )(\mu +\tau +{\gamma }_u)}\\ =\left(1-\frac{1}{{ℛ}_0}\right)\frac{\Lambda }{\mu +\tau +{\gamma }_u},\\ T^{*}=\frac{\tau }{\mu +{\gamma }_t}\frac{\Lambda (1-ϵ\kappa )({\beta }_u(\mu +{\gamma }_t)+{\beta }_t\tau )-\mu (\mu +{\gamma }_t)(\mu +\tau +{\gamma }_u)}{((1-ϵ\kappa )({\beta }_u(\mu +{\gamma }_t)+{\beta }_t\tau )(\mu +\tau +{\gamma }_u)}\\ =\left(1-\frac{1}{{ℛ}_0}\right)\frac{\Lambda \tau }{(\mu +{\gamma }_t)(\mu +\tau +{\gamma }_u)},\\ A^{*}=\frac{\Lambda (1-ϵ\kappa )({\beta }_u(\mu +{\gamma }_t)+{\beta }_t\tau )-\mu (\mu +{\gamma }_t)(\mu +\tau +{\gamma }_u)}{((1-ϵ\kappa )({\beta }_u(\mu +{\gamma }_t)+{\beta }_t\tau )(\mu +\tau +{\gamma }_u)}\frac{(\mu {\gamma }_u+{\gamma }_t(\tau +{\gamma }_u)}{(\mu +{\gamma }_t)(\eta +\mu ))}\\ =\left(1-\frac{1}{{ℛ}_0}\right)\frac{\Lambda \left(\mu {\gamma }_u+{\gamma }_t\left(\tau +{\gamma }_u\right)\right)}{(\mu +{\gamma }_t)(\mu +\tau +{\gamma }_u)(\eta +\mu )}.\end{array}$

From the expression of ${\mathcal{E}}_2$, we have the following theorem.

Theorem 5. The HIV/AIDS endemic equilibrium ${\mathcal{E}}_2$ of the system (1) exist in ${ℝ}_{+}^4$ if ${\mathcal{R}}_0>1$.

In the following theorem, we show the non-existence of backward-bifurcation phenomena of the system (1).

Theorem 6. The system (1) always exhibits a forward bifurcation phenomenon at ${\mathcal{R}}_0=1$.

Proof: To analyze the bifurcation phenomena of system (1), we use the well-known Castillo-Song bifurcation theorem (35). Please see (36–38) for more examples of the implementation of this theorem in another epidemiological model. First, let us define S = x₁, I = x₂, T = x₃, and A = x₄ and g_i for i = 1, ..., 4 represent $\frac{dS}{dt},\frac{dI}{dt},\frac{dT}{dt},$ and $\frac{dA}{dt}$, respectively. Next, we choose β_t as the bifurcation parameter. Solving ${\mathcal{R}}_0=1$ with respect to β_t, we have

$\begin{array}{l}{\beta }_t^{*}=\frac{\mu (\tau +{\gamma }_u+\mu )({\gamma }_t+\mu )}{\Lambda (1-ϵ\kappa )\tau }-\frac{{\beta }_u({\gamma }_t+\mu )}{\tau }.\end{array}$

Next, we calculate the Jacobian matrix of the system (1) and evaluate it at ${\beta }_t^{*}$ and ${\mathcal{E}}_0$, we have

$\begin{array}{l}\mathcal{A}=\left[\begin{array}{cccc}-\mu & -\frac{(1-ϵ\kappa ){\beta }_u\Lambda }{\mu }& \frac{(\Lambda (1-ϵ\kappa ){\beta }_u-{\mu }^2-\mu \tau -\mu {\gamma }_u)(\mu +{\gamma }_t)}{\mu \tau }& 0\\ 0& \frac{(1-ϵ\kappa ){\beta }_u\Lambda }{\mu }-\mu -\tau -{\gamma }_u& -\frac{(\Lambda (1-ϵ\kappa ){\beta }_u-{\mu }^2-\mu \tau -\mu {\gamma }_u)(\mu +{\gamma }_t)}{\mu \tau }& 0\\ 0& \tau & -\mu -{\gamma }_t& 0\\ 0& {\gamma }_u& {\gamma }_t& -\eta -\mu \end{array}\right]\end{array}$

The eigenvalues of $\mathcal{A}$ are

$\begin{array}{c}{\lambda }_1=0,{\lambda }_2=-\mu ,{\lambda }_3=-(\mu +\eta ),\text{ and }{\lambda }_4=\\ \frac{\Lambda (1-\kappa ϵ){\beta }_u-\mu (2\mu +\tau +{\gamma }_t+{\gamma }_u)}{\mu }.\end{array}$

We can see clearly that we have a simple zero eigenvalue, and λ₂ and λ₃ are negative. We have λ₄ < 0 if and only if $\frac{\Lambda (1-\kappa ϵ){\beta }_u}{\mu (2\mu +\tau +{\gamma }_t+{\gamma }_u)}<1$. Since

$\begin{array}{c}{\mathcal{R}}_0=1=\frac{\Lambda ({\beta }_u({\gamma }_t+\mu )+{\beta }_t\tau )(1-ϵ\kappa )}{\mu (\tau +{\gamma }_u+\mu )({\gamma }_t+\mu )}>\frac{c\Lambda {\beta }_u}{\mu (\mu +\tau +{\gamma }_u)}>\\ \frac{\Lambda (1-\kappa ϵ){\beta }_u}{\mu (2\mu +\tau +{\gamma }_t+{\gamma }_u)},\end{array}$

then we also have λ₄ < 0. Hence, we can proceed to the next step to analyze the bifurcation type of our model. The bifurcation type of our model can be determined with the following formula:

$\begin{array}{ll}a={\displaystyle \sum _{k,i,j=1}^4}{v}_k{w}_i{w}_j\frac{{\partial }^2{g}_k}{d{x}_{i}d{x}_j}(\text{0},0),& b={\displaystyle \sum _{k,i=1}^4}{v}_k{w}_i\frac{{\partial }^2{g}_k}{d{x}_{i}d{\beta }_t}(\text{0},0),\end{array}$

where v and w is the left and right eigenvectors of $\mathcal{A}$ with respect to the zero eigenvalue, respectively. The left eigenvalue of $\mathcal{A}$ with respect to the eigenvalue 0 is given by

$\begin{array}{r}\begin{array}{r}{v}_1=0,v_2=1>0,v_3=\\ \frac{(\Lambda (1-\kappa ϵ){\beta }_u-\mu (\mu +\tau +{\gamma }_u))}{\tau \mu }<0,\text{ with }{v}_4=0.\end{array}\end{array}$

On the other hand, the right eigenvector of $\mathcal{A}$ with respect to the eigenvalue 0 is given by

$\begin{array}{r}\begin{array}{r}{w}_1=-\frac{({\gamma }_t+\mu )(\tau +{\gamma }_u+\mu )}{\mu \tau },w_2=\frac{{\gamma }_t+\mu }{\tau },w_3=1,\text{ and}\\ w_4=\frac{\mu {\gamma }_u+\tau {\gamma }_t+{\gamma }_t{\gamma }_u}{(\eta +\mu )}.\end{array}\end{array}$

Hence, we have

$\begin{array}{l}a={\displaystyle \sum _{k,i,j=1}^4}{v}_k{w}_i{w}_j\frac{{\partial }^2{g}_k}{d{x}_{i}d{x}_j}(\text{0},0)\\ =-\frac{2(\tau +{\gamma }_u+\mu )({\beta }_t\tau +{\beta }_u({\gamma }_t+\mu ))({\gamma }_t+\mu )(1-ϵ\kappa )}{{\tau }^2\Lambda },\\ b={\displaystyle \sum _{k,i=1}^4}{v}_k{w}_i\frac{{\partial }^2{g}_k}{d{x}_{i}d{\beta }_t}(\text{0},0)\\ =(1-ϵ\kappa ).\end{array}$

Since a < 0 and b > 0, the system (1) always exhibits a forward bifurcation at ${\mathcal{R}}_0=1$. □

The following corollary is a direct implication from Theorem 6.

Corollary 1. The HIV/AIDS endemic equilibrium ${\mathcal{E}}_2$ of the system (1) is locally stable for ${\mathcal{R}}_0>1$ but close to one.

In epidemiology, bifurcation refers to a qualitative change in the behavior of a dynamical system as a parameter is varied. Forward bifurcation specifically occurs when a stable endemic equilibrium coexists with an unstable disease-free equilibrium at parameter values where the basic reproduction number (${\mathcal{R}}_0$) is larger than one. On the other hand, in cases where the reproduction number falls below one, a stable disease-free equilibrium exists without the existence of the endemic equilibrium. Consequently, the condition where the reproduction number equals one marks the bifurcation point. In the context of the HIV/AIDS model that we developed in Equation (1), forward bifurcation has significant implications. It highlights the importance of not only reducing transmission rates but also addressing factors that contribute to the maintenance of stable disease-free equilibrium, such as the efficacy of condom use and treatment. By considering the implications of forward bifurcation, epidemiologists and policymakers can develop more nuanced and effective approaches to combating the HIV/AIDS epidemic.

2.5 Effect of ${\mathcal{R}}_0$ to endemic size ${\mathcal{E}}_2$

Here in this section, we will perform the elasticity index of the basic reproduction number ${\mathcal{R}}_0$ and the endemic equilibrium size of the HIV/AIDS model in Equation (1) and also the parameter sensitivity to ${\mathcal{R}}_0$. We use parameter values for Lesotho (See table in Appendix 2), except for τ and κ, which is given as follows:

$\begin{array}{c}{\beta }_u=1.4888\times 10^{-6},{\beta }_t=1.14888\times 10^{-7},\mu =\frac{1}{57},\\ {\gamma }_u=0.6,{\gamma }_t=0.3,\tau =0.1,ϵ=0.925,\kappa =0.1,\text{ and}\\ \Lambda =\frac{736680}{57}.& (2)\end{array}$

With this set of parameters, we have ${\mathcal{R}}_0=1.43$, which is larger than one. Hence, the HIV/AIDS endemic equilibrium ${\mathcal{E}}_2$ exists and is given by

${\mathcal{E}}_2=(S,I,T,A)=(514872,5423,1707,15709)$

is locally stable. This confirms the results of Theorems 5 and 6.

Next, we calculate the elasticity index of ${\mathcal{R}}_0$ and the endemic equilibrium size. To perform this simulation, we use the following formula (39):

$\begin{array}{l}{E}_Q^p=\frac{\partial Q}{\partial p}\times \frac{p}{Q},\end{array}$

where p is the model parameter, and Q is the quantity of the model output, such as ${\mathcal{R}}_0$ or endemic equilibrium size ${\mathcal{E}}_2$. Using the above formula, we have

$\begin{array}{l}{E}_{{\mathcal{R}}_0}^{{\beta }_u}=\frac{{\beta }_u(\mu +{\gamma }_t)}{{\beta }_u(\mu +{\gamma }_t)+\tau {\beta }_t}.\end{array}$

Substituting parameter values in Equation (2) yield $E_{{\mathcal{R}}_0}^{{\beta }_u}=+0.969$, which means that increasing β_u by 1% will increase ${\mathcal{R}}_0$ by 0.969%. Therefore, if we increase β_u from 1.4888 × 10⁻⁷ to 1.503 × 10⁻⁷ (increased by 1%), then ${\mathcal{R}}_0$ increases from 1.43 to 1.44, which is an increase of 0.969%. With the same approach, we can calculate the elasticity of each parameter in the system (1) (except Λ that we ignored since the number of recruitment rates cannot be changed in the field) with respect to ${\mathcal{R}}_0$ and ${\mathcal{E}}_2$. The results can be seen in Figure 2.

Figure 2

Figure 2. The elasticity of ${\mathcal{R}}_0$, S^*, I^*, T^*, and A^* with respect to (A) β_u (blue), β_t (cyan), γ_u (green), and γ_t (yellow), and with respect to (B) τ (blue), κ (cyan), and ϵ (yellow).

From Figure 2, we can see that β_u has a significant impact on ${\mathcal{R}}_0$ and ${\mathcal{E}}_0$, more dominant than β_t. Hence, infection from untreated individuals plays a significant role in determining the spread of HIV/AIDS. Increasing β_u and β_t will increase ${\mathcal{R}}_0$ and all infected compartments in ${\mathcal{E}}_0$, but reduce S in ${\mathcal{E}}_0$. Furthermore, we can see that the progression to the A compartment, which is presented by parameters γ_u and γ_t, will reduce ${\mathcal{R}}_0$. We can also see that γ_u is more significant in affecting ${\mathcal{R}}_0$ or all compartments in ${\mathcal{E}}_0$ compared to γ_t. Another important result is that the case detection rate τ is promising in reducing ${\mathcal{R}}_0$ and all infected compartments.

Since we can see that the use of condoms to reduce the infection rates β_u and β_t shows a promising potential to control HIV/AIDS, it is necessary to see the combination of ϵ and κ to reduce ${\mathcal{R}}_0$. By substituting all parameter values in Equation (2), except ϵ and κ, into ${\mathcal{R}}_0$, we have ${\mathcal{R}}_0(ϵ,\kappa )=1.591(1-ϵ\kappa ).$ Drawing this function in the ϵ − κ plane, we have the results in Figure 3. It can be seen that if the efficacy of the condom or the proportion of people who use condoms is <0.37, then the basic reproduction number will always be larger than one (Area 1). Hence, the endemic situation will always appear in the population. On the other hand, if the two above mentioned parameters are >0.37 (Area 2), then there is a chance to eradicate HIV from the population (Area 2b).

Figure 3

Figure 3. Dependency of ${\mathcal{R}}_0$ on the condom efficacy and the proportion of people who use it.

3 A mathematical model of HIV/AIDS considering CD4+T cell number

3.1 Model construction

In this section, we modify our previous model in the system (1) by accommodating the number of CD4+T cells in the human body. To construct the model, we use the following assumptions:

1. We classify each of the I and T compartments in the system (1) into three compartments, which present the class of infected and treated individuals based on the number of CD4+T cells. We denote I₁, I₂, and I₃ as infected individuals who have not yet been detected and have different average numbers of CD4+T cells in their blood. Please see the descriptions of I_i and T_i for i = 1, 2, 3 in Table 1.

2. We assume that since HIV infection itself is not very harmful, there is no additional death rate attributable directly to HIV individuals (compartments I_i and T_i). However, additional deaths occur only among individuals with AIDS (compartment A), and this occurs at a constant rate represented by the parameter η.

3. The case detection test can determine the level of CD4+T cells in the human body. Hence, we have a transition rate τ from I_i to T_i for i = 1, 2, 3 due to successful case detection.

4. Each treated compartment will get antiretroviral treatment which will delay the progression of HIV to AIDS. We assume that successful treatment will increase the number of CD4+T cells. With a duration of treatment is ρ, we have the probability of treatment being unsuccessful given by s, q, and r for individuals in T₁, T₂, and T₃, respectively.

5. The infection rate of individuals in I_i is given by β_u while for individuals in T_i is given by β_t. As an effect of treatment, here we assume β_t < β_u.

Table 1

Table 1. Description of variables of the HIV/AIDS model in the system (3).

Based on this model description and the transmission diagram in Figure 4, the mathematical model of HIV/AIDS considering the level of CD4+T cells, antiretroviral treatment, and case detection is given by the following system of differential equations.

$\begin{array}{ll}\begin{array}{l}\frac{dS}{dt}=\Lambda -((1-ϵ\kappa )\left[{\beta }_u(I_1+I_2+I_3)+{\beta }_t(T_1+T_2+T_3)\right]+\mu )S,\\ \frac{d{I}_1}{dt}=(1-ϵ\kappa )\left[{\beta }_u(I_1+I_2+I_3)+{\beta }_t(T_1+T_2+T_3)\right]S-\\ ({\tau }_1+{\delta }_1+\mu )I_1,\\ \frac{d{I}_2}{dt}={\delta }_1{I}_1-({\tau }_2+{\delta }_2+\mu )I_2,\\ \frac{d{I}_3}{dt}={\delta }_2{I}_2-({\tau }_3+{\gamma }_u+\mu )I_3,\\ \frac{d{T}_1}{dt}={\tau }_1{I}_1+(1-q)\rho T_2-(s\rho +\mu )T_1,\\ \frac{d{T}_2}{dt}={\tau }_2{I}_2+s\rho T_1+(1-r)\rho T_3-((1-q)\rho +q\rho +\mu )T_2,\\ \frac{d{T}_3}{dt}={\tau }_3{I}_3+q\rho T_2-((1-r)\rho +r{\gamma }_t+\mu )T_3,\\ \frac{dA}{dt}=r{\gamma }_t{T}_3+{\gamma }_u{I}_3-(\mu +\eta )A,\end{array}& (3)\end{array}$

completed with a non-negative initial condition S(0) > 0, I₁(0) ⩾ 0, I₂(0) ⩾ 0, I₃(0) ⩾ 0, T₁(0) ⩾ 0, T₂(0) ⩾ 0, T₃(0) ⩾ 0, and A(0) ⩾ 0.

Figure 4

Figure 4. Transmission diagram of the HIV/AIDS model in system (3).

3.2 Preliminary analysis

With a similar approach, as we have given for Theorems 1 and 2, the positivity and boundedness criteria of the system (3) is given in the following theorems.

Theorem 7. All solutions of the HIV/AIDS model in Equation (3) with a non-negative initial condition in ${ℝ}_{+}^8$ remain positive for all time t > 0.

Theorem 8. The region

$\begin{array}{c}𝒟=\{(S,I_1,I_2,I_3,T_1,T_2,T_3,A)\in {ℝ}_{+}^8:0<S+I_1+I_2+I_3+T_1+\\ T_2+T_3+A\le \max \left(\frac{\Lambda }{\mu },N(0)\right)\}\end{array}$

is positively invariant and attracting with respect to the system (1) with a non-negative initial condition in ${ℝ}_{+}^8$.

3.3 The HIV/AIDS-free equilibrium and the control reproduction number

The HIV/AIDS-free equilibrium of the system (3) is given by

$\begin{array}{l}{ℰ}_1=(S^0,I_1^0,I_2^0,I_3^0,T_1^0,T_2^0,T_3^0,A^0)=(\frac{\Lambda }{\mu },0,0,0,0,0,0,0).\end{array}$

Using the next-generation method (40), we determine the expression of the basic reproduction number of the system (3). Similar to the previous section, we first determine the transition (Σ) and transmission (T) matrices of the infected subcompartment of the system (3), which is given by

$\begin{array}{c}T=\left[\begin{array}{cccccc}\frac{\Lambda (1-ϵ\kappa ){\beta }_u}{\mu }& \frac{\Lambda (1-ϵ\kappa ){\beta }_u}{\mu }& \frac{\Lambda (1-ϵ\kappa ){\beta }_u}{\mu }& \frac{\Lambda (1-ϵ\kappa ){\beta }_t}{\mu }& \frac{\Lambda (1-ϵ\kappa ){\beta }_t}{\mu }& \frac{\Lambda (1-ϵ\kappa ){\beta }_t}{\mu }\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right],\end{array}$

and

$\begin{array}{l}\Sigma =\left[\begin{array}{cccccc}-a_1& 0& 0& 0& 0& 0\\ {\delta }_1& -a_2& 0& 0& 0& 0\\ 0& {\delta }_2& -a_3& 0& 0& 0\\ {\tau }_1& 0& 0& -a_4& (1-q)\rho & 0\\ 0& {\tau }_2& 0& \rho s& -a_5& (1-r)\rho \\ 0& 0& {\tau }_3& 0& q\rho & -a_6\end{array}\right],\end{array}$

where a₁ = τ₁ + δ₁ + μ, a₂ = τ₂ + δ₂ + μ, , a₃ = τ₃ + γ_u + μ, a₄ = ρs + μ, a₅ = (1 − q)ρ + qρ + μ, and a₆ = (1 − r)ρ + rγ_t + μ. Since we have only the first row of T is non-zero, while the others are zero, we define

$\begin{array}{l}\text{E}=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right],\end{array}$

and calculate the control reproduction number using the formula of K = − E^tTΣ⁻¹E. Hence, the control reproduction number of the system (3) is given by

$\begin{array}{l}{ℛ}_0=\rho (\text{K})=\frac{b_1}{a_1}+\frac{b_1{\delta }_1}{a_1{a}_2}+\frac{b_1{\delta }_1{\delta }_2}{a_1{a}_2{a}_3}+b_2{K}_1+b_2{K}_2+b_2{K}_3,\end{array}$

where

$\begin{array}{l}\begin{array}{l}{b}_1=\frac{\Lambda (1-ϵ\kappa ){\beta }_u}{\mu },\\ b_2=\frac{\Lambda (1-ϵ\kappa ){\beta }_t}{\mu },\\ K_1=\frac{(-(r-1)((a_2{a}_3+{\delta }_1{\delta }_2)q-{\delta }_1{\delta }_2){\rho }^2+a_6{\delta }_1{a}_3(q-1)\rho -a_2{a}_3{a}_5{a}_6){\tau }_1}{a_1(((a_{6}s+a_4(r-1))q-a_{6}s){\rho }^2+a_4{a}_5{a}_6)a_2{a}_3},\\ K_2=\frac{{\tau }_2((-{\delta }_1{\delta }_2(r-1)a_4+a_{6}s{a}_2{a}_3)\rho +a_3{a}_4{a}_6{\delta }_1)}{a_1((q(r-1)a_4+a_{6}s(q-1)){\rho }^2+a_4{a}_5{a}_6)a_2{a}_3},\text{ and}\\ K_3=\frac{(s((a_2{a}_3+{\delta }_1{\delta }_2)q-{\delta }_1{\delta }_2){\rho }^2+q\rho a_3{a}_4{\delta }_1+a_4{a}_5{\delta }_1{\delta }_2){\tau }_3}{a_1(((a_{6}s+a_4(r-1))q-a_{6}s){\rho }^2+a_4{a}_5{a}_6)a_2{a}_3}.\end{array}\end{array}$

Following (41) , the local stability criteria of ℰ₁ is given in the following theorem.

Theorem 9. The HIV/AIDS-free equilibrium of the system (3) given by ℰ₁ is locally asymptotically stable if ℛ₀ < 1 and unstable if ℛ₀ > 1.

3.4 Data fitting

Here in this section, we fit the model (3) to the data of HIV prevalence (age 15–49 years) from Eswatini, Lesotho, Botswana, and South Africa, which represent the top four countries with the highest prevalence of HIV in the world in 2023. The prevalence data that is shown in Appendix 1 can be download from (42). Some parameters on the model in Equation (3) were estimated using the data-fitting process, while the other parameters were taken from the references. Here is the explanation of how we choose the value of these parameters.

1. The total population N, drawn from individuals aged between 15 and 64 years in each country, is sourced from The World Bank data in (43). According to this data, the populations of Eswatini, Lesotho, Botswana, and South Africa are 736,680, 1,425,560, 1,676,630, and 39,264,160, respectively.

2. The natural death rate is denoted by μ. Using data from the World Bank (44), we have ${\mu }^{\mathrm{\text{Eswatini}}}=\frac{1}{57}{\mathrm{\text{year}}}^{-1},{\mu }^{\mathrm{\text{Lesotho}}}=\frac{1}{53}{\mathrm{\text{year}}}^{-1},{\mu }^{\mathrm{\text{Botswana}}}=\frac{1}{61}{\mathrm{\text{year}}}^{-1},$ and ${\mu }^{\mathrm{\text{South Africa}}}=\frac{1}{62}{\mathrm{\text{year}}}^{-1}$.

3. Recruitment rate is denoted by (Λ). We assume that the total natural death is approximately equal to the total newborn, hence we have Λ = μN.

4. Condom efficacy is denoted by (ϵ). Based on several references (45), the efficacy of the condom usage is between 90 and 95%. Hence, we assume that ϵ = 92.5%.

5. The progression rate due to the decreases in the number of CD4 + T Cells (δ₁, δ₂). Based on (18), we choose δ₁ = 0.33 and δ₂ = 0.34.

6. Based on (27), we use γ_u = 0.6. Since γ_t < γ_u due to treatment that has been followed by T individuals, then we assumed γ_t = 0.3.

7. Since q, r, and s are proportions, we assume that ρ = 1.64, r = 0.5, q = 0.653, and s = 0.653.

8. Due to the treatment undertaken by individuals in T, then we assume that β_t < β_u.

9. We assume that individuals in T₃ are easier to be asked to follow the treatment process for HIV rather than individuals in T₁ and T₂. Hence, we assume τ₃ > τ₂ > τ₁ > 0.

10. Since κ is the proportion of people who use condoms during sexual activity, we assume κ ∈ [0, 1].

Our aim is to estimate β_u, β_t, τ₁, τ₂, τ₃, and ϵ which minimize the following error function

$\begin{array}{l}\text{Error}={\displaystyle \sum _{i=1}^{32}}{\left(\frac{I_1+I_2+I_3+T_1+T_2+T_3+A}{N}\times 100-P_i\right)}^2,\end{array}$

where P_i is the HIV prevalence data at time step−i and N is the total population in 1990 for Eswatini, Lesotho, Botswana, and South Africa as described before. The initial conditions of I_i(0), T_i(0), and A(0), for i = 1, 2, 3 are also estimated here, with $S(0)=N-\sum _{i=1}^3(I_i(0)+T_i(0)+A(0))$. The particle swarm optimization (PSO) is used to find the minimum of error function. By using the PSO algorithm as given in (46), with the number of particles given is 500, the maximum iteration is 1,000, and c₁ = c₂ = 1, we obtain the estimation of the parameters and the initial conditions as shown in Appendix 2. Figure 5 displays the weekly prevalence calculations from the estimated results compared to the HIV prevalence data for Eswatini, Lesotho, Botswana, and South Africa. The model simulation in Figure 5 shows a better agreement with the actual data. The basic reproduction numbers for each country are given as follows: 1.095, 1.682, 1.732, and 1.65 for Eswatini, Lesotho, Botswana, and South Africa, respectively. A 95% confidence interval for the fitting results using bootstrap resampling residual approach is also given in Figure 5. The lower and upper bounds of the confidence interval are computed by sorting 1,000 bootstrap samples and taking the 2.5 and 97.5% percentiles.

Figure 5

Figure 5. A comparison of HIV prevalence based on actual data and the model obtained from the estimated parameters for (A) Eswatini; (B) Lesotho; (C) Botswana; and (D) South Africa.

3.5 Sensitivity analysis of the model and the control reproduction number

The sensitivity methods can be used on infectious disease models to determine which variable or parameter is sensitive to a specific condition. Identifying the key critical parameters is an effective way to study such models more widely and accurately. Recently, we have used the sensitivity methods to identify the critical parameters for some infectious disease models. Suppose that an infectious disease model has m compartments x_i for i = 1, 2, ..., m and n parameters k_j for j = 1, 2, ..., n, then, the local sensitivities can be calculated using three different techniques: non-normalizations, half-normalizations, and full-normalizations. First, the equation of non-normalization sensitivities is given by

$\begin{array}{l}\begin{array}{c}{\mathcal{S}}_{k_j}^{x_i}=\frac{\partial x_i}{\partial k_j},\end{array}\end{array}$

where ${\mathcal{S}}_{k_j}^{x_i}$ is measured as a sensitivity coefficient of each x_i with respect to each parameter k_j. Then, the formula of half-normalization sensitivities is also defined below:

$\begin{array}{l}\begin{array}{c}{\mathcal{S}}_{k_j}^{x_i}=\left(\frac{1}{x_i}\right)\left(\frac{\partial x_i}{\partial k_j}\right).\end{array}\end{array}$

Finally, the equation of full-normalization sensitivities is defined by

$\begin{array}{l}\begin{array}{c}{\mathcal{S}}_{k_j}^{x_i}=\left(\frac{k_j}{x_i}\right)\left(\frac{\partial x_i}{\partial k_j}\right).\end{array}\end{array}$

We have used such estimated parameters and initial variables in computational simulations using MATLAB. The results given in this work provide an important step forward to understand the model dynamics more widely. This helps us to identify critical model parameters and how each model state is affected by the model parameters.

3.5.1 Model sensitivity analysis in Eswatini

In this computational simulation, results from Figure 6 are computed by using incidence data from Eswatini with the model initial populations S(0) = 854, 011, I₁(0) = 9, 529, I₂(0) = 11, I₃(0) = 7, 970, T₁(0) = 63, T₂(0) = 3, 214, T₃(0) = 8, 920, and A(0) = 1, 404, and the estimated model parameters are μ = 0.013, δ₁ = 0.33, δ₂ = 0.34, Λ = 11861.26, γ_u = 0.1, γ_t = 0.018, ρ₁₂ = 0.1462, ρ₂₁ = 0.57, ρ₂₃ = 0.1462, ρ₃₂ = 0.82, β_u = 0.9529, β_t = 0.001, τ = 0.7970, κ = 0.0063, and ϵ = 0.3214

Figure 6

Figure 6. Local sensitivity analysis with non-normalization (A, B), half-normalization (C, D), and full-normalization (E, F) techniques using the best-fit parameter for Eswatini.

When we use the non-normalization technique, it can be seen that the model states S and T₁ are very sensitive to μ, β_t, β_u, κ, and ρ₁₂, while all model states are less sensitive to the model parameters δ₂, Λ, ρ₃₂, and ϵ; see panels (A, B). Furthermore, using the half-normalization method shows that almost all model states are very sensitive to the model parameters μ, β_t, β_u, and τ, whereas they are less sensitive to the other model parameters; see panels (C, D). Interestingly, applying the full-normalization method shows that almost all model variables are sensitive to β_u and τ, while the other variables have different sensitivities to the model parameters, this is clearly seen in panels (E, F).

3.5.2 Model sensitivity analysis in Lesotho

The results from Figure 7 are computed by using incidence data from Lesotho with the model initial populations S(0) = 1, 799, 000, I₁(0) = 3, 676, I₂(0) = 5, 018, I₃(0) = 37, 272, T₁(0) = 3, 420, T₂(0) = 1, 701, T₃(0) = 1, 982, and A(0) = 1, 163, and the estimated model parameters are μ = 0.013, δ₁ = 0.33, δ₂ = 0.34, Λ = 1, 799, 000/72, γ_u = 0.1, γ_t = 0.018, ρ₁₂ = 0.1462, ρ₂₁ = 0.57, ρ₂₃ = 0.1462, ρ₃₂ = 0.82, β_u = 0.8999, β_t = 0.0001, τ = 0.7909, κ = 0.0131, and ϵ = 0.2710

Figure 7

Figure 7. Local sensitivity analysis with non-normalization (A, B), half-normalization (C), and full-normalization (D) techniques using the best-fit parameter for Lesotho.

By using the non-normalization technique, it shows that the model states S and T₁ are very sensitive to μ, β_t, β_u and τ, while the other model states are less sensitive to the model parameters; see panels (A, B). When we use the half-normalization method, it shows that almost all model states are very sensitive to the model parameters μ, β_t, β_u, and τ, whereas they are less sensitive to the other model parameters; see panel (C). Furthermore, by applying the full-normalization method shows that almost all model variables are sensitive to β_u and τ, while we can also observe that there are different levels of sensitivities between model variables and parameters; see panel (D).

3.5.3 Model sensitivity analysis in Botswana

Using the model initial states and estimated model parameters in computational simulations, results from Figure 8 are computed by using incidence data from Botswana with the model initial populations S(0) = 1, 341, 000, I₁(0) = 17, 426, I₂(0) = 3, 379, I₃(0) = 39, 459, T₁(0) = 31, 384, T₂(0) = 7, T₃(0) = 14, 015, and A(0) = 1, 243, and the estimated model parameters are μ = 0.013, δ₁ = 0.33, δ₂ = 0.34, Λ = 1341000/72, γ_u = 0.1, γ_t = 0.018, ρ₁₂ = 0.1462, ρ₂₁ = 0.57, ρ₂₃ = 0.1462, ρ₃₂ = 0.82, β_u = 0.8991, β_t = 0.0001, τ = 0.8689, κ = 0.8974, and ϵ = 0.0067.

Figure 8

Figure 8. Local sensitivity analysis with non-normalization (A, B), half-normalization (C), and full-normalization (D) techniques using the best-fit parameter for Botswana.

The results of the non-normalization technique shows that the model states S and T₁ are very sensitive to μ, β_t, β_u, τ, and ϵ, while there are different levels of sensitivities between other model states and parameters; see panels (A, B). Using the half-normalization method shows that almost all model states are very sensitive to the model parameters μ and β_t, whereas there are sensitivities to the other model parameters; see panel (C). Furthermore, applying the full-normalization method shows that almost all model variables are sensitive to β_u and τ, while we can also see that there are different levels of sensitivities between model variables and parameters; see Figure panel (D).

3.5.4 Model sensitivity analysis in South Africa

Computational results shown in Figure 9 are computed by using incidence data from South Africa with the model initial populations S(0) = 39, 880, 000, I₁(0) = 127, 900, I₂(0) = 145, 386, I₃(0) = 101, 214, T₁(0) = 104, T₂(0) = 197, 709, T₃(0) = 99, 926, and A(0) = 5, 033, and the estimated model parameters are μ = 0.013, δ₁ = 0.33, δ₂ = 0.34, Λ = 39, 880, 000/72, γ_u = 0.1, γ_t = 0.018, ρ₁₂ = 0.1462, ρ₂₁ = 0.57, ρ₂₃ = 0.1462, ρ₃₂ = 0.82, β_u = 0.8999, β_t = 0.0001, τ = 0.7945, κ = 0.6706, and ϵ = 0.0001.

Figure 9

Figure 9. Local sensitivity analysis with non-normalization (A, B), half-normalization (C), and full-normalization (D) techniques using the best-fit parameter for South Africa.

When we use the non-normalization technique, it can be seen that the model states S and T₁ are very sensitive to μ, β_t, β_u.κ, and ϵ, while the other states are less sensitive to the other model parameters; see panels (A, B). Furthermore, using the half-normalization method shows that almost all model states are very sensitive to the model parameters μ, β_t, β_u, and ϵ, whereas they are less sensitive to the other model parameters; see panel (C). Interestingly, applying the full-normalization method shows that almost all model variables are sensitive to β_u and τ, while the other variables have different sensitivities to the model parameters, this is clearly seen in panel (D).

3.6 Autonomous simulation

In this simulation, we conduct the analysis using estimated parameters for Eswatini. We divide our simulation into three scenarios to understand the impact of the infection rate, the quality of prevention (condom use), and the effectiveness of massive detection. We employ bifurcation parameters and autonomous simulations at several sample points on the bifurcation diagram.

3.6.1 Effect of infection rate β_u

To conduct the simulation, we set β_u as the bifurcation parameter, while the other parameters are the best-fit parameter for Eswatini (see Section 3.5.1). The numerical results are presented in Figure 10. It is clear to see that a larger value of β_u will increase ℛ₀ and I₁ in endemic equilibrium. Based on this dataset, we determine that ℛ₀ = 1 when ${\beta }_u=1.23\times 10^{-6}$. At a sample point P₂ (${\beta }_u=0.4\times 10^{-6}$, ℛ₀ = 0.687), we observe that the HIV/AIDS-free equilibrium point is stable, with the following values:

$S=736,680,I_1=0,I_2=0,I_3=0,T_1=0,T_2=0,T_3=0,A=0.$

On the other hand, at P₁ (${\beta }_u=3\times 10^{-6}$, we have ℛ₀ = 1.66), we observe that the HIV/AIDS-endemic equilibrium is stable, with the following values:

$\begin{array}{c}S=443342,I_1=5135,I_2=1361,I_3=142,T_1=15196,\\ T_2=23158,T_3=25493,\text{ and }A=16307.\end{array}$

In panel (A) of Figure 10, it is evident that for ${\beta }_u<1.23\times 10^{-6}$, ℛ₀ < 1, indicating the stability of the HIV/AIDS-free equilibrium. As β_u increases, ℛ₀ also increases (as seen in the cyan curve). Upon reaching the branch point (BP), the HIV/AIDS-free equilibrium becomes unstable, leading to the stable HIV/AIDS-endemic equilibrium, which grows in significance as β_u continues to increase. Panel (B) illustrates how trajectories from different initial conditions converge toward the same stable equilibrium point, namely, the HIV/AIDS-free equilibrium or the HIV/AIDS-endemic equilibrium. Figure 11 shows the dynamic of the system (3) with respect to various values of β_u. It can be seen that larger β_u will increase the number of infected individuals I_i, T_i, and A.

Figure 10

Figure 10. The bifurcation diagram of the system (3) with respect to β_u is presented in (A). In this diagram, “BP" denotes the branch point occurring at ℛ₀ = 1. The cyan, red, and blue curves represent ℛ₀, the endemic equilibrium, and the HIV/AIDS-free equilibrium, respectively. Solid and dotted curves indicate stable and unstable equilibria, respectively. (B) Depicts the trajectories of total infected, total treated, and susceptible individuals tending to HIV/AIDS free-equilibrium point for β_u at P₁ and to HIV/AIDS endemic equilibrium for β_u at P₂.

Figure 11

Figure 11. The impact of β_u on the dynamics of susceptible (up left), total infected (up right), total treated (below left), and AIDS (below right). The blue curve represents the dynamic for a low transmission rate β_u (1.488 × 10⁻⁷), while the red curve represents the dynamic for a high transmission rate β_u (1.488 × 10⁻⁶).

3.6.2 The effect of proportion of the condom use (κ)

This section is dedicated to examining the impact of the proportion of people who use condoms during sexual contact (κ) on reducing the HIV/AIDS infection rate. For this simulation, we use the same parameter values as in Section 3.6.1, with the exception of κ, which is varying. The numerical results are presented in Figure 12. With this data set, we determine that ${\mathcal{R}}_0=1$ at ϵ = 0.691. As indicated by the cyan curve in panel (A), a higher quality of condoms (larger ϵ) leads to a smaller ℛ₀. At sample point P₁, with ϵ = 0.2, we have ℛ₀ = 2.262, resulting in a stable HIV/AIDS-endemic equilibrium at

$\begin{array}{c}S=325,640,I_1=7,195,I_2=1,908,I_3=199,T_1=21,293,\\ T_2=32,450,T_3=35,722,\text{ and }A=22851.\end{array}$

On the other hand, at sample point P₂, when ϵ = 0.9, we have ℛ₀ = 0.464 which gives us a stable HIV/AIDS-free equilibrium point at

$\begin{array}{c}S=736,680,I_1=0,I_2=0,I_3=0,T_1=0,T_2=0,\\ T_3=0,\text{ and }A=0.\end{array}$

Panel (B) shows the trajectories of all solutions from different initial conditions tending toward their stable endemic equilibrium and disease-free equilibrium. Figure 13 shows the dynamic of the system (3) with respect to various values of κ. It can be seen that larger value of κ will reduce the number of infected individuals I_i, T_i, and A.

Figure 12

Figure 12. The bifurcation diagram of the system (3) with respect to κ is presented in (A). In this diagram, “BP” denotes the branch point occurring at ℛ₀ = 1. The cyan, red, and blue curves represent ℛ₀, the endemic equilibrium, and the HIV/AIDS-free equilibrium, respectively. Solid and dotted curves indicate stable and unstable equilibria, respectively. (B) Depicts the trajectories of total infected, total treated, and susceptible individuals tending to HIV/AIDS-endemic equilibrium point for κ at P₁ and to HIV/AIDS-free equilibrium for κ at P₂.

Figure 13

Figure 13. The impact of κ to the dynamics of susceptible (up left), total infected (up right), total treated (below left), and AIDS (below right). The red curve represents the dynamic when no individuals use condom, while the blue curve represents the dynamic when all people use condoms during sexual contact.

3.6.3 The effect of massive case detection from I₁ to T₁ (τ₁)

This section is dedicated to exploring the impact of the rate of case detection (τ₁) on the acceleration of treatment for individuals. In this simulation, we employ data for Eswatini, with the exception of τ₁, which serves as a freely adjustable bifurcation parameter. The numerical findings are presented in Figure 14. Within this dataset, we reveal the refined result that ${\mathcal{R}}_0=1$ when τ₁ = 10.93. As described by the cyan curve in panel (A), a more intensive case detection (characterized by a larger τ₁) corresponds to a smaller ℛ₀. At sample point P₁, τ = 7 gives ℛ₀ = 1.04 and gives us the following endemic equilibrium.

$\begin{array}{c}=708,205,I_1=68,I_2=18,I_3=2,T_1=1,662,T_2=2,343,\\ T_3=2,546,\text{ and }A=1,597.\end{array}$

On the other hand, at sample point P₂, when τ = 15, we have ℛ₀ = 0.98, which gives us a stable HIV/AIDS-free equilibrium point at

$\begin{array}{c}S=736,680,I_1=0,I_2=0,I_3=0,T_1=0,\\ T_2=0,T_3=0,\text{ and }A=0.\end{array}$

Panel (B) shows the trajectory of all solutions from different initial conditions tends toward either how stable endemic equilibrium or disease-free equilibrium. Figure 15 shows the dynamic of the system (3) with respect to various values of τ₁. It can be seen that a larger value of τ will reduce the number of infected individuals I_i, T_i, and A.

Figure 14

Figure 14. The bifurcation diagram of the system (3) with respect to τ₁ is presented in (A). In this diagram, “BP” denotes the branch point occurring at ℛ₀ = 1. The cyan, red, and blue curves represent ℛ₀, the endemic equilibrium, and the HIV/AIDS-free equilibrium, respectively. Solid and dotted curves indicate stable and unstable equilibria, respectively. (B) Depicts the trajectories of total infected, total treated, and susceptible individuals tending to HIV/AIDS-endemic equilibrium point for τ₁ at P₁ and to HIV/AIDS-free equilibrium for τ at P₂.

Figure 15

Figure 15. The impact of τ₁ on the dynamics of susceptible (up left), total infected (up right), total treated (below left), and AIDS (below right). The red curve represents when no early case detection was implemented, while the blue curve represents when it was implemented.

3.7 Minimum proportion on the use of condoms to eliminate HIV/AIDS

In this section, we conduct numerical experiments to understand how the proportion of people using condoms impacts the spread of HIV/AIDS. We use four different datasets representing four countries: Eswatini, Lesotho, Botswana, and South Africa. All parameters are consistent with those in Appendix 2, except for condom quality (ϵ). We calculate the minimum population proportion that needs to use condoms with a specific condom quality (between 70 and 100%) to achieve ℛ₀ < 1. The results are depicted in Figure 16. It is clear that better condom quality requires a smaller proportion of the population to reduce ℛ₀ to less than one. Based on the analysis results, an efficacy of 70% for condoms indicates that Eswatini could eliminate HIV/AIDS if a minimum of 91.39% of the infected population consistently uses condoms during sexual activities. Simultaneously, with the same efficacy, it is evident that Lesotho would need to mandate 100% condom usage to achieve the same goal, while Botswana cannot solely rely on condom use to eliminate HIV/AIDS. If condom efficacy improves, for instance, reaching 90%, the proportion of the infected population required to use condoms during sexual activities decreases. It would be 71% for Eswatini and up to 81% for Botswana.

Figure 16

Figure 16. A comparison of the minimum population percentage in each country required for condom usage, contingent upon the quality of the condom.

4 Extension of the HIV/AIDS model as an optimal control problem

4.1 Optimal control model

In this section, we expand our model from the system (3) to an optimal control problem by introducing time-dependent variables u₁(t) to represent the proportion of condom use (κ) and the case detection rate (τ₁, τ₂, τ₃), which is now denoted by u₂(t), u₃(t), u₄(t), respectively. Consequently, our model now reads as follows:

$\begin{array}{ll}\hfill & \frac{dS}{dt}=\Lambda -((1-ϵu_1(t))[{\beta }_u(I_1+I_2+I_3)+\hfill \\ \hfill & {\beta }_t(T_1+T_2+T_3)]+\mu )S,\hfill \\ \frac{d{I}_1}{dt}\hfill & =(1-ϵu_1(t))[{\beta }_u(I_1+I_2+I_3)+\hfill \\ \hfill & {\beta }_t(T_1+T_2+T_3)]S-(u_2(t)+{\delta }_1+\mu )I_1,\hfill \\ \frac{d{I}_2}{dt}\hfill & ={\delta }_1{I}_1-(u_3(t)+{\delta }_2+\mu )I_2,\hfill \\ \frac{d{I}_3}{dt}\hfill & ={\delta }_2{I}_2-(u_4(t)+{\gamma }_u+\mu )I_3,\hfill \\ \frac{d{T}_1}{dt}\hfill & =u_2(t)I_1+(1-q)\rho T_2-(s\rho +\mu )T_1,\hfill \\ \frac{d{T}_2}{dt}\hfill & =u_3(t)I_2+s\rho T_1+(1-r)\rho T_3-((1-q)\rho +q\rho +\mu )T_2,\hfill \\ \frac{d{T}_3}{dt}\hfill & =u_4(t)I_3+q\rho T_2-((1-r)\rho +r{\gamma }_t+\mu )T_3, \text{and}\hfill \\ \frac{dA}{dt}\hfill & =r{\gamma }_t{T}_3+{\gamma }_u{I}_3-(\mu +\eta )A,\hfill & (4)\end{array}$

The objective of this optimal control approach is to minimize the number of untreated infected individuals I₁, I₂, and I₃, as well as A, by optimizing the intervention of the proportion of condom use (u₁) and the case detection rate (u₂, u₃, and u₄). Therefore, the cost function reads as follows:

$\begin{array}{c}\mathcal{J}(u_1,u_2)={\displaystyle {\int }_0^{t_f}(}{\omega }_1{u}_1^2+{\omega }_2{u}_2^2+{\omega }_3{u}_3^2+{\omega }_4{u}_4^4+\\ {\phi }_1{I}_1+{\phi }_2{I}_2+{\phi }_3{I}_3+{\phi }_{4}A)dt,& (5)\end{array}$

where ω_i for i = 1, 2, 3, and 4 and φ_j for j = 1, 2, 3, and 4 are the positive weight parameter for each component on $\mathcal{J}$.

4.2 Characterization of the problem

We define the Hamiltonian of our problem as follows:

$\begin{array}{l}\mathscr{H}={\omega }_1{u}_1^2+{\omega }_2{u}_2^2+{\omega }_3{u}_3^2+{\omega }_4{u}_4^2+{\phi }_1{I}_1+{\phi }_2{I}_2+{\phi }_3{I}_3+{\phi }_{4}A\hfill \\ +{\lambda }_1(\Lambda -((1-ϵu_1(t))[{\beta }_u(I_1+I_2+I_3)+{\beta }_t(T_1+T_2+T_3)]\hfill \\ +\mu )S)\hfill \\ +{\lambda }_2((1-ϵu_1(t))[{\beta }_u(I_1+I_2+I_3)+{\beta }_t(T_1+T_2+T_3)]S\hfill \\ -(u_2(t)+{\delta }_1+\mu )I_1)\hfill \\ +{\lambda }_3({\delta }_1{I}_1-(u_3(t)+{\delta }_2+\mu )I_2)+{\lambda }_4({\delta }_2{I}_2-(u_4(t)+{\gamma }_u+\mu )\hfill \\ I_3)\hfill \\ +{\lambda }_5(u_2(t)I_1+(1-q)\rho T_2-(s\rho +\mu )T_1)\hfill \\ +{\lambda }_6(u_3(t)I_2+s\rho T_1+(1-r)\rho T_3-((1-q)\rho +q\rho +\mu )T_2)\hfill \\ +{\lambda }_7(u_4(t)I_3+q\rho T_2-((1-r)\rho +r{\gamma }_t+\mu )T_3)\hfill \\ +{\lambda }_8(r{\gamma }_t{T}_3+{\gamma }_u{I}_3-(\mu +\eta )A).\hfill \end{array}$

First, by taking the partial derivative of $\mathcal{H}$ with respect to each variable, the adjoint system of our problem is given as follows:

$\begin{array}{l}\frac{d{\lambda }_1}{dt}=-\frac{\partial \mathcal{H}}{\partial S}\\ =(1-ϵu_1)({\beta }_u(I_1+I_2+I_3)+{\beta }_t(T_1+T_2+T_3))\\ ({\lambda }_1-{\lambda }_2)+\mu {\lambda }_1,\\ \frac{d{\lambda }_2}{dt}=-\frac{\partial \mathcal{H}}{\partial I_1}\\ =-{\phi }_1+(1-ϵu_1){\beta }_{u}S({\lambda }_1-{\lambda }_2)+{\delta }_1({\lambda }_2-{\lambda }_3)+\\ u_2({\lambda }_2-{\lambda }_5)+\mu {\lambda }_2,\\ \frac{d{\lambda }_3}{dt}=-\frac{\partial \mathcal{H}}{\partial I_2}\\ =-{\phi }_2+(1-ϵu_1){\beta }_{u}S({\lambda }_1-{\lambda }_2)+u_3({\lambda }_3-{\lambda }_6)+\\ {\delta }_2({\lambda }_3-{\lambda }_4)+\mu {\lambda }_3,\\ \frac{d{\lambda }_4}{dt}=-\frac{\partial \mathcal{H}}{\partial I_3}\\ =-{\phi }_3+(1-ϵu_1){\beta }_{u}S({\lambda }_1-{\lambda }_2)+u_4({\lambda }_4-{\lambda }_7)+\\ {\gamma }_u({\lambda }_4-{\lambda }_8)+\mu {\lambda }_4,\\ \frac{d{\lambda }_5}{dt}=-\frac{\partial \mathcal{H}}{\partial T_1}\\ =(1-ϵu_1){\beta }_{t}S({\lambda }_1-{\lambda }_2)+s\rho ({\lambda }_5-{\lambda }_6)+\mu {\lambda }_5,\\ \frac{d{\lambda }_6}{dt}=-\frac{\partial \mathcal{H}}{\partial T_2}\\ =(1-ϵu_1){\beta }_{t}S({\lambda }_1-{\lambda }_2)+(1-q)\rho ({\lambda }_6-{\lambda }_5)+\\ q\rho ({\lambda }_6-{\lambda }_7)+\mu {\lambda }_6,\\ \frac{d{\lambda }_7}{dt}=-\frac{\partial \mathcal{H}}{\partial T_3}\\ =(1-ϵu_1){\beta }_{t}S({\lambda }_1-{\lambda }_2)+(1-r)\rho ({\lambda }_7-{\lambda }_6)+\\ r{\gamma }_t({\lambda }_7-{\lambda }_8)+\mu {\lambda }_7,\\ \frac{d{\lambda }_8}{dt}=-\frac{\partial \mathcal{H}}{\partial A}\\ =-{\phi }_4+(\mu +\eta ){\lambda }_8,& (6)\end{array}$

completed with the transversality condition λ_i(t_f) = 0 for i = 1, 2, …8. The optimality condition is taken from $\frac{\partial \mathcal{H}}{\partial u_i}=0$ for i = 1, 2, 3, and 4. Hence, taking this into account with the lower and upper bounds for u_i, we have the optimal control variables that should satisfy as follows:

$\begin{array}{ll}{u}_1^{*}=min\left\{max\left\{{\displaystyle u_1^{min}},\frac{ϵS({\lambda }_2-{\lambda }_1)({\beta }_u(I_1+I_2+I_3)+{\beta }_t(T_1+T_2+T_3))}{2{\omega }_1}\right\},{\displaystyle u_1^{max}}\right\},& (\mathrm{7a})\end{array}$

$\begin{array}{ll}{u}_2^{*}=min\left\{max\left\{{\displaystyle u_2^{min}},\frac{I_1({\lambda }_2-{\lambda }_5)}{2{\omega }_2}\right\},{\displaystyle u_2^{max}}\right\},& (\mathrm{7b})\end{array}$

$\begin{array}{ll}{u}_3^{*}=min\left\{max\left\{{\displaystyle u_3^{min}},\frac{I_2({\lambda }_3-{\lambda }_6)}{2{\omega }_3}\right\},{\displaystyle u_3^{max}}\right\},\text{ and}& (\mathrm{7c})\end{array}$

$\begin{array}{lll}{u}_4^{*}& =min\left\{max\left\{{\displaystyle u_4^{min}},\frac{I_3({\lambda }_4-{\lambda }_7)}{2{\omega }_4}\right\},{\displaystyle u_4^{max}}\right\}.& (\mathrm{7d})\end{array}$

To summarize our problem, we want to minimize the cost function given in Equation (5) subject to the state system in Equation (4) completed with its initial condition, the adjoint system in Equation (6) completed with its transversality condition, and the optimality condition in Equation (7). We use a forward-backward iterative method to solve the problem. We begin by giving an intial guess for the control variables for all t and using it to solve the state system in Equation (4) forward in time. Then, we solve the adjoint system in Equation (6) backward in time with the given transversality condition. Hence, with these results, we can update the optimal control value using the formula in Equation (7). We goback to the first step until the convergence criteria are achieved, which in our case is $|{\mathcal{J}}^{\mathrm{\text{iteration-(i + 1)}}}-{\mathcal{J}}^{\mathrm{\text{iteration-(i)}}}|<10^{-5}$.

4.3 Numerical experiments

The simulation in this section was conducted using parameter values corresponding to the best-fit values for Eswatini. Please refer to Appendix 2 for details. Furthermore, we have set the value for the weight parameter on the cost function as follows:

${\omega }_1=10^5,{\omega }_2={\omega }_3={\omega }_4=5\times 10^5,{\phi }_1={\phi }_2={\phi }_3={\phi }_4=10,$

and the initial condition given by

$\begin{array}{l}\left[S(0),I_1(0),I_2(0),I_3(0),T_1(0),T_2(0),T_3(0),A(0)\right]=\\ \left[3.62,0.1,0.02,0.002,0.358,0.559,0.625,0.4\right]\times 10^5.\end{array}$

In the following, we set three scenarios for the implementation of controls.

• Scenario 1: The implementation of condom use (u₁) only, while case detection (u₂) set to be zero. Hence u₁ ≠ 0 and u₂ = u₃ = u₄ = 0.

• Scenario 2: The implementation of case detection (u₂, u₃, u₄) only, while condom use (u₁) set to be zero. Hence u₁ = 0, u₂ ≠ 0, u₃ ≠ 0, and u₄ ≠ 0.

• Scenario 3: The combination of condom use (u₁) and case detection (u₂, u₃, and u₄). Hence u₁ ≠ 0, u₂ ≠ 0, u₃ ≠ 0, and u₄ ≠ 0.

Scenario 1

The first numerical experiment involves the use of condoms as the sole strategy to prevent the spread of HIV/AIDS. The results are presented in Figure 17 where the dynamics of model output are given in panels (A–D) while the dynamics of control are shown in panels (E–H). It is clear to see that intervention in condom usage should be given almost maximal effort from the beginning of the simulation. Since condom usage can prevent new infections, we can observe an increase in the number of susceptible individuals [panel (A)] and significant decreases in the total number of PLHIV without treatment (total of I_i) and PLHIV with AIDS illness (A). Since the intervention of case detection was not used, we can see that the number of undetected case decreases and tends to zero [see panel (C)]. The total number of infections averted using this strategy is 7.59 × 10⁶ at an optimal cost of 2.06 × 10¹¹.

Figure 17

Figure 17. The dynamic of model output for (A) susceptible, (B) total PLHIV untreated I, (C) total PLHIV receiving treatment T, (D) total PLHIV with AIDS ilness A, and (E–H) for control variables under the scenario when condom use used as a single intervention.

Scenario 2

Figure 18 illustrates the situation when the government relied solely on the implementation of case detection as a single intervention (u₂, u₃, u₄ only). The dynamics of control are depicted in panels (F–H). The intervention is applied with high intensity at the beginning of the simulation and then decreasing to an almost constant for a short period and decreases again when the final time approaches. With this strategy in place, we can observe a decrease in the number of PLHIV without and with treatment [panels (B, C)], resulting in a reduced number of individuals with AIDS [panel (D)]. The total number of infected individuals averted using this strategy is 3.74 × 10⁴ at an optimal cost of 2.52 × 10¹⁶.

Figure 18

Figure 18. The dynamic of model output for (A) susceptible, (B) total PLHIV untreated I, (C) total PLHIV receiving treatment T, (D) total PLHIV with AIDS illness A, and (E–H) for control variables under the scenario when case detection used as a single intervention.

Scenario 3

The last simulation was conducted to assess the impact of combining condom usage and case detection in reducing the spread of HIV/AIDS in Eswatini. The results are presented in Figure 19, with the dynamics of u₁ and u₂ shown in panels (E–H), respectively. We observe that the intervention involving condom use should be initiated at a high rate at the beginning of the simulation while the case detection followed the dynamic of infected individuals. As a result, we witness an initial increase in the number of susceptible individuals [panel (A)] and a significant decrease in the PLHIV without treatment and PLHIV with AIDS illness [panels (B, D), respectively]. On the other hand, since the number of infected PLHIV without treatment is already reduced due to condom use, then the number of treated PLHIV is not significantly different in the case of no control scenario. Finally, the combination of condom use and case detection leads to a significant reduction in the number of PLHIV. The total number of infections averted using this strategy is 7.55 × 10⁶ at an optimal cost of 9.81 × 10¹⁴.

Figure 19

Figure 19. The dynamic of model output for (A) Susceptible, (B) total PLHIV untreated I, (C) total PLHIV receiving treatment T, (D) total PLHIV with AIDS ilness A, and (E–H) for control variables under the scenario when condom use and case detection implemented together.

Cost-effectiveness analysis

In this section, we will calculate the most effective strategy between strategies 1,2, and 3 based on its average cost-effectiveness ratio (ACER) values. Several assumptions need to be elucidated in the cost-effectiveness calculations in this section. In computing the total cost (TC), it is assumed that the total expenses incurred as a consequence of heightened control interventions constitute the overall cost. On the other hand, total infections averted (TIA) are derived from the total number of individuals successfully prevented from infection as a consequence of control interventions. For further clarity, refer to Equations (8) and (9) below.

$\begin{array}{ll}\text{ACER}=\frac{\text{Total cost for intervention (TC)}}{\text{Total number of infections averted (TIA)}},& (8)\end{array}$

where

$\begin{array}{ll}\text{TC}\hfill & ={\displaystyle {\int }_0^{t_f}[{\omega }_1{u}_1(I_1^{‡}+I_2^{‡}+I_3^{‡}+T_1^{‡}+T_2^{‡}+T_3^{‡})+{\omega }_2{u}_2{I}_1^{‡}+}\hfill \\ \hfill & {\omega }_3{I}_2^{‡}+{\omega }_4{I}_3^{‡}]dt \text{and}\hfill \\ \text{TIA}\hfill & ={\displaystyle {\int }_0^{t_f}\left[{\displaystyle \sum _{i=1}^3\left(I_i^{†}-I_i^{‡}\right)}+{\displaystyle \sum _{i=1}^3\left(T_i^{†}-T_i^{‡}\right)}\right]}dt,\hfill & (9)\end{array}$

where symbol † and ‡ represent simulation results, without and with control, respectively. A smaller value of ACER represents the most effective strategy. The results of TC, TIA, and ACER are given in Table 2.

Table 2

Table 2. Simulation result for cost-effectiveness analysis.

Based on the calculations above, we can conclude that strategy 1, which focuses solely on condom use as the single intervention, is the most effective strategy. It is followed by strategy 3, which combines condom use and case detection. Strategy 2, which relies solely on case detection as the single intervention, is the least effective strategy. Figure 20 shows the impact of all possible scenarios on HIV prevalence in Eswatini. We can observe that the HIV prevalence between scenario 1 (condom use only) and scenario 3 (condom use and case detection) is only slightly different. This affirms that case detection is less effective in reducing HIV prevalence when the implementation of condom use is already taking place at maximum effort.

Figure 20

Figure 20. A comparison between control scenarios on the HIV prevalence in Eswatini.

5 Conclusion

In this study, we developed a mathematical model to assess the effectiveness and cost-effectiveness of various strategies aimed at controlling the spread of HIV/AIDS. We considered the impact of case detection, treatment, and condom use in our model. One key aspect of our model is its acknowledgment of the fact that treatment outcomes may not always be successful, making it a more realistic representation of the disease's dynamics. We begin our analysis for a special case model, where we do not consider the number of CD4+T cells. Our mathematical investigation on the expression of the equilibrium points and the reproduction number reveals the potential of case detection to reduce the reproduction number of the virus.

Our comprehensive analysis is extended to the complete model, where we assessed the stability of the HIV/AIDS-free equilibrium point. We found that this equilibrium is stable when the control reproduction number is less than one, indicating the feasibility of disease containment.

To validate our model, we estimate our parameter values using data from four different countries: Eswatini, Lesotho, Botswana, and South Africa. Parameter estimation was performed using incidence data from these regions, and we found that the reproduction of each country is larger than one, which indicates the endemicity of HIV/AIDS in those countries. Furthermore, sensitivity analysis sheds light on the impact of condom use and case detection on HIV spread dynamics. It allowed us to identify the most influential factors in disease control.

Finally, our study delved into optimal control strategies, considering the dynamics of infected individuals when all control variables are time-dependent. From a cost-effectiveness perspective, we found that employing condom use as a sole intervention is the most effective strategy in terms of the average cost for each averted infected individual. It is worth noting that condom use not only serves as a cost-effective approach but is also a crucial tool in preventing the transmission of HIV/AIDS. While the simulation results for case detection in this research indicate its lesser effectiveness in reducing the number of HIV prevalence, it still holds importance. Case detection, despite its limitations, plays a crucial role in assisting the government in mitigating the broad social impact of HIV in the population. In light of the complexities surrounding HIV/AIDS, a comprehensive strategy that combines both condom use and case detection could offer a more robust approach to tackling the multifaceted challenges posed by the HIV epidemic. In conclusion, our research provides valuable insights into the control of HIV/AIDS, offering a comprehensive assessment of strategies and emphasizing the importance of case detection as a highly efficient and cost-effective means of disease containment.

Beyond its potential in reducing HIV spreads, there are several issues about the implementation of condom use and case detection. In societies where discussions about sex and sexual health are taboo or stigmatized, condom use campaigns may face challenges in gaining acceptance. Cultural norms and religious beliefs can influence attitudes toward condom use. Accessibility to condoms may be limited in certain cultural contexts due to factors such as affordability, availability, and distribution channels. Addressing these barriers is crucial for the success of condom use campaigns. On the other hand, HIV-related stigma and discrimination can impede case detection efforts. Fear of being ostracized or discriminated against may discourage individuals from seeking HIV testing and treatment. Privacy concerns are paramount in HIV testing and case detection. Trust in healthcare systems is essential for successful case detection. In some communities, historical distrust or negative experiences with healthcare providers may affect the willingness to engage with HIV testing and treatment services. In summary, the effectiveness of condom use campaigns and case detection for HIV depends on their alignment with social and cultural contexts. Interventions need to be culturally sensitive, addressing barriers related to stigma, discrimination, and accessibility to effectively prevent the spread of HIV and promote early detection and treatment. In diverse social and cultural settings, collaborating with local communities, leaders, and organizations can enhance the relevance and acceptance of these interventions.

Although the research results in this article show in-depth insights into the potential use of condoms and early detection in reducing HIV prevalence in case studies across four countries, there are still some aspects that can be further developed to achieve more satisfying outcomes. One of these aspects is the fact that condom use is not only for suppressing or preventing the spread of HIV but also for other sexually transmitted infections (STIs) such as syphilis, chlamydia, herpes, and human papillomavirus (HPV). Therefore, condom use as an intervention can be applied in models that consider co-infection between two or more STIs. Only few researchers have discussed coinfection models involving HIV and other STIs, as seen in (47–51). The use of condoms as an intervention that can prevent both diseases simultaneously would be intriguing to explore in future research.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

DA: Conceptualization, Formal analysis, Investigation, Methodology, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. RD: Formal analysis, Investigation, Writing – original draft. SK: Methodology, Software, Validation, Visualization, Writing – review & editing. JW: Data curation, Investigation, Software, Validation, Visualization, Writing – review & editing. PK: Formal analysis, Funding acquisition, Investigation, Supervision, Validation, Writing – review & editing. MS: Investigation, Software, Visualization, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was funded by Universitas Indonesia through the PUTI-Q1 research grant scheme (ID:NKB-485/UN2.RST/HKP.05.00/2023).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpubh.2024.1324858/full#supplementary-material

Footnotes

1. ^AIDS by the Numbers (2023). Available online at: https://www.unaids.org/en.

2. ^The Global HIV and AIDS Epidemic (2023). Available online at: https://www.hiv.gov/hiv-basics/overview/data-and-trends/global-statistics/#::text=HIV.

References

1. Kementrian Kesehatan (Kemkes). Ayo Cari Tahu Apa Itu HIV. (2022). Available online at: https://yankes.kemkes.go.id/view_artikel/754/ayo-cari-tahu-apa-itu-hiv (accessed July 20, 2023).

Google Scholar

2. World Population Review. HIV Rates by Country 2024. (2024). Available online at: https://worldpopulationreview.com/country-rankings/hiv-rates-by-country (accessed July 20, 2023).

Google Scholar

3. Universitas Negeri Surabaya (UNESA). Hari AIDS Sedunia 2022. (2022). Available online at: https://www.unesa.ac.id/hari-aids-sedunia-2022-angka-penderita-tinggi-begini-catatan-dosen-unesa (accessed July 20, 2023).

Google Scholar

4. Raphael I, Joern RR, Forsthuber TG. Memory CD4 + T cells in immunity and autoimmune diseases. Cells. (2020) 9:531. doi: 10.3390/cells9030531

PubMed Abstract | Crossref Full Text | Google Scholar

5. Oyovwevotu SO. Mathematical modelling for assessing the impact of intervention strategies on HIV/AIDS high risk group population dynamics. Heliyon. (2021) 7:e07991. doi: 10.1016/j.heliyon.2021.e07991

PubMed Abstract | Crossref Full Text | Google Scholar

6. Li Y, Liu L. The Impact of Wolbachia on Dengue Transmission Dynamics in an SEI-SIS model. Nonlinear Analysis: Real World Applications. Elsevier BV. (2021). 103363 p.

Google Scholar

7. Oke SI, Ojo MM, Adeniyi MO, Matadi MB. Mathematical modeling of malaria disease with control strategy. Commun Math Biol Neurosci. (2020) 2020:43. doi: 10.28919/cmbn/4513

Crossref Full Text | Google Scholar

8. Handari BD, Aldila D, Tamalia E, Khoshnaw SHA, Shahzad M. Assessing the impact of medical treatment and fumigation on the superinfection of malaria: a study of sensitivity analysis. Commun Biomath Sci. (2023) 6:51–73. doi: 10.5614/cbms.2023.6.1.5

Crossref Full Text | Google Scholar

9. Das K, Murthy B, Samad SA, Biswas MHA. Mathematical transmission analysis of SEIR tuberculosis disease model. Sens Int. (2021) 2:100120. doi: 10.1016/j.sintl.2021.100120

Crossref Full Text | Google Scholar

10. Aldila D, Latifa SL, Dumbela PA. Dynamical analysis of mathematical model for Bovine Tuberculosis among human and cattle population. Commun Biomath Sci. (2019) 2:55–64. doi: 10.5614/cbms.2019.2.1.6

Crossref Full Text | Google Scholar

11. Aldila D, Chavez JP, Wijaya KP, Ganegoda NC, Simorangkir GM, Tasman H, et al. A tuberculosis epidemic model as a proxy for the assessment of the novel M72/AS01_E vaccine. Commun Nonlinear Sci Numer Simul. (2023) 20:107162. doi: 10.1016/j.cnsns.2023.107162

Crossref Full Text | Google Scholar

12. Maimunah, D A. Mathematical model for HIV spreads control program with art treatment. J Phys. (2018) 974:012035. doi: 10.1088/1742-6596/974/1/012035

Crossref Full Text | Google Scholar

13. Paul JN, Mbalawata IS, Mirau SS, Masandawa L. Mathematical modeling of vaccination as a control measure of stress to fight COVID-19 infections. Chaos Solit Fract. (2023) 166:30. doi: 10.1016/j.chaos.2022.112920

PubMed Abstract | Crossref Full Text | Google Scholar

14. Balya MA, Dewi BO, Lestari FI, Ratu G, Rosuliyana H, Windyhani T, et al. Investigating the impact of social awareness and rapid test on A COVID-19 transmission model. Commun Biomath Sci. (2021) 4:46–64. doi: 10.5614/cbms.2021.4.1.5

Crossref Full Text | Google Scholar

15. Bentout S, Chekroun A, Kuniya T. Parameter estimation and prediction for coronavirus disease outbreak 2019 (COVID-19) in Algeria. AIMS Public Health. (2020) 7:306–18. doi: 10.3934/publichealth.2020026

PubMed Abstract | Crossref Full Text | Google Scholar

16. Bentout S, Tridane A, Djilali S. Age-structured modeling of COVID-19 epidemic in the USA, UAE and Algeria. Alexandria Eng J. (2021) 60:401–11. doi: 10.1016/j.aej.2020.08.053

Crossref Full Text | Google Scholar

17. Aldila D, Nadya AF, Herdicho FF, Ndii MZ, Chukwu CW. Optimal control of pneumonia transmission model with seasonal factor: learning from Jakarta incidence data. Heliyon. (2023) 9:e18096. doi: 10.1016/j.heliyon.2023.e18096

PubMed Abstract | Crossref Full Text | Google Scholar

18. Rahman SA, Vaidya NK, Zou X. Impact of early treatment programs on HIV epidemics: an immunity-based mathematical model. Math Biosci. (2016) 280:38–49. doi: 10.1016/j.mbs.2016.07.009

PubMed Abstract | Crossref Full Text | Google Scholar

19. Mukandavire Z, Garira W, Tchuenche JM. Modelling effects of public health educational campaigns on HIV/AIDS transmission dynamics. Appl Math Model. (2009) 33:2084–95. doi: 10.1016/j.apm.2008.05.017

PubMed Abstract | Crossref Full Text | Google Scholar

20. Nyabadza F, Mukandavire Z. Modelling HIV/AIDS in the presence of an HIV testing and screening campaign. J Theor Biol. (2011) 280:167–79. doi: 10.1016/j.jtbi.2011.04.021

PubMed Abstract | Crossref Full Text | Google Scholar

21. Zhai X, Li W, Wei F, Mao X. Dynamics of an HIV/AIDS transmission model with protection awareness and fluctuations. Chaos Solit Fract. (2023) 169:113224. doi: 10.1016/j.chaos.2023.113224

Crossref Full Text | Google Scholar

22. Jamil S, Farman M, Akgul A. Qualitative and quantitative analysis of a fractal fractional HIV/AIDS model. Alexand Eng J. (2023) 76:167–77. doi: 10.1016/j.aej.2023.06.021

Crossref Full Text | Google Scholar

23. Xu C, Liu Z, Pang Y, Akgul A, Baleanu D. Dynamics of HIV-TB coinfection model using classical and Caputo piecewise operator: a dynamic approach with real data from South-East Asia, European and American regions. Chaos Solit Fract. (2022) 165:112879. doi: 10.1016/j.chaos.2022.112879

Crossref Full Text | Google Scholar

24. Pinto CMA, Carvalho ARM. New findings on the dynamics of HIV and TB coinfection models. Appl Math Comput. (2014) 242:36–46. doi: 10.1016/j.amc.2014.05.061

Crossref Full Text | Google Scholar

25. Ringa N, Diagne ML, Rwezaura H, Omame A, Tchoumi SY, Tchuenche JM, et al. and COVID-19 co-infection: a mathematical model and optimal control. Inf Med Unlocked. (2022) 31:100978. doi: 10.1016/j.imu.2022.100978

Crossref Full Text | Google Scholar

26. Omame A, Raezah AA, Diala UH, Onuoha C. The optimal strategies to be adopted in controlling the co-circulation of COVID-19, dengue and HIV: insight from a mathematical model. Axioms. (2023) 12:773. doi: 10.3390/axioms12080773

Crossref Full Text | Google Scholar

27. Garba S, Gumel A. Mathematical recipe for HIV elimination in Nigeria. J Niger Math Soc. (2010) 29:51–112.

Google Scholar

28. Greenhalgh D, Hay G. Mathematical modelling of the spread of HIV/AIDS amongst injecting drug users. Math Med Biol. (1997) 14:11–38. doi: 10.1093/imammb/14.1.11

PubMed Abstract | Crossref Full Text | Google Scholar

29. Gumel AB, McCluskey CC, Van Den Driessche P. Mathematical study of a staged-progression HIV model with imperfect vaccine. Bull Math Biol. (2006) 68:2105–28. doi: 10.1007/s11538-006-9095-7

PubMed Abstract | Crossref Full Text | Google Scholar

30. Punyacharoensin N, Edmunds WJ, De Angelis D, White RG. Mathematical models for the study of HIV spread and control amongst men who have sex with men. Eur J Epidemiol. (2011) 26:695–709. doi: 10.1007/s10654-011-9614-1

PubMed Abstract | Crossref Full Text | Google Scholar

31. Ofosuhene OA, Prince PO, Noor AI, Aline C. Analysing the impact of migration on HIV/AIDS cases using epidemiological modelling to guide policy makers. Infect Math Model. (2022) 7:252–61. doi: 10.1016/j.idm.2022.01.002

PubMed Abstract | Crossref Full Text | Google Scholar

32. Fatmawati, Khan MA, Odinsyah HP. Fractional model of HIV transmission with awareness effect. Chaos Solit Fract. (2020) 138:109967. doi: 10.1016/j.chaos.2020.109967

Crossref Full Text | Google Scholar

33. Yusuf A, Mustapha UT, Sulaiman TA, Hincal E, Bayram M. Modeling the effect of horizontal and vertical transmissions of HIV infection with Caputo fractional derivative. Chaos Solit Fract. (2021) 145:110794. doi: 10.1016/j.chaos.2021.110794

Crossref Full Text | Google Scholar

34. Tabassum MF, Saeed M, Akgul A, Farman M, Chaudhry NA. Treatment of HIV/AIDS epidemic model with vertical transmission by using evolutionary Padé-approximation. Chaos Solit Fract. (2020) 2020:1099686. doi: 10.1016/j.chaos.2020.109686

Crossref Full Text | Google Scholar

35. Castillo-Chavez C, Feng Z, Huang W. On the computation of R0 and its role on global stability. Math Approach Emerg Reemerg Infect Dis. (2002) 1:229–50. doi: 10.1007/978-1-4757-3667-0_13

Crossref Full Text | Google Scholar

36. Aldila D, Angelina M. Optimal control problem and backward bifurcation on malaria transmission with vector bias. Heliyon. (2021) 7:e06824. doi: 10.1016/j.heliyon.2021.e06824

PubMed Abstract | Crossref Full Text | Google Scholar

37. Aldila D. Optimal control for dengue eradication program under the media awareness effect. Int J Nonlinear Sci Numer Simul. (2023) 24:95–122. doi: 10.1515/ijnsns-2020-0142

Crossref Full Text | Google Scholar

38. Aldila D, Puspadani CA, Rusin R. Mathematical analysis of the impact of community ignorance on the population dynamics of dengue. Front Appl Math Stat. (2023) 9:1094971. doi: 10.3389/fams.2023.1094971

Crossref Full Text | Google Scholar

39. Chitnis N, Hyman JM, Cushing JM. Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bull Math Biol. (2008) 70:1272. doi: 10.1007/s11538-008-9299-0

PubMed Abstract | Crossref Full Text | Google Scholar

40. Diekmann O, Heesterbeek J, Roberts MG. The construction of next-generation matrices for compartmental epidemic models. J R Soc Interface. (2010) 7:873–85. doi: 10.1098/rsif.2009.0386

PubMed Abstract | Crossref Full Text | Google Scholar

41. Van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. (2002) 180:29–48. doi: 10.1016/S0025-5564(02)00108-6

PubMed Abstract | Crossref Full Text | Google Scholar

42. The World Bank. Prevalence of HIV. (2023). Available online at: https://data.worldbank.org/indicator/SH.DYN.AIDS.ZS (accessed July 20, 2023).

Google Scholar

43. The World Bank. Population Ages 15-64, Total. (2022). Available online at: https://data.worldbank.org/indicator/SP.POP.1564.TO (accessed January 13, 2024).

Google Scholar

44. The World Bank. Life Expectancy at Birth, Total (Years). (2022). Available online at: https://data.worldbank.org/indicator/SP.DYN.LE00.IN (accessed January 13, 2024).

Google Scholar

45. Pinkerton SD, Abramson PR. Effectiveness of condoms in preventing HIV transmission. Soc Sci Med. (1997) 44:1303–12. doi: 10.1016/S0277-9536(96)00258-4

PubMed Abstract | Crossref Full Text | Google Scholar

46. Puspita JW, Fakhruddin M, Fahlena H, Rohim F, Sutimin. On the reproduction ratio of dengue incidence in Semarang, Indonesia 2015-2018. Commun Biomath Sci. (2019) 2:118–26. doi: 10.5614/cbms.2019.2.2.5

Crossref Full Text | Google Scholar

47. Mahiane SG, Nguema EP, Pretorius C, Auvert B. Mathematical models for coinfection by two sexually transmitted agents: the human immunodeficiency virus and herpes simplex virus type 2 case. J R Stat Soc Ser C. (2010) 59:547–72. doi: 10.1111/j.1467-9876.2010.00719.x

Crossref Full Text | Google Scholar

48. Nwankwo A, Okuonghae D. Mathematical analysis of the transmission dynamics of HIV syphilis co-infection in the presence of treatment for syphilis. Bull Math Biol. (2018) 80:437–92. doi: 10.1007/s11538-017-0384-0

PubMed Abstract | Crossref Full Text | Google Scholar

49. David JF, Lima VD, Zhu J, Brauer F. A co-interaction model of HIV and syphilis infection among gay, bisexual and other men who have sex with men. Infect Dis Model. (2020) 5:855–70. doi: 10.1016/j.idm.2020.10.008

PubMed Abstract | Crossref Full Text | Google Scholar

50. Teklu SW, Terefe BB. COVID-19 and syphilis co-dynamic analysis using mathematical modeling approach. Front Appl Math Stat (2023) 8:1–13. doi: 10.3389/fams.2022.1101029

Crossref Full Text | Google Scholar

51. Wang C, Gao S, Li X, Martcheva M. Modeling syphilis and HIV coinfection: a case study in the USA. Bull Math Biol. (2023) 85:20–32. doi: 10.1007/s11538-023-01123-w

PubMed Abstract | Crossref Full Text | Google Scholar