HIV/AIDS Mathematical Model of Triangle Transmission

Abstract

In this paper, a mathematical analysis of the HIV/AIDS deterministic model studied in the paper called Mathematical Model of HIV/AIDS Considering Sexual Preferences Under Antiretroviral Therapy, a case study in the previous works preformed by Espitia is performed. The objective is to gain insight into the qualitative dynamics of the model determining the conditions for the persistence or effective control of the disease in the community through the study of basic properties such as positiveness and boundedness; the calculus of the basic reproduction number; stationary points such as disease-free equilibrium (DFE), boundary equilibrium (BE) and endemic equilibrium (EE); and the local stability (LAS) of disease-free equilibrium. The findings allow us to conclude that the best way to reduce contagion and consequently reach a DFE is thought to be the reduction in the rate of homosexual partners, as they are the most affected population by the virus and are therefore the most likely to become infected and spread it. Increasing the departure rate of infected individuals leads to a decrease in untreated infected heterosexual men and untreated infected women.

1. Introduction

Epidemiological evidence shows that HIV is transmitted only through the exchange of body fluids such as blood, semen, vaginal or anal secretions, and breast milk. As a result, the highly common means of transmission are: unprotected sex, from mother to child during pregnancy, childbirth or breast feeding, injecting drugs with a needle that has come into contact with infected blood, and infected blood donation or organ transplant [1]. There are many myths and misconceptions about how a person can get HIV. It is not transmitted through body fluids such as sweat, tears, or saliva, touching someone who has HIV, mosquito bites, or other transmission methods.

The sexual transmission of HIV is usually considered to be carried out by heterosexual or homosexual men through anal intercourse. Transmission between two women is almost null; however, this form is possible by sharing toys such as sexual vibrators [2,3]. Female homosexual contact has not been demonstrated to pose appreciable HIV transmission risk, and such transmission appears to be rare [4,5]. According to communication with the HIV/AIDS infectious disease specialist Dr. Alexandre Naime Barbosa, the sexual transmission between men can occur through three mechanisms: exclusive homosexual transmission, exclusive heterosexual transmission, or bisexual transmission, while in women, the transmission is almost always heterosexual. The Center for Disease Control and Prevention estimates that HIV rates in men who have sex with men (MSM) are higher than the rates in heterosexual contacts. In part, these differences reflect the fact that an individual MSM can engage in both insertive and receptive sexual roles (versatility), while exclusively heterosexual men and women each engage in only one of these roles [6,7].

When discussing transmission, the term “Discordant Couples” will be used to represent a couple in which one partner has a sexually transmitted disease while the other partner does not. If two participants are infected, the transmission could imply co-infection, which is not the objective in this investigation. However, if the two participants are susceptible, then there is no contagion.

The risk of acquiring HIV is 22 times higher among men who have sex with men (MSM), 22 times higher among individuals who are injectable drug users (IDU) and share needles, 21 times higher for sex workers, and 12 times higher for transgender people compared to the risk of transmission in heterosexual contact [8,9]. One form of measuring how transmissible a disease is the “Basic Reproduction Number”, which describes secondary infections from a first infection; this number depends on the contagion’s form. For example, for HIV/AIDS transmission, the basic reproduction number is 4 in the homosexual population in the United Kingdom, whereas the basic number is 11 for female prostitutes in Kenya [10]. As a result of the variation in these statistics, we consider homosexual transmission to be greater than heterosexual transmission.

In the triangle transmission model, it is assumed that the only way to transmit the HIV virus is through sexual intercourse, and it is commonly considered that the contagion form takes into account heterosexuals and homosexuals in the dynamic of infection. However, can the population be split into heterosexuals and homosexuals and thus the group of bisexuals be ignored? Moreover, what is the contribution of these group in the transmission of HIV? To try to answer these questions, we propose a different mathematical model considering HIV-infected bisexuals under ART. Several articles have also focused on the whole population of constant size when considering force of infection, although some studies such as [11,12] have stressed the importance of variable population size in epidemic dynamics. All these assumptions, such as sexual preference and variable population in force of infection, are considered in our model.

With regard to sexual contact between homosexual men, heterosexual men and women [13,14] say: “There exist individuals that change their sexual behavior depending on the situation or at different stages in their life. A possibly common and transient example of situational sexuality is the person who self-identifies as heterosexual, but will sexually interact with a member of the same sex when lacking other opportunities. Less transient but also possibly common, a person who self-identifies as gay or lesbian (either at the time, or later) may sexually interact with a member of the opposite sex if a same-sex relationship seems unfeasible”. Thus, in our model, we consider bisexual contact.

2. Materials and Methods

The epidemiological model under consideration was studied in [15]. The model contains three population groups: the first being men with homosexual preference in men, the second being men with heterosexual preference, and the third one for women who may be homosexual or heterosexual but engage in sexual relations with homosexual or heterosexual men. We supposed that, eventually, the homosexual men had sexual contact with women and that the heterosexual men had sexual relation with homosexual men. Consequently, we consider bisexual behavior among these groups because the transmission from homosexuals to heterosexual men or women goes through the bisexuals. Female homosexual transmission is not considered in the dynamic of infection. For more information, see references [2,3].

The total population $N\left(t\right)$ is divided into eight classes; $S_h\left(t\right)$ represents susceptible homosexual men, $I_h\left(t\right)$ untreated infected homosexual men, $S_w\left(t\right)$ susceptible women, $I_w\left(t\right)$ untreated infected women, $S_m\left(t\right)$ susceptible heterosexual men, $I_m\left(t\right)$ untreated infected heterosexual men, $T\left(t\right)$ treated individuals on ART, and $A\left(t\right)$ individuals living with AIDS.

