Mathematical Modeling of the HIV/Kaposi’s Sarcoma Coinfection Dynamics in Areas of High HIV Prevalence

Abstract

We formulate a deterministic system of ordinary differential equations to quantify HAART treatment levels for patients co-infected with HIV and Kaposi's Sarcoma in a high HIV prevalence setting. A qualitative stability analysis of the equilibrium states is carried out and we find that the disease-free equilibrium is globally attracting whenever the reproductive number ℛ_𝑘<1. A unique endemic equilibrium exists and is locally stable whenever ℛ_𝑘>1. Therefore, reducing ℛ_𝑘 to below unity should be the goal for disease eradication. Provision of HAART is shown to provide dual benefit of reducing HIV spread and the risk of acquiring another fatal disease for HIV/AIDS patients. By providing treatment to 10% of the HIV population, about 87% of the AIDS population acquire protection against coinfection with HIV and Kaposi's Sarcoma (KS). Most sub-Sahara African countries already have programmes in place to screen HIV. Our recommendation is that these programmes should be expanded to include testing for HHV-8 and KS counseling.

1. Introduction

Kaposi's Sarcoma is a cancer that occurs mostly in humans with suppressed immune systems [1]. The development of this cancer depends upon prior infection to the human herpesvirus-8 (HHV-8) [2], a virus which is usually transmitted either sexually or via saliva [3]. For HIV-related Kaposi's Sarcoma development, immunosuppression is a necessary causal factor [2]. Because HIV is an immunosuppressive virus, it promotes the development of Kaposi's Sarcoma in individuals dually infected with both viruses (HIV and HHV-8) and this combination has proved to be fatal and has made Kaposi's Sarcoma the fourth largest killer of people living with HIV/AIDS in sub-Sahara Africa [4].

In competent immune systems, acquisition of HHV-8 does not guarantee the development of KS; in fact, most individuals with a strong immune response could remain latently infected with HHV-8 throughout their lifetime [2]. The HIV-1 growth factors stimulate the immune cells including the healthy and infected B-cells to proliferate. The activation of latently infected B-cells to proliferate only leads to production of more HHV-8 that, according to the theory proposed by Foreman et al. [2], may be responsible for infection of progenitor cells of endothelial origin, which once infected with HHV-8 develop into Kaposi's Sarcoma cells [2].

Recent research findings [12] show that highly active antiretroviral therapy (HAART) for HIV significantly decreases (and, in some instances, completely forces the cancer into remission) KS activity in a patient, but such treatment is only effective against persons who have seroconverted. Lungu et al. studied the within-host dynamics of HIV and KS and came to the conclusion that HIV infection accentuates the potential for infected individuals to develop KS conditions and that administration of HAART on KS individuals who have seroconverted results in the reversal of KS conditions. However, the same result could not be demonstrated if HAART was administered to HIV negative individuals.

From what has been described above, we formulate a mathematical model which includes the following classes: a susceptible class, 𝑆, two infected classes of individuals infected with HIV-1 only, 𝐼, individuals coinfected with HHV-8 and HIV-1, 𝐼_𝑘, two classes of asymptomatic infectives with HIV only, 𝑃, asymptomatic coinfected, 𝑃_𝑘, a class of treated individuals, 𝑇, and two AIDS classes, namely, individuals with full-blown AIDS only, 𝐴, and individuals coinfected with full-blown AIDS and Kaposi's Sarcoma.

The structure of the paper is as follows. In Section 2 we formulate the model that describes the coinfection transmission dynamics. Existence of solutions is proved in Section 3, but our main interest is to simulate and study the dynamic behavior of the steady state solution. In particular, we want to quantify the proportion of HAART-treated individuals who never develop KS at various treatment rates 𝜙₁ and 𝜙₂.

AIDS is the progressive stage at which an HIV-infected individual loses competency of his immune system and becomes prone to opportunistic infections. For this reason, HIV coinfection models have found a lot of space in HIV epidemiology studies. HIV and tuberculosis coinfection have been studied by Cohen et al. [13], Roeger et al. [14], and Ramkissoon et al. [15] while Barley et al. [9], Chiyaka et al. [16], and Mukandavire et al. [17] studied HIV-malaria coinfection dynamics with each study giving results pertinent to the coinfection under investigation. The study of HIV-KS coinfection dynamics is still in infancy and, to the best of our knowledge, no mathematical modeling study has been carried out to assess the coinfection dynamics of HIV and KS at the population level in sub-Sahara Africa. Because our model incorporates HAART administration to all infective classes, it can provide insights into treatment strategies and, in particular, decide whether the current policy based on a CD4 count threshold to access treatment should be continued as a strategy.

2. Model Formulation

We begin with a human population of susceptibles which is free of both KS and HIV denoted by 𝑆. This population is replenished at constant rate Λ through sexual maturity or immigration. Upon effective contact with individuals infected with the HIV virus, the new infectives progress into the infected class 𝐼 at rate 𝜆_ℎ, where 𝜆_ℎ denotes the force of infection.

For simplicity of the model, we assume that all individuals in the class 𝑆 are latently infected with the HHV-8 virus which causes the KS infection. Therefore, in our model every class, except for the susceptible class, is assumed to be at risk of developing KS. HIV infection is known to promote or enhance the development of KS condition [18, 19], and so we assume that individuals in the class 𝐼 can develop mild KS at rate 𝜖₁ and move into the class of coinfected individuals KS denoted by 𝐼_𝑘.

HIV-only and coinfected individuals are assumed to progress to the asymptomatic pre-AIDS classes (𝑃 and 𝑃_𝑘) at the same constant rate 𝜓₁. Furthermore, persons in the pre-AIDS class are deemed to be sexually interacting and that individuals with mild KS in the 𝐼_𝑘 class can develop acute KS, which manifests in the form of visible lesions [19] and severe debilitation. Due to the nature of these symptoms, we assume that individuals who develop acute KS will no longer be sexually interacting. Those individuals in the 𝐼_𝑘 class, who develop acute KS, die at a disease induced rate 𝜏₁.