Figure 1 represents the transmission dynamics between the three studied sexual preferences. Each vertex of the triangle represents one population, and the sides of the triangle denote the different forms of transmission between the populations involved. To begin, the exclusive transmission among homosexual men is illustrated by the upper circular dotted arrow labeled as ${\lambda }_h$. Then, the transmission between homosexual and heterosexual men and the transmission between homosexual men and women are represented by dashed lines identified as ${\lambda }_{hm}$ and ${\lambda }_{hw}$, respectively. Finally, heterosexual transmission between men and women is a continuous line represented by ${\lambda }_{m,w}$. The direction of the arrows represents the sense of the analyzed contagion; nonetheless, contagions can biologically occur in all directions. Consequently, the two following hypotheses are assumed: the only form of contagion among homosexuals is among themselves, and heterosexual people become infected due to the contact with homosexual men or heterosexual partners of the opposite sex. Thus, dashed lines have only one direction, while the continuous line between heterosexual men and women has two directions. The following assumed hypotheses in the model were evaluated by HIV/AIDS specialist Dr. Alexandre Naime Barbosa from Stadual University of Sao Paulo, UNESP, Botocatu, Brazil.

Assumed Hypotheses in the Model.

H1 

Constant recruitment in all susceptible classes is assumed.

H2 

Sexual transmission in discordant couples is considered.

H3 

Homosexual individuals become infected among themselves. HIV transmission in the susceptible female population happens through sexual relations with infected heterosexual men or with infected homosexual men. Susceptible heterosexual men can become infected by infected women or infected homosexual men.

H4 

There is no gender differentiation in either sexual preference in treated individuals or individuals living with AIDS.

H5 

Individuals living with AIDS could be treated or untreated, noting that an individual that developed AIDS during a hospital treatment will be diagnosed and enrolled in ART.

H6 

It is considered both natural mortality in all classes and induced mortality in individuals living with AIDS.

Parameters in the Model

The constant recruitment in all susceptible classes is denoted by $\mathsf{\Psi }$. The male proportion is labeled by $\theta$, $0\le \theta \le 1$. The heterosexual proportion is represented by $\gamma$, $0\le \gamma \le 1$. The proportion of initially treated individuals is p, $0\le p\le 1$; consequently, $(1-p)$ denotes the proportion of untreated individuals. Natural mortality rate is symbolized by $\mu$. Induced mortality rate in individuals living with AIDS is d. AIDS development rate in treated individuals is $\delta$. Departure rate of infected individuals is $\alpha$. Subscripts $s,h,hw,hm$ mean sexual contact between heterosexual men and women, among homosexual men, between homosexual men and women, and, finally, between homosexual men and heterosexual men, respectively; thus, ${\beta }_{s,h,hw,hm}$ represents the probability of transmission and $c_{s,h,hw,hm}$ mean rate of sexual partners in the aforementioned contacts. $B_h=c_h{\beta }_h,B_s=c_s{\beta }_s,B_{hm}=c_{hm}{\beta }_{hm},B_{hw}=c_{hw}{\beta }_{hw}$ rates will be considered for parameter simplification. All parameters are non-negatives and are listed in

Table 1. Description of Parameters. Adapted from [15].

Parameter	Description
$\mathsf{\Psi }$	Constant Recruitment
$\theta$	Male Proportion
$\gamma$	Heterosexual Proportion
p	Proportion of Initially Treated Individuals
$\mu$	Natural Mortality Rate
d	Induced Disease Mortality Rate
$\delta$	AIDS Development Rate in Treated Individuals
$\alpha$	Departure Rate of infected individuals
${\beta }_{s,h,hw,hm}$	Sexual Transmission Probability
$c_{s,h,hw,hm}$	Sexual Partners Rate

.

Initially treated individuals and individuals living with AIDS receiving ART are disregarded from the transmission because their viral load is negligible. In Figure 1, we assume that susceptible homosexual men only become infected by infected homosexual men, and susceptible women (or men) become infected by infected men (or women) or infected homosexual men. This means that susceptible homosexual men select their partner randomly from the infected homosexual population, while women or men select their partners randomly from the infected heterosexual or infected homosexual population [15].

The force of infection or disease incidence function measures the susceptible person’s risk of becoming infected. In some epidemic models, this function is assumed to be bilinear in both the infected individuals and the susceptible individuals. In addition, a bilinear force of infection or mass action law incidence may not yield appropriate results for several reasons. In particular, this force of infection does not permit one to consider the difference among infected individuals. Thus, we decided that since this function represents the contact between an infected person and a susceptible one, the denominator would have to only be formed by susceptible individuals and those who can transmit the disease. We excluded both treated individuals and people living with AIDS under ART since their viral charge is negligible; in addition, people living with AIDS are too sick, and their sexual life can be considered as almost null. Therefore, the following infection forces by sexual contact are:

$$\begin{array}{cccc}\hfill {\lambda }_h& =B_h\frac{I_h}{S_h+I_h}\hfill & \hfill & \mathrm{Exclusive}\mathrm{Homosexual}\mathrm{Contact},\hfill \\ \hfill {\lambda }_{hw}& =B_{hw}\frac{I_h}{S_m+S_h+I_m+I_h}\hfill & \hfill & \mathrm{Contact}\mathrm{between}\mathrm{Homosexual}\mathrm{Men}\mathrm{and}\mathrm{Women},\hfill \\ \hfill {\lambda }_{hm}& =B_{hm}\frac{I_h}{S_m+S_h+I_m+I_h}\hfill & \hfill & \mathrm{Contact}\mathrm{between}\mathrm{Homosexual}\mathrm{Men}\mathrm{and}\mathrm{Heterosexual}\mathrm{Men},\hfill \\ \hfill {\lambda }_{m,w}& =B_s\frac{I_{m,w}}{S_{m,w}+I_{m,w}}\hfill & \hfill & \mathrm{Heterosexual}\mathrm{Contact}.\hfill \end{array}$$