We assume that all infected individuals in the classes 𝐼 and 𝐼_𝑘 are tested and if they are tested positive for HIV-1 they receive treatment at rate 𝜙₁ and move into the class of treated individuals, 𝑇. Individuals in the pre-AIDS classes 𝑃 and 𝑃_𝑘 have higher viral loads and therefore possess less competent immune systems. These individuals will then present for medical attention at a higher rate 𝜙₂>𝜙₁. HAART is known to reverse KS conditions for people with HIV [19], and so in the model we assume that the therapy is perfect in reversing mild KS conditions. Upon accessing HAART, individuals in the infective class 𝐼_𝑘 and pre-AIDS class 𝑃_𝑘 achieve full KS recovery and progress to the treated class 𝑇. Additionally, individuals in the pre-AIDS class 𝑃_𝑘 progress to acute KS at rate 𝜏₂>𝜏₁. Individuals in the pre-AIDS class 𝑃 are at risk of developing KS at a constant rate 𝜖₂. Individuals in the pre-AIDS classes 𝑃 and 𝑃_𝑘 develop clinical symptoms and progress to full-blown AIDS 𝐴 and 𝐴_𝑘, respectively, at rate 𝜃₁. Progression to full-blown AIDS by persons in the class 𝑇 represents treatment failure, and they progress to the class 𝐴 at the same rate 𝜃₂. Full-blown AIDS individuals 𝐴 can also develop KS at rate 𝜖₃>𝜖₂>𝜖₁ due to weakened immune systems. Persons in the full-blown AIDS classes have additional AIDS-induced mortality 𝛿₁. Progression rate to acute KS for persons with full-blown AIDS is 𝜏₃, which we assume to be greater than 𝜏₂ and 𝜏₁. Individuals in both AIDS classes are also subject to a natural mortality rate of 𝜇.

The force of infection 𝜆_ℎ depends on the probability of transmission per contact 𝛽, the proportion of infected individuals in each category (𝐼 and 𝐼_𝑘), the proportion of infected individuals in receipt of HAART (𝑇), and the pre-AIDS classes (𝑃 and 𝑃_𝑘). Individuals in receipt of HAART have reduced viral load [20] and are therefore assumed to be less infectious relative to infectives not in receipt of HAART. This reduced infectiousness is modeled by 𝛾₁<1. Coinfected infectives and pre-AIDS individuals (𝐼_𝑘 and 𝑃_𝑘) have weaker immune systems and are therefore likely to carry higher viral loads than their counterparts in the (𝐼 and 𝑃) classes. This added infectiousness is modeled by 𝛾₂>1 and 𝛾₄>1. Pre-AIDS infectives 𝐴 are more infectious because of increased viral load, and this increased infectiousness is modeled by 𝛾₃ where 𝛾₁<1<𝛾₂<𝛾₃<𝛾₄. Persons with full-blown AIDS exhibit symptoms related to HIV and therefore are assumed to be noninteracting. The total sexually active variable population at time 𝑡 is given by 𝑁(𝑡)=𝑆(𝑡)+𝐼(𝑡)+𝐼_𝑘(𝑡)+𝑃(𝑡)+𝑃_𝑘(𝑡)+𝑇(𝑡). Assuming homogeneous mixing, the time dependent force of infection for HIV is given by 𝜆_ℎ𝛾=𝛽¹𝑇+𝐼+𝛾₂𝐼_𝑘+𝛾₃𝑃+𝛾₄𝑃_𝑘𝑆+𝑇+𝐼+𝐼_𝑘+𝑃(𝑡)+𝑃_𝑘. (1)

The model flow diagram depicting this biological system is illustrated in Figure 1.

Figure 1 

Model flow diagram.

The above assumptions and formulations lead us to this nonlinear system of differential equations: 𝑆^′(𝑡)=Λ−𝜆_ℎ𝑆(𝑡),𝐼+𝜇_′(𝑡)=𝜆_ℎ𝑆(𝑡)−𝜖₁+𝜇+𝜓₁+𝜙₁𝐼(𝑡),𝐼^′_𝑘(𝑡)=𝜖₁𝐼(𝑡)−𝜇+𝜓₁+𝜙₁+𝜏₁𝐼_𝑘(𝑡),𝑃^′(𝑡)=𝜓₁𝐼(𝑡)−𝜙₂+𝜖₂+𝜃₁𝑃(𝑡),𝑃+𝜇_′𝑘(𝑡)=𝜖₂𝑃(𝑡)+𝜓₁𝐼_𝑘−𝜙₂+𝜃₁+𝜇+𝜏₂𝑃_𝑘(𝑡),𝑇^′(𝑡)=𝜙₁𝐼(𝑡)+𝐼_𝑘(𝑡)+𝜙₂𝑃(𝑡)+𝑃_𝑘(𝑡)−𝜃₂𝑇(𝑡),𝐴+𝜇_′(𝑡)=𝜃₁𝑃(𝑡)+𝜃₂𝑇(𝑡)−𝜖₃+𝛿₁𝐴(𝑡),𝐴+𝜇_′𝑘(𝑡)=𝜃₁𝑃_𝑘(𝑡)+𝜖₃𝐴(𝑡)−𝛿₁+𝛿₂+𝜇𝐴_𝑘(𝑡). (2)

All parameters for the model system (2) are assumed to be nonnegative for all time 𝑡>0.

3. Basic Properties

3.1. Positivity and Boundedness of Solutions

We denote by 𝑅⁸₊ the set of points 𝑥_𝑡=(𝑥₁,𝑥₂,𝑥₃,𝑥₄,𝑥₅,𝑥₆,𝑥₇,𝑥₈) in 𝑅⁸ with positive coordinates and consider the system (2) with initial values 𝑥⁰=𝑥⁰₁,𝑥⁰₂,𝑥⁰₃,𝑥⁰₄,𝑥⁰₅,𝑥⁰₆,𝑥⁰₇,𝑥⁰₈. (3)

In this section, we prove the following lemma.