It is important to note that in exclusive homosexual contact, the fraction denotes untreated infected homosexual men among susceptible and untreated infected homosexual men. However, in the contact between homosexual men and women, the fraction denotes untreated infected homosexual men among susceptible heterosexual men and untreated homosexual men because this contact is considered bisexual behavior. The same reasoning should be applied to the contact between homosexual men and heterosexual men. For heterosexual contact, the fraction denotes untreated infected heterosexual men (women) among susceptible heterosexual men (women) and untreated infected heterosexual men (women). The compartmental model is presented in Figure 2. The dynamic is governed by the system of nonlinear ordinary differential Equations (1)–(8), where a dot represents differentiation with respect to t.

$$\begin{array}{cc}\hfill \dot{S_h}& =\mathsf{\Psi }\theta (1-\gamma )-B_h\frac{I_h}{S_h+I_h}{S}_h-\mu S_h,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \dot{I_h}& =B_h\frac{I_h}{S_h+I_h}{S}_h-(\alpha +\mu )I_h,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \dot{S_w}& =\mathsf{\Psi }(1-\theta )-B_s\frac{I_m}{S_m+I_m}{S}_w-B_{hw}\frac{I_h}{S_m+S_h+I_m+I_h}{S}_w-\mu S_w,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \dot{I_w}& =B_s\frac{I_m}{S_m+I_m}{S}_w+B_{hw}\frac{I_h}{S_m+S_h+I_m+I_h}{S}_w-(\alpha +\mu )I_w,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \dot{S_m}& =\mathsf{\Psi }\theta \gamma -B_s\frac{I_w}{S_w+I_w}{S}_m-B_{hm}\frac{I_h}{S_m+S_h+I_m+I_h}{S}_m-\mu S_m,\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \dot{I_m}& =B_s\frac{I_w}{S_w+I_w}{S}_m+B_{hm}\frac{I_h}{S_m+S_h+I_m+I_h}{S}_m-(\alpha +\mu )I_m,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \dot{T}& =\alpha p(I_h+I_w+I_m)-(\delta +\mu )T,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \dot{A}& =\alpha (1-p)(I_h+I_w+I_m)+\delta T-(d+\mu )A.\hfill \end{array}$$

(8)

With initial conditions

$$\begin{array}{cc}\hfill S_h\left(0\right)& >0,I_h\left(0\right)\ge 0S_w\left(0\right)>0,I_w\left(0\right)\ge 0,S_m\left(0\right)>0,\hfill \\ \hfill I_m\left(0\right)& \ge 0,T\left(0\right)\ge 0,A\left(0\right)\ge 0.\hfill \end{array}$$

Explanation of Equations

Susceptible individuals such as homosexual men, women, and heterosexual men $S_h\left(t\right),S_w\left(t\right)$, and $S_m\left(t\right)$, grow in number with recruitment $\mathsf{\Psi }\theta (1-\gamma ),\mathsf{\Psi }(1-\theta )$, and $\mathsf{\Psi }\theta \gamma$, respectively, where $\mathsf{\Psi }$ is a constant recruitment, $\theta$ is the male proportion, and $\gamma$ is the heterosexual proportion; these susceptible populations decrease due to contagion with the virus in contact rates ${\lambda }_h,{\lambda }_m$, and ${\lambda }_w$, respectively. Women and heterosexual men additionally acquire the virus with rate ${\lambda }_{hw}$ and ${\lambda }_{hm}$. Finally, they can die from natural causes with rate $\mu$.

The number of infected individuals such as homosexual men, women, and heterosexual men, $I_h\left(t\right),$$I_w\left(t\right)$, and $I_m\left(t\right)$, grows with the rates of infection ${\lambda }_h,{\lambda }_m$, and ${\lambda }_w$. However, women and heterosexual men grow with rates ${\lambda }_{hw}$ and ${\lambda }_{hm}$, respectively. This infected population reduces because its individuals become treated or as a result of people living with AIDS in rates $\alpha$ and $\alpha (1-p)$, respectively. Finally, they die from natural causes with rate $\mu$.

The number of treated individuals, $T\left(t\right)$, grows because infected ones enroll in ART, develop AIDS with a rate $\delta$, or die from natural causes with rate $\mu$.

The number of individuals living with AIDS, $A\left(t\right)$, grows due of the entrance of infected people with or without treatment whom develop AIDS; they die from natural causes with rate $\mu$ and from induced disease death with rate d.

The correspondent mathematical analysis of this ordinary differential equations system is developed as follows.

2.1. Positiveness and Boundedness

Theorem 1.

Let the initial conditions be $S_h\left(0\right)>0,I_h\left(0\right)\ge 0,S_w\left(0\right)\ge 0,I_w\left(0\right)\ge 0,S_m\left(0\right)\ge 0,I_m\left(0\right)\ge 0,T\left(0\right)\ge 0,A\left(0\right)\ge 0$. Then, the solutions $S_h\left(t\right),I_h\left(t\right),S_w\left(t\right),$ $I_w\left(t\right),S_m\left(t\right),I_m\left(t\right),T\left(t\right),A\left(t\right)$ of the system (1) to (8) will be positive for all time $t>0$.

Proof. 

Let $t_1=sup\{t>0:S_h\left(t\right)>0,I_h\left(t\right)>0,S_w\left(t\right)>0,I_w\left(t\right)>0,S_m\left(t\right)>0,I_m\left(t\right)>0,T\left(t\right)>0,A\left(t\right)>0\}$. From the first Equation (1), we have

$${\displaystyle \frac{d{S}_h}{dt}}\left(t\right)=\mathsf{\Psi }\theta (1-\gamma )-B_h\frac{I_h\left(t\right)}{S_h\left(t\right)+I_h\left(t\right)}{S}_h\left(t\right)-\mu S_h\left(t\right)=\mathsf{\Psi }\theta (1-\gamma )-({\lambda }_h\left(t\right)+\mu )S_h\left(t\right),$$