Lemma 1. The system (2) can be written as a system of differential inequalities 𝑑𝑥_𝑖𝑑𝑡≥𝐴_𝑖𝑥_𝑖+^𝑛_𝑗𝐵_𝑖𝑗𝑥_𝑗+𝜖(𝑖=1,…,𝑛), (4) where 𝐵_𝑖𝑗≥0,𝜖≥0. (5) If 𝑥_𝑖(0)≥𝜖 then 𝑥_𝑖(𝑡)≥0 for all 𝑡>0 and 1≤𝑖≤𝑛.

Without loss of generality we may assume that 𝜖>0, since the case 𝜖=0 follows by approximating the system with a sequence 𝜖=𝜖_𝑘↓0.

Proof. Suppose that the assertion 𝑥_𝑖(0)≥𝜖>0, for 1≤𝑖≤𝑛, is not true. Then there exists the smallest number 𝑡₀>0, such that 𝑥_𝑖(𝑡)>0for1≤𝑖≤𝑛,0≤𝑡≤𝑡₀,𝑥_𝑖𝑡₀=0foratleastone𝑖,say𝑖=𝑖₀. (6)
Then 𝑥_𝑖0 is a decreasing function at 𝑡=𝑡₀, so that 𝑑𝑥_𝑖0𝑡𝑑𝑡⁰≤0. (7)
From the differential inequality (4) for 𝑥_𝑖0(𝑡) we get 𝑑𝑥_𝑖0𝑡𝑑𝑡⁰≥^𝑛_𝑗𝐵_𝑖𝑗𝑥_𝑗+𝜖≥𝜖>0 (8) which is a contradiction.
For the state variables in our model, we always take 𝑆(0)≥0,𝐼(0)≥0,𝐼_𝑘(0)(0)𝑃≥0,𝑃≥0,_𝑘(0)≥0,𝑇(0)≥0,𝐴(0)≥0,𝐴_𝑘(0)≥0. (9)
From Lemma 1 we conclude that 𝑆(𝑡)≥0,𝐼(𝑡)≥0,𝐼_𝑘(𝑡)(𝑡)𝑃≥0,𝑃≥0,_𝑘(𝑡)≥0,𝑇(𝑡)≥0,𝐴(𝑡)≥0,𝐴_𝑘(𝑡)≥0. (10)
Thus, in the region 𝑅⁸₊ the model is epidemiologically and mathematically well posed and we can use it to study the dynamics of HIV-KS coinfection.

3.1.1. Equilibrium States, Reproductive Number, and Stability Analysis

The model (2) possesses a disease-free equilibrium, 𝜉₀, given by 𝜉₀=𝑆,𝐼,𝐼_𝑘,𝑃,𝑃_𝑘,𝑇,𝐴,𝐴_𝑘=Λ𝜇,0,0,0,0,0,0,0, (11) and at least one endemic equilibrium state whose existence is discussed in Section 3.3. Following Van Den Driessche and Watmough [21], the KS-induced reproductive number of system (2), ℛ_𝑘, is given by the spectral radius of the matrix 𝐹𝑉⁻¹ where the matrices 𝐹 and 𝑉 are given by 𝐹=𝛽𝛾₃𝛽𝛾₂𝛽𝛾₄𝛽𝛾₁𝛽,00000000000000000000𝑉=−𝐾₁𝜓0000₁−𝐾₂𝜖000₁0−𝐾₃000𝜖₂𝜓₁−𝐾₄0𝜙₁𝜙₂𝜙₁𝜙₂−𝐾₅, (12) where 𝐾₁=𝜇+𝜖₁+𝜙₁+𝜓₁,𝐾₂=𝜇+𝜖₂+𝜃₁+𝜙₂,𝐾₃=𝜇+𝜏₁+𝜙₁+𝜓₁,𝐾₄=𝜇+𝜃₁+𝜏₂+𝜙₂,𝐾₅=𝜇+𝜃₂. (13)

From (12) we have calculated the reproduction number ℛ_𝑘=ℛ_𝐼+ℛ_𝐼𝑘+ℛ_𝑃+ℛ_𝑃𝑘+ℛ_𝑇, (14) where ℛ_𝐼=𝛽𝜇+𝜖₁+𝜙₁+𝜓₁,ℛ_𝐼𝑘=𝛽𝛾₂𝜖₁𝜇+𝜖₁+𝜙₁+𝜓₁𝜇+𝜏₁+𝜙₁+𝜓₁,ℛ_𝑃=𝛽𝛾₃𝜓₁𝜇+𝜖₂+𝜃₁+𝜙₂𝜇+𝜖₁+𝜙₁+𝜓₁,ℛ_𝑃𝑘=𝛽𝛾₄𝜓₁×𝜖₁𝜇+𝜖₂+𝜃₁+𝜙₂+𝜖₂𝜇+𝜏₁+𝜙₁+𝜓₁×𝜇+𝜖₂+𝜃₁+𝜙₂𝜇+𝜃₁+𝜏₂+𝜙₂×𝜇+𝜖₁+𝜙₁+𝜓₁𝜇+𝜏₁+𝜙₁+𝜓₁⁻¹,ℛ_𝑇=ℛ_𝑇𝐼+ℛ_𝑇2+ℛ_𝑇3+ℛ_𝑇4+ℛ_𝑇5,ℛ_𝑇1=𝛽𝛾₁𝜙₁𝜇+𝜃₂𝜇+𝜖₁+𝜙₁+𝜓₁,ℛ_𝑇2=𝛽𝛾₁𝜙₂𝜓₁𝜇+𝜃₂𝜇+𝜖₂+𝜃₁+𝜙₂𝜇+𝜖₁+𝜙₁+𝜓₁,ℛ_𝑇3=𝛽𝛾₁𝜖₁𝜙₁𝜇+𝜃₂𝜇+𝜖₁+𝜙₁+𝜓₁𝜇+𝜏₁+𝜙₁+𝜓₁,ℛ_𝑇4=𝛽𝛾₁𝜖₁𝜙₂𝜓₁×𝜇+𝜃₂𝜇+𝜃₁+𝜏₂+𝜙₂×𝜇+𝜖₁+𝜙₁+𝜓₁𝜇+𝜏₁+𝜙₁+𝜓₁⁻¹,ℛ_𝑇5=𝛽𝛾₁𝜖₂𝜙₂𝜓₁×𝜇+𝜃₂𝜇+𝜖₂+𝜃₁+𝜙₂×𝜇+𝜃₁+𝜏₂+𝜙₂𝜇+𝜖₁+𝜙₁+𝜓₁⁻¹. (15)

This number ℛ_𝑘 is a threshold such that if ℛ_𝑘<1 the disease clears from the population. If ℛ_𝑘>1 the steady state 𝜉₀ becomes unstable and the disease establishes itself into the population. This number is comprehensively analyzed further in Section 3.4 to reveal the impact of treatment.

3.2. Global Stability of the Disease-Free Equilibrium

If ℛ_𝑘<1 then fixed point ℰ₀ is locally asymptotically stable. We now determine conditions which guarantee global asymptotic stability of the disease-free state [22]. Rewrite model system (2) as 𝑑𝑋(),𝑑𝑡=𝐹𝐱,𝑍𝑑𝑍()()𝑑𝑡=𝐺𝑋,𝑍,𝐺𝐱,0=0, (16) where 𝑋∈ℝ^𝑚 denotes the number of uninfected individuals and 𝑍∈ℝ^𝑛 denotes the number of infected individuals including those latently infected and those who are infectious. The disease-free equilibrium state can now be written as 𝑈₀=𝑋^∗,0, (17) where 𝑋^∗=((Λ/𝜇),0,0,0,0,0,0,0,0). To guarantee global stability, the following conditions must be satisfied: (i)for 𝑑𝑋/𝑑𝑡=𝐹(𝑋,0), 𝑋^∗ is globally asymptotically stable; (ii)𝐺(𝑋,𝑍)=A𝑍−𝐺(𝑋,𝑍),𝐺(𝑋,𝑍)≥0 for (𝑋,𝑍)∈Ω; (iii)𝐴=𝐷_𝑍𝐺(𝑋^∗,0) is an 𝑀-matrix (the off diagonal elements of 𝐴 are nonnegative) and Ω is the region where the model makes biological sense.

Lemma 2. The fixed point 𝑈₀=(𝑋^∗,0) is a globally asymptotically stable point of (2) provided that ℛ_𝑘<1 and conditions stated above are satisfied.

Proof. Consider 𝐹(𝑋,0)=Λ−𝜇𝑆_𝑓0000,𝐺(𝑋,𝑍)=𝐴𝑍−𝐺(𝑋,𝑍),𝐴=𝛽−𝐾₁𝛾₃𝛽𝛾2𝛽𝛾₄𝛽𝛾₁𝛽𝜓₁−𝐾₂𝜖000₁0−𝐾₃000𝜖₂𝜓₁−𝐾₄0𝜙₁𝜙₂𝜙₁𝜙₂−𝐾₅. (18)
Then 𝐺(𝑋,𝑍)=𝛽𝛾₁𝑇+𝐼_𝑚+𝛾₂𝐼_𝑘+𝛾₃𝑃+𝛾₄𝑃_𝑘×𝑆1−𝑆+𝐼+𝐼_𝑘+𝑃+𝑃_𝑘0000. (19)
Since 0≤𝑆≤𝑆+𝐼+𝐼_𝑘+𝑃+𝑃_𝑘, it is clear that 𝐺(𝑋,𝑍)≥0 and 𝐴 is an 𝑀-matrix whenever its spectral radius is less than unity. Then, by Lemma 2, 𝜉₀ is a globally asymptotically stable equilibrium of (2).

The biological interpretation of Lemma 2 is that this population is observing the one partner policy or abstinence from unprotected sex resulting in HIV being eliminated from the population regardless of the size of the initial subpopulations [23].

3.3. Existence of Endemic Equilibria and Bifurcation Analysis

In this section following the approach in [7], we establish conditions for the existence of endemic equilibria and investigate their stability. Denote the arbitrary endemic equilibrium of model (2) by ℰ^∗ where ℰ^∗=𝑆^∗,𝐼^∗,𝐼^∗_𝑘,𝑇^∗𝑃^∗,𝑃^∗,𝐴^∗,𝐴^∗_𝑘. (20)

Let 𝜆^∗_ℎ𝛾=𝛽¹𝑇^∗+𝐼^∗+𝛾₂𝐼^∗_𝑘+𝛾₃𝑃^∗+𝛾₄𝑃^∗_𝑘𝑆^∗+𝑇^∗+𝐼^∗+𝐼^∗_𝑘+𝑃^∗(𝑡)+𝑃^∗_𝑘, (21) be the force of infection. By setting the right hand side of model (2) to zero, we obtain, in terms of the force of infection, 𝑆^∗=Λ𝜆^∗_ℎ+𝜇,𝐼^∗=𝜆^∗_ℎ𝐾₁,𝐼^∗_𝑘=𝜖₁𝐼^∗𝐾₂,𝑇^∗=𝜓₁𝐼^∗𝐾₃,𝑃^∗=𝜓₁𝐼^∗_𝑘+𝜖₂𝑃^∗𝐾₄,𝐴^∗=𝜃₁𝐼^∗+𝜃₂𝑇^∗𝜖₂+𝜇+𝛿₁,𝐴_𝑘=𝜃₁𝐼_𝑘+𝜖₃𝐴_𝑘𝜇+𝛿₁+𝛿₂. (22)