which can be re-written as:

$$\begin{array}{cc}\hfill \frac{d}{dt}(S_h\left(t\right)exp[\mu t+{\int }_0^t{\lambda }_h\left(\tau \right)d\tau \left]\right)& =\mathsf{\Psi }\theta (1-\gamma )exp[\mu t+{\int }_0^t{\lambda }_h\left(\tau \right)d\tau ]\hfill \\ \hfill S_h\left(t_1\right)exp[\mu t_1+{\int }_0^{t_1}{\lambda }_h\left(\tau \right)d\tau ]& -S_h\left(0\right)=\mathsf{\Psi }\theta (1-\gamma ){\int }_0^{t_1}exp[\mu y+{\int }_0^y{\lambda }_h\left(\tau \right)d\tau ]dy\hfill \end{array}$$

$$\begin{array}{cc}\hfill S_h\left(t_1\right)& =S_h\left(0\right)exp[-\mu t_1-{\int }_0^{t_1}{\lambda }_h\left(\tau \right)d\tau ]\hfill \\ \hfill & +exp[-\mu t_1-{\int }_0^{t_1}{\lambda }_h\left(\tau \right)d\tau ]\mathsf{\Psi }\theta (1-\gamma ){\int }_0^{t_1}exp[\mu y+{\int }_0^y{\lambda }_h\left(\tau \right)d\tau ]dy\ge 0.\hfill \end{array}$$

Similarly, it can be shown that $I_h\left(t\right),S_w\left(t\right),I_w\left(t\right),S_m\left(t\right),I_m\left(t\right),T\left(t\right),A\left(t\right)$ are non-negatives for all time $t>0$. In this way, all solutions of the system remain positive for all non-negative initial conditions. □

Theorem 2.

All the solutions of the system (1) to (8) are uniformly bounded. It means any trajectory that starts in ${\mathbb{R}}_8^{+}$ remains in ${\mathbb{R}}_8^{+}$ for all time $t\ge 0$.

Proof. 

Adding all eight equations from (1) to (8) gives:

$$\begin{array}{cc}\hfill \frac{dN}{dt}& =\mathsf{\Psi }-\mu N-dA\hfill \\ \hfill & \le \mathsf{\Psi }-\mu N.\hfill \end{array}$$

Solving the differential in-equation, we have:

$$N\left(t\right)\le (N\left(0\right)-\frac{\mathsf{\Psi }}{\mu })exp\left(-\mu t\right)+\frac{\mathsf{\Psi }}{\mu }.$$

(9)

Therefore, all solutions of the system will enter into the region:

$${\mathsf{\Omega }}_{III}=\left\{(S_h\left(t\right),I_h\left(t\right),S_w\left(t\right),I_w\left(t\right),S_m\left(t\right),I_m\left(t\right),T\left(t\right),A\left(t\right))\in {\mathbb{R}}_8^{+}:N\left(t\right)\le \frac{\mathsf{\Psi }}{\mu }\right\}.$$

(10)

In Equation (9), if $N\left(0\right)\le \frac{\mathsf{\Psi }}{\mu }$, then $N\left(t\right)\le \frac{\mathsf{\Psi }}{\mu }$; if $N\left(0\right)\ge \frac{\mathsf{\Psi }}{\mu }$ then either the solution enters in ${\mathsf{\Omega }}_{III}$ in finite time or $N\left(t\right)$ approaches $\frac{\mathsf{\Psi }}{\mu }$ asymptotically. Therefore, ${\mathsf{\Omega }}_{III}$ attracts all solutions in ${\mathbb{R}}_{+}^8$. □

The previous theorems allow us to conclude that the region ${\mathsf{\Omega }}_{III}$ is a positively invariant set.

2.2. Basic Reproduction Number

The basic reproduction number, $R_0$, determines the ability of the virus to develop and persist in the population. It is the average number of individuals that a single infected individual can infect during their life time when introduced into a wholly susceptible population. If $R_0<1$, then, on average, a few infected individuals brought into a fully susceptible population will not be able to replace themselves and the disease will not spread. If $R_0>1$, then the number of infected individuals will increase with each generation and the disease will spread.

In this research, we use the next generation matrix method as presented in [16]. This method is as follows:

Let $x={(x_1,x_2,\dots ,x_n)}^T$ be the number of individuals in each compartment, where the first $m<n$ compartments contain infected individuals. Consider these equations written in the form:

$$\dot{x_i}=f_i\left(x\right)={\mathcal{F}}_i\left(x\right)-{\mathcal{V}}_i\left(x\right),\mathrm{for}i=1,\dots ,m.$$

(11)

In this splitting, ${\mathcal{F}}_i\left(x\right)$ is the rate of appearance of new infections in compartment i and ${\mathcal{V}}_i\left(x\right)={\mathcal{V}}_i^{-}\left(x\right)-{\mathcal{V}}_i^{+}\left(x\right)$, where ${\mathcal{V}}_i^{+}\left(x\right)$ is the rate of transfer of individuals into compartment i by all others, and ${\mathcal{V}}_i^{-}\left(x\right)$ is the rate of transfer of individuals out of the i compartment.

Note that ${\mathcal{F}}_i\left(x\right)$ includes only infections that are newly arising, but does not include terms which describe the transfer of infectious individuals from one compartment to another. Let $X_s=\{x\ge 0|x_i=0,i=1,\dots ,m\}$ be the DFE. Assume that ${\mathcal{F}}_i$ and ${\mathcal{V}}_i$ satisfy the following axioms outlined by [16]:

$\left(A_1\right)$

If $x\ge 0$, then ${\mathcal{F}}_i,{\mathcal{V}}_i^{+},{\mathcal{V}}_i^{-}\ge 0$ for $i=1,\dots ,m.$

$\left(A_2\right)$

If $x_i=0$, then ${\mathcal{V}}_i^{-}=0$. In particular, if $x\in X_s$, then ${\mathcal{V}}_i^{-}=0$ for $i=1,\dots ,m.$

$\left(A_3\right)$

${\mathcal{F}}_i=0$ if $i>m$;

$\left(A_4\right)$

If $x\in X_s$, then ${\mathcal{F}}_i\left(x\right)=0$ and ${\mathcal{V}}_i^{+}=0$ for $i=1,\dots ,m.$

$\left(A_5\right)$

All eigenvalues of $Df\left(x_0\right)$ have negative real parts, where $Df\left(x_0\right)$ is the Jacobian matrix evaluated at the disease free equilibrium $x_0$.

Theorem 3

(Exposed in [16]). If $x_0$ is the disease free equilibrium (DFE) and $f_i\left(x\right)$ satisfies $\left(A_1\right)-\left(A_5\right)$, then the derivatives $D\mathcal{F}\left(x_0\right)$ and $D\mathcal{V}\left(x_0\right)$ are partitioned as:

$$D\mathcal{F}\left(x_0\right)=\left(\begin{array}{cc}F& 0\\ 0& 0\end{array}\right),D\mathcal{V}\left(x_0\right)=\left(\begin{array}{cc}V& 0\\ J_3& J_4\end{array}\right).$$

where F and V are the $m\times m$ matrices defined by:

$$F=\left[\frac{\partial {\mathcal{F}}_i}{\partial x_j}\left(x_0\right)\right],V=\left[\frac{\partial {\mathcal{V}}_i}{\partial x_j}\left(x_0\right)\right]\mathit{with}1\le i,j\le m.$$

Furthermore, F is non-negative, V is a non-singular $M-$matrix, and all eigenvalues of $J_4$ have a positive real part.

According to [17], $F{V}^{-1}$ is called the next generation matrix for model (11), and the spectral radius (dominant eigenvalue) is the basic reproduction number:

$$R_0=\rho \left(F{V}^{-1}\right).$$

(12)

Theorem 4

(Exposed in [16]). Consider the disease transmission model given by (11) with $f\left(x\right)$ satisfying conditions $\left(A_1\right)$ to $\left(A_5\right)$. If $x_0$ is a DFE of the model, then $x_0$ is locally asymptotically stable if $R_0<1$, but unstable if $R_0>1$, where $R_0$ is defined by Equation (12).

The basic reproduction number is defined as the spectral radius of the matrix $F{V}^{-1}$ and denoted by:

$$R_0=max\left\{\frac{B_s}{\alpha +\mu },\frac{B_h}{\alpha +\mu }\right\}=max\{R_0^{het},R_0^{hom}\}.$$

(13)

Details are presented in Appendix A.

2.3. Stationary Points

To calculate stationary points, we solve the associated homogeneous system (14)–(21), state variables with a star (*) superscript will be assumed to be an equilibrium value:

$$\begin{array}{cc}\hfill 0& ={\mathsf{\Lambda }}_h-B_h\frac{I_h^{*}}{S_h^{*}+I_h^{*}}{S}_h^{*}-\mu S_h^{*},\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill 0& =B_h\frac{I_h^{*}}{S_h^{*}+I_h^{*}}{S}_h^{*}-(\alpha +\mu )I_h^{*},\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill 0& ={\mathsf{\Lambda }}_w-B_s\frac{I_m^{*}}{S_m^{*}+I_m^{*}}{S}_w^{*}-B_{hw}\frac{I_h^{*}}{S_m^{*}+S_h^{*}+I_m^{*}+I_h^{*}}{S}_w^{*}-\mu S_w^{*},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill 0& =B_s\frac{I_m^{*}}{S_m^{*}+I_m^{*}}{S}_w^{*}+B_{hw}\frac{I_h^{*}}{S_m^{*}+S_h^{*}+I_m^{*}+I_h^{*}}{S}_w^{*}-(\alpha +\mu )I_w^{*},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill 0& ={\mathsf{\Lambda }}_m-B_s\frac{I_w^{*}}{S_w^{*}+I_w^{*}}{S}_m^{*}-B_{hm}\frac{I_h^{*}}{S_m^{*}+S_h^{*}+I_m^{*}+I_h^{*}}{S}_m^{*}-\mu S_m^{*},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill 0& =B_s\frac{I_w^{*}}{S_w^{*}+I_w^{*}}{S}_m^{*}+B_{hm}\frac{I_h^{*}}{S_m^{*}+S_h^{*}+I_m^{*}+I_h^{*}}{S}_m^{*}-(\alpha +\mu )I_m^{*},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill 0& =\alpha p(I_h^{*}+I_w^{*}+I_m^{*})-(\delta +\mu )T^{*},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill 0& =\alpha (1-p)(I_h^{*}+I_w^{*}+I_m^{*})+\delta T^{*}-(d+\mu )A^{*}.\hfill \end{array}$$

(21)

where

$$\begin{array}{cc}\hfill {\mathsf{\Lambda }}_h& =\mathsf{\Psi }\theta (1-\gamma ),{\mathsf{\Lambda }}_w=\mathsf{\Psi }(1-\theta ),{\mathsf{\Lambda }}_m=\mathsf{\Psi }\theta \gamma ,B_h=c_h{\beta }_h,B_s=c_s{\beta }_s,\hfill \\ \hfill B_{hm}& =c_{hm}{\beta }_{hm},B_{hw}=c_{hw}{\beta }_{hw}.\hfill \end{array}$$

Thus, stationary points are:

• Disease-Free Equilibrium
  This happens when $I_h^{*}=I_w^{*}=I_m^{*}=0$ and represents absence of infection. It is:

  $$E_0=\left(\frac{{\mathsf{\Lambda }}_h}{\mu },0,\frac{{\mathsf{\Lambda }}_w}{\mu },0,\frac{{\mathsf{\Lambda }}_m}{\mu },0,0,0\right).$$

  (22)