Substituting the values of 𝑆^∗,𝐼^∗,𝐼^∗_𝑘,𝑇^∗𝑃^∗,𝑃^∗,𝐴^∗, and 𝐴^∗_𝑘 into (21) we obtain the polynomial 𝑎₀𝜆²_ℎ+𝐾₁𝜆_ℎ=0, (23) where 𝑎₀=𝐾₂𝐾₄𝐾₃+𝜖₁𝜇+𝜃₂+𝜙₁+𝐾₂𝜖₁+𝐾₃𝐾₄+𝜖₂𝜇+𝜃₂+𝜙₂𝜓₁,𝑎₁=−𝛽𝛾₁𝐾₃+𝜖₁+𝐾₂𝜖₁+𝐾₃𝐾₄+𝜖₂𝜙₂𝜓₁+𝜇+𝜃₂−𝛽𝐾₃𝐾₄𝛾₃+𝛾₄𝜖₂𝜓₁−𝐾₂𝛽−𝐾₁𝐾₃𝐾₄+𝛽𝜖₁𝐾₄𝛾₂+𝛾₄𝜓₁. (24) Equation (23) has two solutions given by 𝜆_ℎ0𝜆=0,_ℎ1=𝐾₁𝐾₂𝐾₃𝐾₄𝐾₂𝐾₄ℛ_𝐾−1𝐾₆𝐾₈+ℛ_𝐾𝛾₁𝐾₉𝜙₁+ℛ_𝐾−1𝐾₁₀𝐾₈+ℛ_𝐾𝛾₁𝐾₇𝜙₂𝜓₁×𝐾₂𝐾₄𝐾₉𝐾₈+𝜙₁+𝐾₇𝐾₈+𝜙₂𝜓₁×𝐾₅𝜓₁+𝐾₂𝐾₄𝐾₆+𝛾₄𝜖₁𝜓₁⁻¹, (25) where 𝐾₅=𝐾₃𝐾₄𝛾₃+𝛾₄𝜖₂,𝐾₆=𝐾₃+𝛾₂𝜖₁,𝐾₈=𝜇+𝜃₂,𝐾₉=𝐾₃+𝜖₁,𝐾₁₀=𝐾₂𝛾₄𝜖₁+𝐾₅. (26)

From the solution 𝜆_ℎ1, the condition ℛ_𝐾>1 is necessary for existence of the nontrivial endemic equilibrium and we summarize this as follows.

Lemma 3. The model (2) has a unique endemic equilibrium whenever ℛ_𝑘>1.

Having proved the existence of the endemic equilibrium point we now investigate its stability using the Centre Manifold Theory [24] as described by Castillo-Chavez and Song [25]. Consider the following general system of ordinary differential equations with a parameter 𝜙𝑑𝑥𝑑𝑡=𝑓𝑥,𝜙,𝑓∶ℝ^𝑛×ℝ⟶ℝ,𝑓∈𝐶²(ℝ×ℝ). (27)

Without loss of generality, it is assumed that 0 is an equilibrium for system (27) for all values of the parameter 𝜙, (i.e., 𝑓(𝜙,0)≡0 for all 𝜙).

Assume the following.(A1):𝐴=𝐷_𝑥𝑓(0,0)=((𝜕𝑓_𝑖/𝜕𝑥_𝑗)(0,0)) is the linearise matrix of system (27) around the equilibrium 0 with 𝜙 evaluated at 0. Zero is a simple eigenvalue of 𝐴 and all other eigenvalues of 𝐴 have negative real parts. (A2): Matrix 𝐴 has a nonnegative right eigenvector 𝑤 and a left eigenvector 𝑣 corresponding to the zero eigenvalue.

Let 𝑓_𝑘 be the 𝑘th component of 𝑓 and 𝑎=^𝑛_{𝑘,𝑖,𝑗=1}𝑣_𝑘𝑤_𝑖𝑤_𝑗𝜕²𝑓_𝑘𝜕𝑥_𝑖𝜕𝑥_𝑗(0,0),𝑏=^𝑛_𝑘,𝑖=1𝑣_𝑘𝑤_𝑖𝜕²𝑓_𝑘𝜕𝑥_𝑖().𝜕𝜙0,0 (28)

The local dynamics of system (27) around 0 are totally determined by 𝑎 and 𝑏. (i)𝑎>0, 𝑏>0. When 𝜙<0 with |𝜙|≪1, 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0<𝜙≪1, 0 is unstable and there exists a negative and locally asymptotically stable equilibrium. (ii)𝑎<0, 𝑏<0. When 𝜙<0 with |𝜙|≪1, 0 is unstable; when 0<𝜙≪1, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium. (iii)𝑎>0, 𝑏<0. When 𝜙<0 with |𝜙|≪1, 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when 0<𝜙≪1, 0 is stable, and a positive unstable equilibrium appears. (iv)𝑎<0, 𝑏>0. When 𝜙 changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly, a negative unstable equilibrium becomes positive and locally asymptotically stable. Particularly, if 𝑎>0 and 𝑏>0, then, a backward bifurcation occurs at 𝜙=0.

We make the following change of variables: 𝑆=𝑥₁, 𝐼=𝑥₂, 𝐼_𝑘=𝑥₃, 𝑃=𝑥₄, 𝑃_𝑘=𝑥₅, 𝑇=𝑥₆, 𝐴=𝑥₇, and 𝐴_𝑘=𝑥₈ so that 𝑁=∑⁸_𝑛=1𝑥_𝑛. We now use the vector notation 𝑋=(𝑥₁,𝑥₂,𝑥₃,𝑥₄,𝑥₅,𝑥₆,𝑥₇,𝑥₈)^𝑇. Then, model system (2) can be written in the form 𝑑𝑋/𝑑𝑡=𝐹=(𝑓₁,𝑓₂,𝑓₃,𝑓₄,𝑓₅,𝑓₆,𝑓₇,𝑓₈)^𝑇, such that 𝑥^′₁(𝑡)=𝑓₁=Λ−𝜆_ℎ+𝜇𝑥₁(𝑡),𝑥^′₂(𝑡)=𝑓₂=𝜆_ℎ𝑥₁(𝑡)−𝜖₁+𝜇+𝜓₁+𝜙₁𝑥₂(𝑡),𝑥^′₃(𝑡)=𝑓₃=𝜖₁𝑥₂(𝑡)−𝜇+𝜓₁+𝜙₁+𝜏₁𝑥₃(𝑡),𝑥^′₄(𝑡)=𝑓₄=𝜓₁𝑥₂(𝑡)−𝜙₂+𝜖₂+𝜃₁+𝜇𝑥₄(𝑡),𝑥^′₅(𝑡)=𝑓₅=𝜖₂𝑥₄(𝑡)+𝜓₁𝑥₃−𝜙₂+𝜃₁+𝜇+𝜏₂𝑥₅(𝑡),𝑥^′₆(𝑡)=𝑓₆=𝜙₁𝑥₂(𝑡)+𝑥₃(𝑡)+𝜙₂𝑥₄(𝑡)+𝑥₅(𝑡)−𝜃₂+𝜇𝑥₆(𝑡),𝑥^′₇(𝑡)=𝑓₇=𝜃₁𝑥₄(𝑡)+𝜃₂𝑥₆(𝑡)−𝜖₃+𝛿₁+𝜇𝑥₇(𝑡),𝑥^′₈(𝑡)=𝑓₈=𝜃₁𝑥₅(𝑡)+𝜖₃𝑥₇(𝑡)−𝛿₁+𝛿₂+𝜇𝑥₈(𝑡), (29) where 𝜆_ℎ=𝛾₁𝑥₆+𝑥₂+𝛾₂𝑥₃+𝛾₃𝑥₄+𝛾₄𝑥₅𝑥₁+𝑥₂+𝑥₃+𝑥₄+𝑥₅+𝑥₆. (30)

The Jacobian matrix of system (29) at the disease-free equilibrium is given by 𝐽𝜉₀=𝛽−𝐾₁𝛾₃𝛽𝛾2𝛽𝛾₄𝛽𝛾₁𝛽𝜓₁−𝐾₂𝜖000₁0−𝐾₃000𝜖₂𝜓₁−𝐾₄0𝜙₁𝜙₂𝜙₁𝜙₂−𝐾₅, (31) from which it can be shown that ℛ_𝑘 is the same as in (14). If 𝛽^∗ is taken as a bifurcation point and if we consider the case ℛ₁=1 and solve for 𝛽^∗=𝛽 gives 𝛽^∗=1ℛ_𝑘. (32)

The linearised system of the transformed equations (29) at 𝛽^∗ has a simple zero eigenvalue. Hence, the Centre Manifold Theory [24] can be used to analyse the dynamics of system (29) near 𝛽^∗. It can be shown that the Jacobian of (29) at 𝛽^∗ has a right eigenvector associated with the zero eigenvalue given by 𝑢=[𝑢₁,𝑢₂,𝑢₃,𝑢₄,𝑢₅,𝑢₆,𝑢₇,𝑢₈]^𝑇, where 𝑢₂=𝑢₂>0,𝑢₃=𝑢₂𝜓₁𝐾₂𝑢>0,₄=𝑢₂𝜖₁𝐾₃>0,𝑢₆=𝑢₂𝜙₁+𝑢₄𝜙₁+𝑢₃+𝑢₅𝜙₂𝐾₅𝑢>0,₇=𝑢₂𝜃₁+𝑢₆𝜃₂𝐾₆>0,𝑢₈=𝑢₇𝜖₃+𝑢₄𝜃₁𝐾₇𝑢>0,₁𝛽𝑢=−²+𝑢₆𝛾₁+𝑢₄𝛾₂+𝑢₃𝛾₃+𝑢₅𝛾₄𝜇<0. (33)

The left eigenvector of 𝐽(𝜉₀) associated with the zero eigenvalue at 𝛽=𝛽^∗ is given by 𝑣=[𝑣₁,𝑣₂,𝑣₃,𝑣₄,𝑣₅,𝑣₆,𝑣₇,𝑣₈] where 𝑣₁=𝑣₇=𝑣₈=0,𝑣₂=𝑣₂𝑣>0,₃=𝛽𝑣₂𝛾₃+𝑣₅𝜖₂+𝑣₆𝜙₂𝐾₂>0,𝑣₄=𝛽𝑣₂𝛾₂+𝑣₆𝜙₁+𝑣₅𝜓₁𝐾₃,𝑣₅=𝛽𝑣₂𝛾₄+𝑣₆𝜙₂𝐾₄,𝑣₆=𝛽𝑣₂𝛾₁𝐾₅. (34)

For the model system (2), the associated nonzero partial derivatives of 𝐹 at the disease-free equilibrium are given by 𝜕𝑓²₂𝜕𝑥₂𝜕𝑥₂=−2𝛽𝜇Λ,𝜕𝑓²₂𝜕𝑥₂𝜕𝑥₃=𝜕𝑓²₂𝜕𝑥₂𝜕𝑥₃=−𝛽𝜇1+𝛾₂Λ,𝜕𝑓²₂𝜕𝑥₂𝜕𝑥₄=𝜕𝑓²₂𝜕𝑥₄𝜕𝑥₂=−𝛽𝜇1+𝛾₃Λ,𝜕𝑓²₂𝜕𝑥₂𝜕𝑥₅=−𝜕𝑓²₂𝜕𝑥₅𝜕𝑥₂=𝛽𝜇1+𝛾₄Λ,𝜕𝑓²₂𝜕𝑥₂𝜕𝑥₆=𝜕𝑓²₂𝜕𝑥₆𝜕𝑥₂=−𝛽𝜇𝑢₆1+𝛾₁Λ,𝜕𝑓²₂𝜕𝑥₃𝜕𝑥₃=−2𝛽𝜇𝛾₂Λ,𝜕𝑓²₂𝜕𝑥₃𝜕𝑥₄=𝜕𝑓²₂𝜕𝑥₄𝜕𝑥₃𝛾=−𝛽𝜇₂+𝛾₃Λ,𝜕𝑓²₂𝜕𝑥₃𝜕𝑥₅=𝜕𝑓²₂𝜕𝑥₅𝜕𝑥₃𝛾=−𝛽𝜇₂+𝛾₄Λ,𝜕𝑓²₂𝜕𝑥₃𝜕𝑥₆=𝜕𝑓²₂𝜕𝑥₆𝜕𝑥₃=𝛽𝜇𝛾₁+𝛾₂Λ,𝜕𝑓²₂𝜕𝑥₄𝜕𝑥₄=−2𝛽𝜇𝛾₃Λ,𝜕𝑓²₂𝜕𝑥₄𝜕𝑥₅=𝜕𝑓²₂𝜕𝑥₅𝜕𝑥₄=𝛽𝜇𝛾₃+𝛾₄Λ,𝜕𝑓²₂𝜕𝑥₄𝜕𝑥₆=𝜕𝑓²₂𝜕𝑥₆𝜕𝑥₄𝛾=−𝛽𝜇₁+𝛾₃Λ,𝜕𝑓²₂𝜕𝑥₅𝜕𝑥₅=2𝛽𝜇𝛾₄Λ,𝜕𝑓²₂𝜕𝑥₅𝜕𝑥₆=𝜕𝑓²₂𝜕𝑥₆𝜕𝑥₅=𝛽𝜇𝛾₁+𝛾₄Λ,𝜕𝑓²₂𝜕𝑥₆𝜕𝑥₆=2𝛽𝜇𝛾₁Λ. (35)