• Boundary Equilibrium
  This occurs when $I_h^{*}=0$, the male homosexual population is null, and $I_w^{*},I_m^{*}$ are non-zero. The subscript ($\stackrel{¯*}{ }$) means the boundary equilibrium coordinate, which is:

  $$E_1=\left(\frac{{\mathsf{\Lambda }}_h}{\mu },0,\overline{S_w^{*}},\overline{I_w^{*}},\overline{S_m^{*}},\overline{I_m^{*}},\overline{T^{*}},\overline{A^{*}}\right),\mathrm{where}$$

  $$\begin{array}{cccc}\hfill \overline{S_w^{*}}& =\frac{{\mathsf{\Lambda }}_w}{B_s-\alpha },\hfill & \hfill \overline{S_m^{*}}& =\frac{{\mathsf{\Lambda }}_m}{B_s-\alpha },\hfill \\ \hfill \overline{I_w^{*}}& =\frac{{\mathsf{\Lambda }}_w}{B_s-\alpha }\left[R_0^{het}-1\right],\hfill & \hfill \overline{I_m^{*}}& =\frac{{\mathsf{\Lambda }}_m}{B_s-\alpha }\left[R_0^{het}-1\right],\hfill \\ \hfill \overline{T^{*}}& =\frac{\mathsf{\Psi }p\alpha \left[1-\theta (1-\gamma )\right]}{(B_s-\alpha )(\delta +\mu )}\left[R_0^{het}-1\right],\hfill & \hfill \overline{A^{*}}& =\frac{\mathsf{\Psi }\alpha \left[1-\theta (1-\gamma )\right]\left[\delta +\mu (1-p)\right]}{(B_s-\alpha )(\delta +\mu )(d+\mu )}\left[R_0^{het}-1\right].\hfill \end{array}$$

  (23)
  Note that the boundary equilibrium only exists when $R_0^{het}>1$ (implying $B_s>\alpha$).
• Endemic Equilibrium
  This represents persistence of the infection, it is:

  $$E_2=\left(S_h^{*},I_h^{*},S_w^{*},I_w^{*},S_m^{*},I_m^{*},T^{*},A^{*}\right),\mathrm{where}$$

  (24)

$$\begin{array}{cc}\hfill S_h^{*}=& \frac{{\mathsf{\Lambda }}_h}{B_h-\alpha },I_h^{*}=\frac{{\mathsf{\Lambda }}_h}{B_h-\alpha }\left[R_0^{hom}-1\right],\hfill \\ \hfill S_w^{*}=& \frac{{\mathsf{\Lambda }}_w-(\alpha +\mu )I_w^{*}}{\mu },S_m^{*}=\frac{{\mathsf{\Lambda }}_m-(\alpha +\mu )I_m^{*}}{\mu },\hfill \\ \hfill T^{*}=& \frac{\alpha p\left[{\mathsf{\Lambda }}_h(R_0^{hom}-1)+(B_h-\alpha )(I_w^{*}+I_m^{*})\right]}{(B_h-\alpha )(\delta +\mu )},\hfill \\ \hfill A^{*}=& \frac{\alpha \left(\delta +\mu (1-p)\right)\left[{\mathsf{\Lambda }}_h(R_0^{hom}-1)+(B_h-\alpha )(I_w^{*}+I_m^{*})\right]}{(B_h-\alpha )(\delta +\mu )(d+\mu )}.\hfill \end{array}$$

(25)

$$I_w^{*}=\frac{{\mathsf{\Lambda }}_w(\alpha I_m^{*}-{\mathsf{\Lambda }}_m)[B_s{I}_m^{*}(B_h-\alpha )+B_{hw}{\mathsf{\Lambda }}_h(R_0^{hom}-1)]-B_s{\mathsf{\Lambda }}_h{\mathsf{\Lambda }}_w\mu R_0^{hom}{I}_m^{*}}{(\alpha +\mu )\left\{\right(\alpha I_m^{*}-{\mathsf{\Lambda }}_m\left)\right[(B_h-\alpha )(I_m^{*}(B_s-\alpha )+{\mathsf{\Lambda }}_m)+B_{hw}{\mathsf{\Lambda }}_h(R_0^{hom}-1)]-{\mathsf{\Lambda }}_h\mu R_0^{hom}(I_m^{*}(B_s-\alpha )+{\mathsf{\Lambda }}_m\left)\right\}}.$$

(26)

$I_m^{*}$ is given by the roots of fourth degree polynomial:

$$p\left(I_m^{*}\right)=a_4{\left(I_m^{*}\right)}^4+a_3{\left(I_m^{*}\right)}^3+a_2{\left(I_m^{*}\right)}^2+a_1\left(I_m^{*}\right)+a_0.$$

(27)

Coefficients $a_0,a_1,a_2,a_3$, and $a_4$ are shown in Appendix B.

EE exists when $R_0^{hom}>1$. For infected males and females, the following inequalities (28) must be satisfied. Otherwise, the populations of susceptible men $S_m^{*}$ and women $S_w^{*}$ will be negatives:

$$0<I_m^{*}<\frac{{\mathsf{\Lambda }}_m}{\alpha +\mu }\mathrm{and}0<I_w^{*}<\frac{{\mathsf{\Lambda }}_w}{\alpha +\mu }.$$

(28)

Figure 3 shows the existence of equilibrium points, such as Disease-Free Equilibrium (DFE), Boundary Equilibrium (BE), and Endemic equilibrium (EE) in function of $R_0^{het}$ and $R_0^{hom}$. The figure shows two important aspects. First, the DFE is the only stationary point that exists when $R_0^{hom}$ or $R_0^{het}$ are less than one; it gives an idea of how stability can be. Second, for existence of EE, the $R_0^{hom}$ is more important that $R_0^{het}$ because when $R_0^{hom}$ is greater than 1 the EE exit, whereas when $R_0^{het}$ is greater than 1 is necessary that $R_0^{hom}$ will be greater than 1.