From (28), it follows that 1𝑎=−Λ𝐾²₂𝐾²₃𝐾²₄𝐾²₅×2𝛽𝜇𝑢²₂𝑣₂𝐾₂𝐾₄𝐾₃+𝜖₁𝐾₅+𝜙₁+𝐾₂𝜖₁+𝐾₃𝐾₄+𝜖₂𝐾₅+𝜙₂𝜓₁×𝐾₃𝐾₅𝑎₄𝛾₂+𝛾₄𝜖₂+𝛾₁𝐾₄+𝜖₂𝜙₂𝜓₁+𝐾₂𝐾₄𝐾₅𝐾₃+𝛾₃𝜖₁+𝛾₁𝐾₃+𝜖₁𝜙₁+𝜖₁𝐾₅𝛾₄+𝛾₁𝜙₂𝜓₁<0. (36)

For the sign of 𝑏, it is associated with the nonvanishing partial derivatives of 𝐹, 𝜕²𝑓₂𝜕_𝑥2𝜕_𝛽𝑚𝜕=1,²𝑓₂𝜕_𝑥3𝜕_𝛽𝑚=𝛾₂,𝜕²𝑓₂𝜕_𝑥4𝜕_𝛽𝑚=𝛾₃,𝜕²𝑓₂𝜕_𝑥5𝜕_𝛽𝑚=𝛾₄,𝜕²𝑓₂𝜕_𝑥6𝜕_𝛽𝑚=𝛾₁. (37)

From (28), it follows that 𝑏=𝑣₂𝑢₂𝜕²𝑓₂𝜕_𝑥2𝜕_𝛽𝑚+𝑢₃𝜕²𝑓₂𝜕_𝑥3𝜕_𝛽𝑚+𝑢₄𝜕²𝑓₂𝜕_𝑥4𝜕_𝛽𝑚+𝑢₅𝜕²𝑓₂𝜕_𝑥5𝜕_𝛽𝑚+𝑢₆𝜕²𝑓₂𝜕_𝑥6𝜕_𝛽𝑚>0. (38)

By Theorem 11 in [25] we establish result in Lemma 4.

Lemma 4. Since 𝑎<0, then model system (2) has a forward transcritical bifurcation and a unique locally stable endemic equilibrium (ℰ^∗) guaranteed by Lemma 3 exists for ℛ_𝑘>1 but sufficiently close to 1.

3.4. Analysis of the Reproductive Number ℛ_𝑘

To analyze the coinfection dynamics of KS and HIV, we investigate the KS-induced reproductive number in (14). In (14), we note that the partial reproductive numbers ℛ_𝐼, ℛ_𝐼𝐾, ℛ_𝑃, ℛ_𝑃𝐾, and ℛ_𝑇 represent the contribution to the reproduction number ℛ_𝑘 from the following classes of infectives HIV only (𝐼), coinfected (𝐼_𝑘), Pre-AIDS (𝑃), coinfected pre-AIDS (𝑃_𝑘), and those receiving HAART treatment (𝑇). The partial reproductive number ℛ_𝑇 representing the contribution of groups in receipt of ART is split into ℛ_𝑇1, ℛ_𝑇2, ℛ_𝑇3, ℛ_𝑇4, and ℛ_𝑇5 and represents the contribution of HIV infectives, infectives dually infected with KS, Pre-AIDS infectives, and Pre-AIDS infectives dually infected with KS, respectively, who receive anti-retroviral therapy.

If we define ℛ₀ to be the basic HIV reproductive number in the absence of KS and HIV treatment and set parameters related to KS and treatment to zero in (14) we obtain ℛ_{𝐾∣𝜖1=𝜖2=𝜙1=𝜙2=0}=ℛ₀=ℛ_𝐼+ℛ_𝑃. (39)

If we define the partial reproductive number without treatment to be ℛ_𝑁𝑇 and set parameters related to treatment to zero then we obtain ℛ_{𝐾∣𝜙1=𝜙2=0}=ℛ_𝑁𝑇,=ℛ_𝐼+ℛ_𝐼𝐾+ℛ_𝑃+ℛ_𝑃𝐾>ℛ₀. (40)

If ℛ₀>1, then not providing treatment to all infective classes worsens the HIV in this population. Note that even if ℛ₀<1, the option not to provide treatment to all infective classes could increase the reproduction number above 1. Treatment to all infective classes is necessary to effectively control the HIV disease spread. Using sensitivity analysis on ℛ_𝐾, we have reinforced some of these findings as shown in Figure 2.

(a)

(b)

(c)

(d)

Figure 2 

(a) 3D plot of ℛ𝐾 against the rates of KS manifestation for infected and pre-AIDS classes for variable values of 𝜖1 and 𝜖2, all other parameters constant as in Table 1, (b) contour plot map of ℛ𝐾 against 𝜖1 and 𝜖2, all other parameters constant as in Table 1, (c) 3D plot of ℛ𝐾 against the rates of ART administration for infected and pre-AIDS classes for variable values of 𝜙1 and 𝜙2, all other parameters constant as in Table 1, and (d) contour plot map of ℛ𝐾 against 𝜙1 and 𝜙2, all other parameters constant as in Table 1.

From Figures 2(a) and 2(b) wem note that an increase in the rate of infected people who develop KS will lead to an increase in ℛ_𝐾 and consequently an increase in the HIV epidemic. However, an increase in the rate of progress to dual infection for people in the pre-AIDS class will result in a marginal increase of ℛ_𝐾 and hence the HIV epidemic.

Figures 2(c) and 2(d) show that treatment of infected individuals and dually infected individuals will have a much higher positive impact on ℛ_𝐾 and the epidemic than treatment at the pre-AIDS stage and this leads us to conclude that early administration of HAART on HIV infectives to curtail the growth of opportunistic infections, such as KS, will have a more positive impact than delayed therapy at the pre-AIDS stage.