2.4. Local Stability of Disease-Free Equilibrium

Theorem 5.

The DFE $E_0=\left(\frac{{\mathsf{\Lambda }}_h}{\mu },0,\frac{{\mathsf{\Lambda }}_w}{\mu },0,\frac{{\mathsf{\Lambda }}_m}{\mu },0,0,0\right)$ is LAS if $R_0^{hom}<1$ and $R_0^{het}<1$ and is unstable when $R_0^{hom}>1$.

Proof. 

LAS will be demonstrated with the eigenvalues of the Jacobian matrix related to the system (1) to () evaluated in $E_0$, it is:

$$J\left(E_0\right)=\left(\begin{array}{cccccccc}-\mu & -B_h& 0& 0& 0& 0& 0& 0\\ 0& B_h-(\alpha +\mu )& 0& 0& 0& 0& 0& 0\\ 0& -B_{hw}\frac{1-\theta }{\theta }& -\mu & 0& 0& -B_s\frac{1-\theta }{\theta \gamma }& 0& 0\\ 0& B_{hw}\frac{1-\theta }{\theta }& 0& -(\alpha +\mu )& 0& B_s\frac{1-\theta }{\theta \gamma }& 0& 0\\ 0& -B_{hm}\gamma & 0& -B_s\frac{\theta \gamma }{1-\theta }& -\mu & 0& 0& 0\\ 0& B_{hm}\gamma & 0& B_s\frac{\theta \gamma }{1-\theta }& 0& -(\alpha +\mu )& 0& 0\\ 0& \alpha p& 0& \alpha p& 0& \alpha p& -(\delta +\mu )& 0\\ 0& \alpha (1-p)& 0& \alpha (1-p)& 0& \alpha (1-p)& \delta & -(d+\mu )\end{array}\right)$$

(29)

The characteristic polynomial is

$$p\left(\lambda \right)={(\lambda +\mu )}^3(\lambda +d+\mu )(\lambda +\delta +\mu )(\lambda +\alpha +\mu -B_h)\left[{(\alpha +\lambda +\mu )}^2-B_s^2\right].$$

Eigenvalues are:

$$\begin{array}{cccccc}\hfill {\lambda }_1& =-\mu \hfill & \hfill {\lambda }_2& =-\mu \hfill & \hfill {\lambda }_3& =-\mu ,\hfill \\ \hfill {\lambda }_4& =-(d+\mu ),\hfill & \hfill {\lambda }_5& =-(\delta +\mu ),\hfill & \hfill {\lambda }_6& =-(\alpha +\mu +B_s),\hfill \\ \hfill {\lambda }_7& =-(\alpha +\mu )\left[1-R_0^{het}\right],\hfill & \hfill {\lambda }_8& =-(\alpha +\mu )\left[1-R_0^{hom}\right].\hfill \end{array}$$

(30)

$R_0=max\{R_0^{hom},R_0^{het}\}<1$ imply $R_0^{hom}<1$ and $R_0^{het}<1$; thus, if all eigenvalues are negatives, it follows that $E_0$ is LAS. On the other hand, if $R_0>1$, then $R_0^{hom}>1$ or $R_0^{het}>1$, implying that ${\lambda }_7$ or ${\lambda }_8$, respectively, will be positive and in this case $E_0$ is unstable. □

2.5. Global Sensitivity Analysis

A sensitivity analysis will help us better understand which of the parameters in the model we should focus on estimating most precisely. We answer the following questions: Which parameters contribute most to output variability? Which parameters require additional research or are insignificant? These questions can be answered by performing an analysis with Latin Hypercube Sampling (LHS) and Partial Rank Coefficient (PRCC). We use Matlab to solve the system of ordinary differential equations and to implement most of the SA functions described throughout the manuscript; the functions are available at http://malthus.micro.med.umich.edu/lab/usanalysis.html (accessed on 20 April 2020).

LHS is a statistical sampling method that allows for an efficient analysis of parameter variations across simultaneous uncertainty ranges in each parameter [18]. PRCC shows which parameters have the largest influence on model outcomes [19]. To summarize, we can say that LHS is a sample method, and PRCC conducts the statistical treatment of each sample.

The model contains 12 parameters; however, to perform sensitivity analysis, only parameters related to HIV infection and related to a basic reproduction number are considered. They are: $\gamma ,p,\mu ,\alpha ,B_h,B_s,$$B_{hw},B_{hm}$. According to [20], a uniform distribution was chosen over a Gaussian (normal) one because we have no evidence of the ends of the ranges and we carry out multiple runs (NR = 300); parameters, baselines, ranges, and probability density functions (PDF) are listed in

Table 2. Parameters used in Sensitivity Analysis through Latin Hypercube Sampling and Partial Rank Correlation Coefficients (LHS/PRCC).

Parameter	Baseline	Range	PDF	Source
$\gamma$	0.9	$[0.3678,1]$	Uniform	Assumed
p	0.85	$[0.1353,1]$	Uniform	Assumed
$\mu$	0.0140	$[0.01,0.02]$	Uniform	[21]
$\alpha$	0.3333	$[0.1353,1]$	Uniform	[22]
$B_h$	2.64	$[0.05,3.95]$	Uniform	Assumed
$B_s$	0.04	$[0.0497,0.5]$	Uniform	[23]
$B_{hw}$	0.04	$[0.0497,0.5]$	Uniform	Assumed
$B_{hm}$	0.3	$[0.0497,0.5]$	Uniform	Assumed

. A Partial Rank Correlation Coefficient was created for each infected population. In addition, scatterplots for each of the aforementioned parameters are presented in Figure 4, Figure 5 and Figure 6.

3. Discussion

This analysis focuses on identifying the main parameters that play a dominant role in three different response outputs such as $I_h$, untreated infected homosexual men; $I_w$, untreated infected women; and $I_m$, untreated infected heterosexual men. The more sensitive parameters are: the departure rate of infected individuals, $\alpha$, and the infection rates in homosexuals and heterosexuals, $B_h$ and $B_s$, respectively. Scatterplots show the variation in the infected populations size with changes in parameters when examined, thus providing specific qualitative information on the relationship between an infected population and a parameter. Parameters with positive PRCCs will increase $I_{h,w,m}$ when their value is increased, whereas parameters with negative PRCCs will decrease $I_{h,w,m}$ when their value is increased. PRCC values are represented in Figure 4, Figure 5 and Figure 6.

It follows from Figure 4 that untreated infected homosexuals, $I_h$, have a negative correlation with the $\alpha$ parameter; in fact, the PRCC = $-0.98023$, which allows us to conclude that an increase in the $\alpha$ parameter means a decrease in the number of untreated infected homosexual men. In Figure 5, the number of untreated infected women, $I_w$, has a positive correlation with rate of infection in heterosexuals $B_s$. In fact, PRCC = $0.68725$, thus an increase in heterosexual contact implies an increase in women being infected; analogously, this population has a negative correlation with the $\alpha$ parameter. Figure 6 allows us to conclude that the $\gamma$ parameter does not influence untreated infected heterosexual men, $I_m$. In fact, PRCC = $-0.0063685$, showing that an increase in this parameter has little influence on the number of infected heterosexual men. In addition, infection rates such as $B_s$ and $B_h$ have similar behavior in infected heterosexual men, $I_m$, as for the infected women population, $I_w$.

4. Conclusions

Models will be a tool for understanding the disease dynamics and for predicting possible trends. Obviously, more accurate predictions require more complex models with more classes and compartments. Although such models are relatively easy to formulate, their mathematical analysis is difficult, and obtaining the necessary social and sexual behavior data is more complicated. Several key features could be included to create more realistic HIV/AIDS models in human populations, such as by looking infectious classes or transmission among injectable drug users through needle sharing.

Sensitivity analysis for the eight parameters related to infection population allow us to conclude that the most influential parameter in the HIV dynamic is the departure rate for infected individuals, $\alpha$, because it presents the highest PRCC coefficient. This behavior can be explained because the $\alpha$ parameter is present in the basic reproduction number and governs how those infected people are emerging from untreated status to obtain treatment or to develop AIDS.

Bisexual parameters, such as those of the probability of infection via sexual contact between homosexual men and heterosexual men, ${\beta }_{hm}$, and between homosexual men and women, ${\beta }_{hw}$, allow us to conclude that higher values of ${\beta }_{hm}$ and ${\beta }_{hw}$ imply a high infection rate in untreated infected women and heterosexual men.

Mathematics can provide information for the decision maker about how to promote awareness campaigns aimed at specific populations. This research allowed us to conclude that the best way to reduce contagion and consequently to reach a DFE is thought to be the reduction in homosexual partners rate, as they are the population most affected by the virus and are therefore the most likely to become infected and to spread it. Increasing the departure rate of infected individuals leads to a decrease in untreated infected heterosexual men and untreated infected women. However, it is not the only was to prevent and curb the rate of contagion in San Juan de Pasto. Consequently, it is also necessary to increase anti-retroviral treatment.

With the population parameters of San Juan de Pasto, several numerical simulations were performed by modifying parameters that make the basic reproduction number greater than or less than one. This seems to suggest that when $R_0^{het}<1$ and $R_0^{hom}>1$, there is a general decline in the rate of HIV infection over the next few years, but the infection persists. As a result, we can conclude that the most important observation from our findings is that, in the population, there is a short-term rise in HIV infection in which there exists a significant increase in new HIV infections, followed by a decline in the generation of new infections.

The dynamic of the HIV/AIDS epidemic, to a large extent, depends on changes in the basic reproduction number among homosexual men, $R_0^{hom}$, which was also evidenced by modifying several parameters in the scenarios above. In the background section, it was mentioned that the probability of HIV infection in homosexual men is great than in heterosexual people; thus, the basic reproduction number in heterosexual people, $R_0^{het}$, is less influential. In addition, investigations such as [24,25] permit us to conclude that the rate of sexual partners in homosexual men is greater than the rate of sexual partners in heterosexuals; thus, $R_0^{hom}>R_0^{het}$. This suggests that HIV infection can be controlled or eliminated from the community if control programs are directed towards reducing $R_0^{hom}$ to values less than one. The model shows the persistence of the disease when $R_0^{hom}>1$.

The dynamics of HIV/AIDS are, in general, too complex to allow for intuitive predictions and require the support of mathematical modeling for qualitatively and quantitatively assessing and understanding the functioning system. Furthermore, one of the most difficult tasks of mathematical modeling is obtaining parameters for a chosen model. Moreover, by using real parameter values to study and analyze the diverse sexual behavior in San Juan de Pasto, the proposed HIV/AIDS model tries to be as approximate as possible to the current situation of this infection. The emphasis was not on the accuracy of the scenarios, but on the actions that can be taken as a result of comprehending the state of the epidemic in the future. For example, scenario 5 shows that when the number of sexual partners is high, the basic reproduction number is greater than 1 and the infection spreads more easily, implying that more and more people are being treated with higher public health costs, and therefore, it is better and more economically efficient to invest in educational campaigns. These actions can involve, among other things, the prevention of new infections, the provision and delivery of anti-retroviral therapy, and educational campaigns such as those that aim to reduce the number of sexual partners or the use of condoms for self-protection.

This application in San Juan de Pasto shows the effects of modifying the parameters related to infected populations. These variations imply huge social and economic expenses which can and should be avoided through government actions such as educational campaigns. In this way, this research aims to be a useful tool in the design of establishing strategies for implementing valid public health policies and introducing efficient public health campaigns.