4. Numerical Simulations

Using the R programming environment, we ran numerical simulations of the model. The following data were input as initial conditions: 𝑁₀=𝐼,𝐼_𝑘,𝑃,𝑃_𝑘,𝑇,𝐴,𝐴_𝑘,𝑆=[1,0,0,0,0,0,0,800]. (41)

Parameter values used in the numerical simulations of model system (2) are shown in Table 1. HIV/KS coinfection study is still in its infancy. A number of numeric values for the parameters shown in Table 1 are reasonable estimates. Some of the parameters were estimated by modifying baseline values from published literature and these are denoted with an asterisk in Table 1.

Table 1 

Parameter values and their estimates.

Our value for the rate of HIV transmission is the average of the minimum (0.011) and maximum values (0.95) for the same parameter in Baggaley et al. [5] and Boily et al. [6]. We assumed that a death rate due to AIDS and KS coinfection would be 0.4; in our model, 𝛿₁ and 𝛿₂ are summed. Our recruitment rate is calculated by incorporating an annual recruitment rate (0.029) from Malunguza et al. [7]. In this study, we take the population of sub-Sahara Africa to be 767 million as in Mukandawire et al. [16] and Barley et al. [13].

Figure 3(a) shows that in the absence of treatment, the number of individuals with the AIDS/KS coinfection exceeds the number of those with just AIDS; that is, 54.0% of the total AIDS population has KS by the end of the model. Merely affecting a treatment rate of 1%, is enough to reverse this relation as depicted in Figure 3(b) where 61.7% of all AIDS patients are devoid of KS. If the treatment rate is increased to 10%, the AIDS population without the coinfection increases to 85.7% as depicted in Figure 3(c). Using parameter values in Table 1 causes the model to tend towards some endemic equilibrium over time. In other words, the DFE is unstable, and so we expect ℛ₀>1. In fact, carrying out these computations with 𝜙₁=𝜙₂=0.1, we obtain ℛ₀≈4.55. This value is slightly higher than estimated values of ℛ₀ for HIV/AIDS in European countries, which is consistent with our expectations [26].

(a)

(b)

(c)

(d)

Figure 3 

Epidemic curves for our model with the following treatment scenarios: (a) no treatment (𝜙1=𝜙2=0), (b) 1% treatment level (𝜙1=𝜙2=0.01), and (c) 10% treatment level (𝜙1=𝜙2=0), (d) treatment rate versus AIDS population without coinfection. A treatment level of 10% is sufficient at minimizing endemic populations in all the infective, untreated classes at the equilibrium.

Figure 3(d) provides an interesting insight into the predictions of our model. By providing treatment to 100% of all infectives and pre-AIDS individuals the model fails to eradicate the coinfection, and approximately an 8% prevalence of the coinfection among the full-blown AIDS cohort exists. This result means that we need to be cognisant of two objectives, namely, that we should have benchmark target goals in the provision of treatment and that work should be done to decrease the number of patients entering the pre-AIDS coinfection class. It is entirely unrealistic to presume that the coinfection can be eliminated entirely based on providing treatment to infectives and pre-AIDS individuals alone. In fact, Figure 3(d) shows that providing a treatment rate of 𝜙₁=𝜙₂=0.2 results in about 90% of the AIDS population losing the coinfection. Providing any further increases to 𝜙₁ and 𝜙₂ results in only a 2% increase.

Our model makes no attempt to consider the cost associated with providing the kind of treatment we have discussed. However, it is noteworthy to consider that the gains associated with providing treatment to any more than 20% of all infectives and pre-AIDS individuals are minimal. It is clear that there is an 8% gap between where we would like to see the percentage of the full-blown AIDS cohort without coinfection and where our model currently estimates that percentage to be. Assuming that it is unrealistic to completely eradicate the coinfection, steps should be taken to shift the curve in Figure 3(d) closer to 100%. We believe that a good way to do this is by screening patients for the coinfection once the HIV infection is first diagnosed. In general, this will serve to decrease the total number of patients who become pre-AIDS with the coinfection and, in turn, this will decrease the total number of individuals presenting with the coinfection.

5. Discussion

In numerous HIV-positive cohorts, susceptibility to KS development is incredibly high. Recent findings that antiretroviral treatment for HIV clinical symptoms can reverse KS in most patients suggest that the only obstacle in preventing KS-related complications (or even deaths) is one's ability to access the treatment itself. This is especially the case in sub-Saharan Africa, where KS prevalence is higher than anywhere else in the world, yet access to proper treatment remains low. We have shown that providing treatment to just 10% of HIV-infected individuals, regardless of their KS status, can offer a significant reduction in the overall number of individuals who end up with full-blown AIDS and KS. Certainly, we expect that, as we increase the treatment rate, the disparity will become more and more favorable. However, our results only apply to situations in which individuals develop KS on their own—that is, in the absence of interaction with HHV-8-infected individuals. While such a generalization may be perfectly valid in sub-Saharan Africa, it will not apply to other regions where HHV-8 needs to be acquired sexually or through saliva. In other words, this model is very specific to the sub-Saharan region.

In the future, we hope to extend this model to consider all modes of HHV-8 transmission. One way to do this is by introducing separate classes of HHV-8 infectious individuals. This will increase the overall reach of the model and hopefully provide a positive influence on policy making in areas affected by KS. For now, supporting KS education and awareness needs to become as important as providing universal access to antiretroviral treatment for HIV patients. Most southern African countries already have programs in place to encourage HIV screening; this is the perfect time to also test for HHV-8 and provide KS counseling. Regardless of how long the HIV/AIDS epidemic lasts in Africa, there is no reason people should have to fight KS as well.

Acknowledgments

This work was done under SAMSA Masamu Program with support from the National Science Foundation under Grant no. DMS-1050259. Lungu acknowledges the support of the University of Botswana, Massaro would like to acknowledge partial support from the NSF SCALE-IT graduate fellowship program at the University of Tennessee, and Malunguza acknowledges the support of the National University of Science and Technology.