Effect of innate and adaptive immune mechanisms on treatment regimens in an AIDS-related Kaposi's Sarcoma model

Abstract

Kaposi Sarcoma (KS) is the most common AIDS-defining cancer, even as HIV-positive people live longer. Like other herpesviruses, human herpesvirus-8 (HHV-8) establishes a lifelong infection of the host that in association with HIV infection may develop at any time during the illness. With the increasing global incidence of KS, there is an urgent need of designing optimal therapeutic strategies for HHV-8-related infections. Here we formulate two models with innate and adaptive immune mechanisms, relevant for non-AIDS KS (NAKS) and AIDS-KS, where the initial condition of the second model is given by the equilibrium state of the first one. For the model with innate mechanism (MIM), we define an infectivity resistance threshold that will determine whether the primary HHV-8 infection of B-cells will progress to secondary infection of progenitor cells, a concept relevant for viral carriers in the asymptomatic phase. The optimal control strategy has been employed to obtain treatment efficacy in case of a combined antiretroviral therapy (cART). For the MIM we have shown that KS therapy alone is capable of reducing the HHV-8 load. In the model with adaptive mechanism (MAM), we show that if cART is administered at optimal levels, that is, 0.48 for protease inhibitors, 0.79 for reverse transcriptase inhibitors and 0.25 for KS therapy, both HIV-1 and HHV-8 can be reduced. The predictions of these mathematical models have the potential to offer more effective therapeutic interventions in the treatment of NAKS and AIDS-KS.

Keywords:

• HIV-1
• HHV-8
• mathematical model
• optimal control

1. Introduction

Despite significant progress made in ending the HIV/AIDS epidemic, an estimated 38 million people were living with HIV at the end of 2018, resulting in about 2$\mathrm{\%}$ deaths. The African region remains to be the most affected, accounting for two-third of the people living with HIV worldwide (https://www.who.int/gho/hiv/en/). Although HIV-positive people who start antiretroviral therapy (HAART) have the same life expectancy as their HIV negative peers, they develop co-morbidities on average 16 years earlier than HIV negative people (http://www.natap.org/2020/CROI/croi_134). KS is one of the most common malignancies causing co-morbidity in patients with human immunodeficiency virus-1 (HIV-1) infection, especially at the later HIV stage (AIDS). Most of AIDS-related cancers are caused by oncogenic viruses such as Epstein Barr virus (EBV), human herpesvirus 8 (HHV-8) and Human papillomavirus (HPV) [7].

There are four different forms of KS: Classic or sporadic KS, African or Endemic KS, AIDS-associated or epidemic KS and Transplant or Immunosuppression-associated or Iatrogenic KS [13]. The development of each of these forms is dependent on prior infection with HHV-8. However, HHV-8 infection alone is insufficient for the development of KS and some form of immunodeficiency is necessary for disease progression [31].

Most individuals infected with African KS and Classic KS but with strong immune responses have remained latently infected with HHV-8 throughout their lifetime [13,14] The co-factors involved in the development of Classic and Endemic KS are not fully understood although environmental and genetic factors such as age, sex, malnutrition and so on have been implicated [13]. Progression from HHV-8 infection to KS is a complex process. For instance, not every AIDS patient develops KS even in the face of profound immunosuppression, only a minority of HHV-8-infected transplant recipients develop iatrogenic KS, and that people with Classic or Endemic KS are not typically immunosuppressed [17,20].

Whether HHV-8 infection develops into an asymptomatic or symptomatic KS, depends on the interplay between HHV-8 and the host immune system. When HHV-8 infection occurs, the immune system promotes an environment where cellular proliferation, cell migration, angiogenesis and cytokine/chemokine production are enhanced [13]. The immune response occurs in two stages: first by triggering the innate response and second, if the infection persists, the adaptive response [30]. A review by Foreman et al. [13] has suggested how infection of progenitor cells by HHV-8 can initiate the development of all forms of KS. For individuals dually infected with both HIV-1 and HHV-8, the HHV-8 infection is enhanced by the HIV-1 growth factors which stimulate both uninfected and infected B-cells to proliferate in response to T-cell signals [13]. The T-cell signals stimulate the latently infected B cells. These cells that were dormant are now capable to proliferate and increase the population of HHV-8 producing cells.

With regard to mitigating the spread of the disease, especially in the case for childhood diseases and malaria, preventive measures are given priority over treatment. Current protocol advises individuals going to malaria endemic areas to take malaria prophylaxis drugs 1 week before departure to prepare their immune system to fight and clear the infection before it develops into active disease. In this study, our objective is to demonstrate how administration of HAART to individuals co-infected with HIV-1 and HHV-8 can prevent the occurrence of KS by ensuring low HIV-1 viremia which prevents reactivation of latently infected B cells [13].

When a pathogen invades the body, the body triggers an innate, non-specific immune response to clear the infection. This response consists of cellular (immune cells) and chemical (e.g. cytokines) defenses to reduce the growth of the population of infected cells and to eliminate the pathogens. The innate immune response may be viewed as a way to suppress and control HHV-8 infection before the adaptive immune response characterized by the clonal expansion of lymphocytes is activated. Using mathematical modelling, we show that a dynamic motif in Figure 1 comprising of interactions between infected B cells, infected progenitor cells, KS cells, HHV-8 virions and the innate immune response, is able to prevent a potentially NAKS from developing into a clinical disease. This has twofold implications. First, key innate immune signalling molecules induced by viral infection lead to the production of a broad range of antiviral proteins and cytokines. Uncontrolled release of these cytokines can lead to cytokine storm, causing tissue damage or indirectly causing pathology even before the initiation of HAART [23]. Second, antibody test can detect HIV infection as early as 1–2 weeks after exposure and testing after 2 or 3 weeks is not very useful (http://i-base.info/guides/testing/what-is-the-window-period). Hence, it is essential to understand at what level to deem the innate immune response or cytokine therapy to be safe.

Figure 1. Schematic diagram of the MIM describing interactions for NAKS.

If the infection progresses despite the innate response, the immune system mounts a more robust, longer lasting adaptive or acquired immune response. Hence, we construct a second model that mimics the body's adaptive immune response by including the interactions as in Figure 4 between HIV-1 virions, HIV-1- and HHV-8-specific effector cells, infected CD4 T cells, uninfected B- and CD4 T cells. The initial condition of MAMl will be determined by the equilibrium states of the MIM, assuming an advanced stage of HIV-1 and HHV-8 co-infection. Importantly, we find an infectivity threshold that will be critical for the primary HHV-8 infection to develop into an advanced KS. To determine the drug efficacy level of HAART alone or combined HAART and chemotherapy, we will take an optimal control approach. This is motivated by the fact that HAART should be the first step therapy in optimal control of HIV infection for AIDS-KS. However, patients with high-risk KS rarely respond to HAART alone and hence, chemotherapy is recommended which requires balancing the immunosuppressive effects of chemotherapy with its potential benefit.

To develop these models, we shall apply the Foreman et al. [13] approach. According to this hypothesis, AIDS-KS arises from the erroneous infection of progenitor cells by HHV-8 which is enhanced by action of HIV-1 infected host cells. These HIV-1 infected cells produce cytokines and growth factors that stimulate the progenitors of the KS cells which makes them susceptible to HHV-8 infection. In summary, we will show that early HHV-8-specific intervention is important as it can control the HHV-8 infection from developing into a progressive KS. We also determine efficacy levels for cART therapy at which HIV-1 and HHV-8 co-infection can be kept under control, thus providing valuable testable predictions for clinical researchers.

2. Model with innate mechanism (MIM)

2.1. Model formulation and description

We formulate a model based on Foreman et al. [13] representing two subsystems as follows: the first subsystem representing the primary infection of B cells leading to the production of HHV-8 and the second subsystem representing the erroneous infection of progenitor cells leading to the development of KS.

The MIM includes infected B-cells, $X_1\left(t\right)$, HHV-8 virions, $X_2\left(t\right),$ infected progenitor cells, $X_3\left(t\right)$, KS cells, $X_4\left(t\right)$, and the innate immune response, $X_5\left(t\right)$. The interaction among the different classes are illustrated in Figure 1 and described by the following system of ordinary differential equations: (1) $\begin{array}{r}{\dot{X}}_1\left(t\right)=K_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_1\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}{X}_1.\end{array}$(1) Equation (1(1) $\begin{array}{r}{\dot{X}}_1\left(t\right)=K_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_1\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}{X}_1.\end{array}$(1) ) describes the dynamics of the infected B cells, $X_1$. This class is assumed to grow logistically but regulated by the efficacy threshold of the innate immune response. The dependence on $X_1$ itself rather than the HHV-8 load is plausible, since no correlation has been observed between the B-cell subsets and the HHV-8 viremia [6]. The last term accounts for natural death of infected B cells at a constant rate ${\mu }_{x_1}$. (2) $\begin{array}{r}{\dot{X}}_2\left(t\right)=N_{x_2}{\mu }_{x_1}{X}_1-{\mu }_{x_2}{X}_2.\end{array}$(2) Equation (2(2) $\begin{array}{r}{\dot{X}}_2\left(t\right)=N_{x_2}{\mu }_{x_1}{X}_1-{\mu }_{x_2}{X}_2.\end{array}$(2) ) represents the dynamics of HHV-8, $X_2$. The first term represents the production of these virions from the bursting of infected B cells, where $N_{x_2}$ is the carrying capacity or maximum number of virions that can be contained within an infected B cell. The last term represents natural clearance of HHV-8 at a constant rate ${\mu }_{x_2}$. (3) $\begin{array}{r}{\dot{X}}_3\left(t\right)=K_1\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_2\left(1-\frac{X_2}{x_{2_{max}}}\right)-{\mu }_{x_3}{X}_3.\end{array}$(3) Equation (3(3) $\begin{array}{r}{\dot{X}}_3\left(t\right)=K_1\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_2\left(1-\frac{X_2}{x_{2_{max}}}\right)-{\mu }_{x_3}{X}_3.\end{array}$(3) ) describes the dynamics of the infected progenitor cells, $X_3$. The first term represents a source term which grows logistically with respect to the viral level, $X_2$, and moderated by the innate immune response, $X_5$. The effect of the innate immune response, $X_5$, is moderated by the saturation parameter, ${\theta }_{x_5}$, which is significant in this study as it mimics how administration of vaccines or drugs can alter the progression of the infection [19]. The logistic growth term is expressed in terms of HHV-8 viremia to emphasize that progenitor cells and in general stem cells proliferate in response to infectious stimuli [3,5,8,15,19,24,27]. This formulation can assist to make decisions on viral load dependent intervention measures depending on the viremia reservoir levels. The second term represents the blanket death of these cells at a constant rate, ${\mu }_{x_3}$. (4) $\begin{array}{r}{\dot{X}}_4\left(t\right)={\hat{\mu }}_{x_3}{X}_3-{\mu }_{x_4}{X}_4.\end{array}$(4) Equation (4(4) $\begin{array}{r}{\dot{X}}_4\left(t\right)={\hat{\mu }}_{x_3}{X}_3-{\mu }_{x_4}{X}_4.\end{array}$(4) ) represents the concentration of KS, $X_4$. The first term designates the growth of KS, as the infected progenitor cells transform into cancerous cells, at a constant rate, ${\hat{\mu }}_{x_3}<{\mu }_{x_3}$. It is assumed that ${\hat{\mu }}_{x_3}<{\mu }_{x_3}$ as not all infected progenitor cells progress to KS [13]. The second term is natural death of KS at a constant rate ${\mu }_{x_4}$. (5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) Equation (5(5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) ) represents the innate immune response, $X_5$. The first term represents the stimulation of the innate immunity due to the presence of infected B cells, $X_1$. It is assumed that $X_5$ is stimulated by the infected B cells, $X_1$, in a saturable manner with the scaling constant, ${\theta }_{x_1},$ and decays at a constant rate $K_3$ [4]. The model (1(1) $\begin{array}{r}{\dot{X}}_1\left(t\right)=K_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_1\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}{X}_1.\end{array}$(1) )–(5(5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) ) is developed to demonstrate how infection errors committed by HHV-8 by erroneously infecting progenitor cells lead to a more serious problem of KS. This model can then be used to demonstrate that externally administered drugs or immune boosters can alter the infection and stop the development of KS.

2.2. Analysis of the model

2.2.1. Positivity and boundedness of solutions

We denote by ${\mathbb{R}}_{+}^5$ the set of points $X_t=\left(X_1\right(t),X_2(t),X_3(t),X_4(t),X_5(t\left)\right)$ in ${\mathbb{R}}^5$ with positive coordinates and consider the system (1(1) $\begin{array}{r}{\dot{X}}_1\left(t\right)=K_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_1\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}{X}_1.\end{array}$(1) )–(5(5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) ) with initial values $X^0=\left(X_1^0,X_2^0,X_3^0,X_4^0,X_5^0\right)\in {\mathbb{R}}_{+}^5.$In this section, we prove the following theorem.

Theorem 2.1

If $X_i^0\ge 0,$ then $X_i\left(t\right)\ge 0$ for all $t>0,i=1,\dots ,5$.

Before we prove Theorem 2.1, we rearrange the system into a subsystem of infected progenitor and infected B cells, $W=\left(X_1\right(t),X_3(t){)}^T,$ written in matrix form as (6) $\begin{array}{r}\dot{W}=M_{w}W+Q,\end{array}$(6) where (7) $\begin{array}{r}\left\{\begin{array}{l}{M}_w=\left(\begin{array}{cc}-{\mu }_{x_1}& 0\\ 0& -{\mu }_{x_3}\end{array}\right),F=\left(\begin{array}{c}{X}_1{f}_1\left(X_1\right)\\ X_2{f}_2\left(X_2\right)\end{array}\right),\mathrm{a}\mathrm{n}\mathrm{d}W=\left(\begin{array}{c}{X}_1\\ X_3\end{array}\right),\\ K=\left(\begin{array}{cc}{K}_2& 0\\ 0& K_1\end{array}\right),f_j\left(X_j\right)=\left(1-\frac{X_j}{x_{j_{max}}}\right),j\in \{1,2\},\\ r\left(X_5\right)=1-\frac{X_5}{{\theta }_{x_5}+X_5},Q=r\left(X_5\right)KF,\\ F_j\left(X_j\right)=X_j{f}_j\left(X_j\right)\text{ is the }\mathit{\text{j}}\text{th entry of the vector},F,j=1,2\end{array}\right.\end{array}$(7) Let $Q_j$ denote the jth entry of $Q,j=1,2,$ where $Q_j$ represents the jth source term in (6(6) $\begin{array}{r}\dot{W}=M_{w}W+Q,\end{array}$(6) ) and $r\left(X_5\right)$ measures the efficacy of the innate immune response. Define a subsystem consisting of HHV-8 and KS, $Y=\left(Y_2\right(t),Y_4(t){)}^T=(X_2\left(t\right),X_4\left(t\right){)}^T,$ which can be expressed in matrix form as (8) $\begin{array}{r}\dot{Y}=M_{y}Y+EW,\end{array}$(8) where $M_y=\left(\begin{array}{cc}-{\mu }_{y_2}& 0\\ 0& -{\mu }_{y_4}\end{array}\right),Y=\left(\begin{array}{c}{Y}_2\\ Y_4\end{array}\right),\mathrm{a}\mathrm{n}\mathrm{d}E=\left(\begin{array}{cc}{N}_{x_2}{\mu }_{x_1}& 0\\ 0& {\hat{\mu }}_{x_3}\end{array}\right).$

The proof is done in three steps: First, we prove that the source terms in (6(6) $\begin{array}{r}\dot{W}=M_{w}W+Q,\end{array}$(6) ) are nonnegative, i.e. $Q:=r\left(X_5\right)KF\ge 0.$ Second, we want to show that $X_i\left(t\right)\ge 0,i=1,3,t>0$ and finally, we conclude that $Y_i\left(t\right)\ge 0$ for $t>0,i=2,4.$

Proof.

From Equation (5(5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) ), we can deduce that $X_5\left(t\right)\ge X_5\left(0\right)\exp (-K_{3}t)\ge 0$. Notice that $r\left(X_5\right)=1-\frac{X_5}{{\theta }_{x_5}+X_5}=\frac{{\theta }_{x_5}}{{\theta }_{x_5}+X_5}>0$, for ${\theta }_{x_5}>0$.

The function $F_j\left(X_j\right)$ in (7(7) $\begin{array}{r}\left\{\begin{array}{l}{M}_w=\left(\begin{array}{cc}-{\mu }_{x_1}& 0\\ 0& -{\mu }_{x_3}\end{array}\right),F=\left(\begin{array}{c}{X}_1{f}_1\left(X_1\right)\\ X_2{f}_2\left(X_2\right)\end{array}\right),\mathrm{a}\mathrm{n}\mathrm{d}W=\left(\begin{array}{c}{X}_1\\ X_3\end{array}\right),\\ K=\left(\begin{array}{cc}{K}_2& 0\\ 0& K_1\end{array}\right),f_j\left(X_j\right)=\left(1-\frac{X_j}{x_{j_{max}}}\right),j\in \{1,2\},\\ r\left(X_5\right)=1-\frac{X_5}{{\theta }_{x_5}+X_5},Q=r\left(X_5\right)KF,\\ F_j\left(X_j\right)=X_j{f}_j\left(X_j\right)\text{ is the }\mathit{\text{j}}\text{th entry of the vector},F,j=1,2\end{array}\right.\end{array}$(7) ) has zeros at $X_j=0$ and $X_j=x_{j_{max}},$ has a peak at $X_j=\frac{x_{j_{max}}}{2}$ and is positive in the interval $0<X_j<x_{j_{max}}$. Since $K\ge 0,$ we have $Q_j\ge 0,j=1,2$.

The matrix $M_w$ in (6(6) $\begin{array}{r}\dot{W}=M_{w}W+Q,\end{array}$(6) ) is a Mertzler matrix and $f_j\left(X_j\right)\ge 0$, the solution of (6(6) $\begin{array}{r}\dot{W}=M_{w}W+Q,\end{array}$(6) ) is nonnegative for all $t>0.$ The matrix $E\ge 0$ since $X_i\ge 0$ for $i=1,3.$ The matrix $M_y$ is a Mertzler matrix and hence, the solution of (8(8) $\begin{array}{r}\dot{Y}=M_{y}Y+EW,\end{array}$(8) ) is nonnegative for $t>0.$ We conclude that if $X_i^0\ge 0,i=1,2,3,4,5$, then the solution $X_i\left(t\right)$, of the system (1(1) $\begin{array}{r}{\dot{X}}_1\left(t\right)=K_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_1\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}{X}_1.\end{array}$(1) )–(5(5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) ) remains in ${\mathbb{R}}_{+}^5$.

2.3. Steady states and the basic reproduction number

The virus free equilibrium of the MIM given by Equations (1(1) $\begin{array}{r}{\dot{X}}_1\left(t\right)=K_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_1\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}{X}_1.\end{array}$(1) )–(5(5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) ) is ${ϵ}^0=(X_1^0,X_2^0,X_3^0,X_4^0,X_5^0)=(0,0,0,0,0)$. In what follows, we will calculate the basic reproduction number of the system (1(1) $\begin{array}{r}{\dot{X}}_1\left(t\right)=K_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_1\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}{X}_1.\end{array}$(1) )–(5(5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) ) using the next generation operator method [9]. The basic reproduction number is determined by the number of newly infected B cells. Using this approach, we first assume that the model system (1(1) $\begin{array}{r}{\dot{X}}_1\left(t\right)=K_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_1\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}{X}_1.\end{array}$(1) )–(5(5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) ) can be written in the form (9) $\begin{array}{r}\begin{array}{rlr}\frac{\mathrm{d}X}{\mathrm{d}t}& =& f\left(\mathbf{X}\mathbf{,}\mathbf{Y}\mathbf{,}\mathbf{Z}\right),\\ \frac{\mathrm{d}Y}{\mathrm{d}t}& =& g\left(\mathbf{X}\mathbf{,}\mathbf{Y}\mathbf{,}\mathbf{Z}\right),\\ \frac{\mathrm{d}Z}{\mathrm{d}t}& =& h\left(\mathbf{X}\mathbf{,}\mathbf{Y}\mathbf{,}\mathbf{Z}\right),\end{array}\end{array}$(9) where $\mathbf{X}\in \mathbb{R},\mathbf{Y}\in \mathbb{R}$ and $\mathbf{Z}\in {\mathbb{R}}^3$, and $h(\mathbf{X},0,\mathbf{0})=\mathbf{0}.$ Assuming that the equation $g\left({\mathbf{X}}^{\mathbf{\ast }}\mathbf{,}\mathbf{Y}\mathbf{,}\mathbf{Z}\right)=0$ implicitly determines a function $Y=\tilde{g}\left({\mathbf{X}}^{\mathbf{\ast }}\mathbf{,}\mathbf{Z}\right).$ We let $\mathbf{A}=D_{Z}h({\mathbf{X}}^{\mathbf{\ast }},\tilde{g}({\mathbf{X}}^{\mathbf{\ast }}\mathbf{,}\mathbf{0}),0)$ and further assume that $\mathbf{A}$ can be written in the form $\mathbf{A}\mathbf{=}\mathbf{C}\mathbf{-}\mathbf{D},$ with $\mathbf{C}\ge \mathbf{0}$ (that is $m_{ij}\ge 0$) and $\mathbf{D}\ge \mathbf{0}$ is a diagonal M-matrix.

In system (9(9) $\begin{array}{r}\begin{array}{rlr}\frac{\mathrm{d}X}{\mathrm{d}t}& =& f\left(\mathbf{X}\mathbf{,}\mathbf{Y}\mathbf{,}\mathbf{Z}\right),\\ \frac{\mathrm{d}Y}{\mathrm{d}t}& =& g\left(\mathbf{X}\mathbf{,}\mathbf{Y}\mathbf{,}\mathbf{Z}\right),\\ \frac{\mathrm{d}Z}{\mathrm{d}t}& =& h\left(\mathbf{X}\mathbf{,}\mathbf{Y}\mathbf{,}\mathbf{Z}\right),\end{array}\end{array}$(9) ), $\mathbf{X}$ denotes the innate immune response, $\mathbf{Y}$ represents the HHV-8 virions and the components of $\mathbf{Z}$ represent the HHV-8-associated cells, i.e. $\mathbf{X}=X_5,\mathbf{Y}=X_2,\mathbf{Z}=(X_1,X_3,X_4)$. Let ${\mathbf{U}}_0=({\mathbf{X}}^{\mathbf{\ast }},0,0)$ denote the virus free equilibrium, that is, $f\left({\mathbf{X}}^{\mathbf{\ast }}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{0}\right)=g\left({\mathbf{X}}^{\mathbf{\ast }}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{0}\right)=0,\mathrm{a}\mathrm{n}\mathrm{d}h\left({\mathbf{X}}^{\mathbf{\ast }}\mathbf{,}\mathbf{0}\mathbf{,}\mathbf{0}\right)=\mathbf{0},\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}Y=\tilde{g}\left({\mathbf{X}}^{\mathbf{\ast }}\mathbf{,}\mathbf{Z}\right),$where $\tilde{g}\left({\mathbf{X}}^{\mathbf{\ast }}\mathbf{,}\mathbf{Z}\right)=\frac{N_{x_2}{\mu }_{x_1}{X}_1}{{\mu }_{x_2}}\mathrm{a}\mathrm{n}\mathrm{d}\mathbf{D}=diag\left({\mu }_{x_1},{\mu }_{x_3},{\mu }_{x_4}\right).$We compute $\mathbf{A}=D_{Z}h({\mathbf{X}}^{\mathbf{\ast }},\tilde{g}({\mathbf{X}}^{\mathbf{\ast }}\mathbf{,}\mathbf{0}),0)$ and get $\mathbf{C}=\left(\begin{array}{ccc}{K}_2& 0& 0\\ \frac{K_1{N}_{x_2}{\mu }_{x_1}}{{\mu }_{x_2}}& 0& 0\\ 0& {\hat{\mu }}_{x_3}& 0\end{array}\right)\mathrm{a}\mathrm{n}\mathrm{d}{\mathbf{D}}^{-1}=\left(\begin{array}{ccc}\frac{1}{{\mu }_{x_1}}& 0& 0\\ 0& \frac{1}{{\mu }_{x_3}}& 0\\ 0& 0& \frac{1}{{\mu }_{x_4}}\end{array}\right)$The reproduction number is given by the next generation spectral radius $\rho \left(\mathbf{C}{\mathbf{D}}^{\mathbf{-}\mathbf{1}}\right)$ to be ${\mathcal{R}}_0=\frac{K_2}{{\mu }_{x_1}}.$

2.4. Local stability of the virus free equilibrium, ${ϵ}^0$

Lemma 2.1

The virus free equilibrium point ${ϵ}^0$ is locally asymptotically stable if ${\mathcal{R}}_0<1$.

Proof.

We consider the Jacobian matrix of the system (1(1) $\begin{array}{r}{\dot{X}}_1\left(t\right)=K_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_1\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}{X}_1.\end{array}$(1) )–(5(5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) ), evaluated at the virus free steady state denoted by $J\left({ϵ}^0\right)$. $J\left({ϵ}^0\right)=\left(\begin{array}{ccccc}{K}_2-{\mu }_{x_1}& 0& 0& 0& 0\\ N_{x_2}{\mu }_{x_1}& -{\mu }_{x_2}& 0& 0& 0\\ 0& K_1& -{\mu }_{x_3}& 0& 0\\ 0& 0& {\hat{\mu }}_{x_3}& -{\mu }_{x_4}& 0\\ 0& 0& \frac{K_3}{{\theta }_{x_1}}& 0& -K_3\end{array}\right)$We note that $J\left({ϵ}^0\right)$ is a lower triangular matrix. Hence, the corresponding eigenvalues are the entries in the main diagonal. In other words, ${\lambda }_1=K_2-{\mu }_{x_1},{\lambda }_2=-{\mu }_{x_2},{\lambda }_3=-{\mu }_{x_3},{\lambda }_4=-{\mu }_{x_4},{\lambda }_5=-K_3.$For local stability of ${ϵ}^0,{\lambda }_1=K_2-{\mu }_{x_1}={\mu }_{x_1}({\mathcal{R}}_0-1)<0.$ Hence, all the eigenvalues are negative and the result follows.

2.5. Existence of the KS present equilibrium, ${ϵ}_1^{\ast }$

Setting the system (1(1) $\begin{array}{r}{\dot{X}}_1\left(t\right)=K_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_1\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}{X}_1.\end{array}$(1) )–(5(5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) ) to zero and solving the resulting system simultaneously yields: $X_1=0\mathrm{o}\mathrm{r}{K}_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}=0,$Suppose $X_1\ne 0$. Then, (10) $\begin{array}{r}{K}_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}=0,\end{array}$(10) Using (5(5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) ) to solve for $X_5$ and replacing in (10(10) $\begin{array}{r}{K}_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}=0,\end{array}$(10) ) yields $A_2{X}_1^2+A_1{X}_1+A_0=0,$where $A_2=K_2{\theta }_{x_5}$, $A_1={\mu }_{x_1}{x}_{1_{max}}+K_2{\theta }_{x_1}{\theta }_{x_5}+{\mu }_{x_1}{\theta }_{x_5}{x}_{1_{max}}(1-{\mathcal{R}}_0)$ and $A_0={\mu }_{x_1}{\theta }_{x_1}{\theta }_{x_5}{x}_{1_{max}}(1-{\mathcal{R}}_0)$. To establish the existence of a positive root for $g\left(X_1\right)=A_2{X}_1^2+A_1{X}_1+A_0$, say $X_1^{\ast },$ we argue as follows:

Note that $g\left(0\right)={\mu }_{x_1}{\theta }_{x_1}{\theta }_{x_5}{x}_{1_{max}}(1-{\mathcal{R}}_0)<0\mathrm{i}\mathrm{f}{\mathcal{R}}_0>1$ and by continuity of g, we have $\underset{X_1\to \mathrm{\infty }}{lim}g\left(X_1\right)=+\mathrm{\infty }.$This implies that there is a positive number, $X_1^{\ast }\in (0,+\mathrm{\infty })$ such that $g\left(X_1^{\ast }\right)=0$. In particular, we one can show that $X_1^{\ast }=\frac{-A_1+\sqrt{A_1^2-4A_0{A}_2}}{2A_2}$. Hence, we obtain the following coordinates for the KS present equilibrium, ${ϵ}_1^{\ast }=(X_1^{\ast },X_2^{\ast },X_3^{\ast },X_4^{\ast },X_5^{\ast })$, where (11) $\begin{array}{rl}{X}_2^{\ast }& =\frac{{\mu }_{x_1}{x}_{2_{max}}{\mathcal{R}}_f}{K_2},\end{array}$(11) (12) $\begin{array}{rl}{X}_3^{\ast }& ={\left(x_{2_{max}}{\mathcal{R}}_f\right)}^2\frac{{\mu }_{x_1}}{{\mu }_{x_3}{x}_{2_{max}}}\left(\frac{{\theta }_{x_5}\left({\theta }_{x_1}+X_1^{\ast }\right)}{{\theta }_{x_5}\left({\theta }_{x_1}+X_1^{\ast }\right)+X_1^{\ast }}\right)\left(\frac{1}{{\mathcal{R}}_f}-\frac{1}{{\mathcal{R}}_0}\right),\end{array}$(12) (13) $\begin{array}{rl}{X}_4^{\ast }& ={\left(x_{2_{max}}{\mathcal{R}}_f\right)}^2\frac{{\mu }_{x_1}{\hat{\mu }}_{x_3}}{{\mu }_{x_3}{\mu }_{x_4}{x}_{2_{max}}}\left(\frac{{\theta }_{x_5}\left({\theta }_{x_1}+X_1^{\ast }\right)}{{\theta }_{x_5}\left({\theta }_{x_1}+X_1^{\ast }\right)+X_1^{\ast }}\right)\left(\frac{1}{{\mathcal{R}}_f}-\frac{1}{{\mathcal{R}}_0}\right)\end{array}$(13) (14) $\begin{array}{rl}{X}_5^{\ast }& =\frac{X_1^{\ast }}{{\theta }_{x_1}+X_1^{\ast }},\end{array}$(14) (15) $\begin{array}{rl}{\mathcal{R}}_f& =\frac{K_2{N}_{x_2}{X}_1^{\ast }}{{\mu }_{x_2}{x}_{2_{max}}}.\end{array}$(15)

Theorem 2.2

The endemic equilibrium, ${ϵ}_1^{\ast },$ exists if ${\mathcal{R}}_0>max\{1,{\mathcal{R}}_f\}$.

We define ${\mathcal{R}}_f$ as the infectivity resistance threshold which must be exceeded for the infection of progenitor cells to occur. The concept of pathogen load in relationship to infectivity is discussed in many studies (see , e.g.  [21]).

Remark 2.1

From (11(11) $\begin{array}{rl}{X}_2^{\ast }& =\frac{{\mu }_{x_1}{x}_{2_{max}}{\mathcal{R}}_f}{K_2},\end{array}$(11) )–(15(15) $\begin{array}{rl}{\mathcal{R}}_f& =\frac{K_2{N}_{x_2}{X}_1^{\ast }}{{\mu }_{x_2}{x}_{2_{max}}}.\end{array}$(15) ), we deduce the following scenarios:

1. If ${\mathcal{R}}_f<1<{\mathcal{R}}_0,$ then the endemic steady state ${ϵ}_1^{\ast }=(X_1^{\ast },X_2^{\ast },X_3^{\ast },X_4^{\ast },X_5^{\ast })\mathrm{e}\mathrm{x}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{e},X_i^{\ast }\ne 0,i=1,2,3,4,5$. For ${\mathcal{R}}_f<1,$ the risk of developing KS exists for any ${\mathcal{R}}_0>1.$ The review article by Jeffrey et al. [21] has summarized the circulating levels of infectious agents and the likelihood of infectivity from these levels. We identify in this study the state ${ϵ}_1^{\ast }$ as one of the levels of infectivity of KS.

2. Despite the inequality in Theorem 2.2, that is, ${\mathcal{R}}_0>max\{1,{\mathcal{R}}_f\}$, it is interesting to note that for ${\mathcal{R}}_0={\mathcal{R}}_f$ the components $X_3^{\ast }$ and $X_4^{\ast }$ vanish but the components $X_1^{\ast }$ and $X_2^{\ast }$ are nonzero giving rise to the KS free equilibrium, ${ϵ}_2^{\ast }=(X_1^{\ast },X_2^{\ast },0,0,X_5^{\ast })$. We can calculate the critical value $X_{2f}^{\ast }$ for specified parameter values in (11(11) $\begin{array}{rl}{X}_2^{\ast }& =\frac{{\mu }_{x_1}{x}_{2_{max}}{\mathcal{R}}_f}{K_2},\end{array}$(11) ) below which the HHV-8 viral load is sufficient to maintain the replication of HHV-8 virions only but is not high enough to support the secondary infection of progenitor cells which can lead to the development of KS. The endemic point ${ϵ}_2^{\ast }$ is, however, the starting point for the next KS state $,{ϵ}_3^{\ast },$ discussed below.

3. For ${\mathcal{R}}_0\in [1,\mathrm{\infty })\setminus [1,{\mathcal{R}}_f],X_2\left(t\right)>X_{2f}^{\ast },X_i^{\ast }\ne 0,1,2,3,4,5,$ giving rise to the KS present equilibrium, ${ϵ}_3^{\ast }$. In this case like, in (a), the HHV-8 viral load is sufficient to support both the primary infection of B cells and the secondary infection of progenitor cells making the development of KS real.

4. Note that for $1<{\mathcal{R}}_0<{\mathcal{R}}_f,$ the endemic equilibrium point does not exist by virtue of Theorem 2.1.

5. We conclude that KS does not necessarily develop because ${\mathcal{R}}_0>1,$ but it is sufficient that ${\mathcal{R}}_0>{\mathcal{R}}_f$. (see (b) and (d)).

   We can summarize the results for the innate model as follows:

Lemma 2.2

Consider the system (1(1) $\begin{array}{r}{\dot{X}}_1\left(t\right)=K_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_1\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}{X}_1.\end{array}$(1) )–(5(5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) ). The following statements hold:

1. If ${\mathcal{R}}_0<1,$ then the virus free equilibrium, ${ϵ}^0,$ is the only equilibrium point.

2. If ${\mathcal{R}}_0>{\mathcal{R}}_f,$ then there exist three possible equilibria: the KS-present equilibrium, ${ϵ}_1^{\ast },$ for ${\mathcal{R}}_f<1<{\mathcal{R}}_0,$ the KS-free equilibrium, ${ϵ}_2^{\ast },$ for ${\mathcal{R}}_0={\mathcal{R}}_f$ and the KS-present equilibrium, ${ϵ}_3^{\ast },$ for ${\mathcal{R}}_0\in [1,\mathrm{\infty })\setminus [1,{\mathcal{R}}_f]$.

3. For $1<{\mathcal{R}}_0<{\mathcal{R}}_f,$ no equilibrium point exists by virtue of Theorem 2.1.

3. Numerical simulations of the MIM

The parameter values used in Figure 2 are given in Table A1. Figure 2 shows the sensitivity analysis demonstrating how the model parameters are correlated to the reproduction number, ${\mathcal{R}}_0$. We have found that the rate of infected B-cell proliferation, $K_2$, is positively and significantly correlated with the reproduction number, a conclusion supported by experimental observations by [13]. The other model parameters are not significantly correlated to ${\mathcal{R}}_0$, and their effect on disease progression is peripheral.

Figure 2. PRCCs for parameters of the MIM and $\log \left({\mathcal{R}}_0\right)$ as a function of the most sensitive parameter, $K_2$: (a) scatter plot for ${\mathcal{R}}_0$ and (b) PRCCs for the model.

The parameter values used in Figure 3 are given in Table 1. Figure 3 demonstrates the effect of the parameter ${\theta }_{x_5}$ on disease progression. In particular, we have found a threshold value for ${\theta }_{x_5}$ given by ${\theta }_{x_5}^{\ast }\approx 0.0205$, below which the HHV-8 infection clears even if ${\mathcal{R}}_0>1$. This condition suggests that a potential anti-KS therapy can be found, probably involving pro-inflammatory cytokines such as IL-2 that have already demonstrated the potential to stimulate type I immunity [28,29]. We recommend that experimental and clinical studies should be conducted to assess the therapeutic effects associated with the parameter, ${\theta }_{x_5}$, and quantify its effect in reducing the KS load. We believe that clinical studies are necessary to establish the severity of the infections in (a) and (c) and the possible location of the cancer [19].

Figure 3. Dynamics of the individual components of the MIM for different values of the efficacy threshold ${\theta }_{x_5}$. The HHV-8 infection clears for ${\theta }_{x_5}<{\theta }_{x_5}^{\ast }$ even if ${\mathcal{R}}_0>1$: (a) infected progenitor cells, (b) KS cell dynamics, (c) infected B-cell dynamics, (d) HHV-8 dynamics, (e) Innate immune response.

Table 1. Parameters of the MIM and their definitions.

Variable	Definition	Initial Value	Reference
$X_1$	Infected B cells	0 cell mm${}^{-3}$	Estimated
$X_2$	HHV-8 viral load	0 cell mm${}^{-3}$	Estimated
$X_3$	Infected progenitor cells	0 cell mm${}^{-3}$	[25]
$X_4$	KS cells	0 cell mm${}^{-3}$	Estimated
$X_5$	Innate immune response	${10}^{-2}$	Estimated
Parameter	Definition	Value	Reference
$K_1$	Infected progenitor cell proliferation rate	0.2 day${}^{-1}$	Estimated
$K_2$	Infected B-cell proliferation rate	0.36 day${}^{-1}$	[31]
$K_3$	Innate immune response activation rate	0.01 day${}^{-1}$	[31]
${\theta }_{x_1}$	Innate immune response activation threshold	200 cell mm${}^{-3}$	Estimated
${\theta }_{x_5}$	Efficacy threshold for innate immune response	0.08	Estimated
$x_{1_{max}}$	Infected B-cell carrying capacity	400 cell mm${}^{-3}$	[31]
$x_{2_{max}}$	HHV-8 carrying capacity	$5\times {10}^5$ virions mm${}^{-3}$	[31]
${\mu }_{x_1}$	Death rate of infected B cells	0.33 day${}^{-1}$	[31]
${\mu }_{x_2}$	Clearance rate of HHV-8	0.57 day${}^{-1}$	[31]
${\mu }_{x_3}$	Blanket death rate of infected progenitor cells	0.1 day${}^{-1}$	[31]
${\mu }_{x_4}$	Death rate of KS cells	0.03 day${}^{-1}$	[31]
${\hat{\mu }}_{x_3}$	Progression rate of KS	0.09 day${}^{-1}$	[31]
$N_{x_2}$	Maximum carrying capacity of infected B cells	700	[31]

The MIM gave a very important result regarding the development of KS. First, the model has identified three possible equilibria: ${ϵ}_1^{\ast }$, which exists for ${\mathcal{R}}_f<1<{\mathcal{R}}_0$, the KS-free equilibrium, ${ϵ}_2^{\ast }$, which exists for ${\mathcal{R}}_0={\mathcal{R}}_f$ and ${ϵ}_3^{\ast }$ which exists for ${\mathcal{R}}_0\in [1,\mathrm{\infty })\setminus [1,{\mathcal{R}}_f]$. We recommend clinical studies to establish the severity and location of KS for the two equilibria, ${ϵ}_1^{\ast }$ and ${ϵ}_3^{\ast }$. The equilibrium state, ${ϵ}_3^{\ast }$, is possibly the most common in HHV-8 infected individuals as most of them never develop KS as a result of a high infectivity resistance threshold.

For ${\mathcal{R}}_0={\mathcal{R}}_f,$ the populations of infected progenitor cells and KS cells vanish. For this value of ${\mathcal{R}}_0$, only the primary infection of B cells and replication of HHV-8 take place. Specifically, KS cannot develop, while for ${\mathcal{R}}_f<{\mathcal{R}}_0$ both the production of HHV-8 and KS cells take place.

The MIM will be used to extract the initial conditions for the MAM in the next section for the state variables, $X_i\left(t\right),i=1,2,\dots ,5$ and the parameters ${\theta }_{x_5}$ and ${\mathcal{R}}_f$ for which KS can occur. We want to study the efficacy of the externally administered drugs that can clear/reduce the KS load.

4. MODEL WITH ADAPTIVE MECHANISM (MAM)

4.1. Model formulation and description

In case of HIV-1 and HHV-8 co-infection, a more robust adaptive immune response is developed which includes virus-specific effector cells. Their interaction with the infected and uninfected cell and virus populations is modelled below and depicted in Figure 4. For simplicity of notation, we denote (16) $\begin{array}{r}{x}_i:=x_i\left(t\right),\mathrm{a}\mathrm{n}\mathrm{d}{x}_i^0:=x_i\left(0\right),i=1,2,\dots ,10,\end{array}$(16) where $x_1,x_2,x_3,x_4$ in that order denote the uninfected CD4 T, B-cell populations, HHV-8 and HIV-1-specific effector cells, $x_5,x_6,x_7,x_8$ denote the infected CD4 T cells, B-cell, progenitor cell populations and KS cells, respectively and finally, $x_9,x_{10}$ are the HIV-1 and HHV-8 virions.

We have formulated the adaptive immune response described by the following system of equations: (17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) Equation (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) ) describes the dynamics of the susceptible CD4 T cells, $x_1$. The first term in (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) ) represents the constant natural replacement, ${\mathrm{\Pi }}_1$, of the CD4 T cells, $x_1$, the second term represents proliferation of a proportion of circulating $x_1$ cells due to the presence of HIV-1 virions, at the constant rate, ${\alpha }_5$, and the third term represents proliferation of a proportion of circulating $x_1$ cells due to the presence of HHV-8 virions, at the constant rate, ${\alpha }_6$. There are two proliferation terms because proliferation is pathogen dependent [16,18]. The fourth term is the natural death of these cells at a constant rate, ${\mu }_1$, and the fifth term represents infection of $x_1$ cells by HIV-1 at a constant infection rate ${\text{β}}_1,$ and $S_9$ and $S_{10}$ are half saturation constants of proliferation for HIV-1 and HHV-8, respectively. (18) $\begin{array}{r}\dot{x_2}={\mathrm{\Pi }}_2+\frac{{\mathrm{\Pi }}_2{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_2{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_2{x}_2-{\text{β}}_2{x}_2{x}_{10}.\end{array}$(18) Equation (18(18) $\begin{array}{r}\dot{x_2}={\mathrm{\Pi }}_2+\frac{{\mathrm{\Pi }}_2{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_2{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_2{x}_2-{\text{β}}_2{x}_2{x}_{10}.\end{array}$(18) ) describes the dynamics of the susceptible B cells, $x_2$. The first term represents the constant natural replacement, ${\mathrm{\Pi }}_2$, of the B cells, $x_2$, the second term represents the proliferation of a proportion of circulating $x_2$ cells due to the presence of HIV-1 virions, at the constant rate, ${\alpha }_5$, and the third term represents proliferation of a proportion of circulating $x_2$ cells due to the presence of HHV-8 virions, at the constant rate, ${\alpha }_6$. As in (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) ), there are two pathogen-dependent proliferation terms [16,18]. The fourth term is the natural death of these cells at a constant rate, ${\mu }_2$, and the fifth term represents infection of $x_2$ cells by HHV-8 at a constant infection rate, ${\text{β}}_2$. (19) $\begin{array}{r}{\dot{x}}_3={\mathrm{\Pi }}_3+\frac{{\mathrm{\Pi }}_3{c}_{10}{x}_{10}}{x_{10}+f_{10}}-{\mu }_3{x}_3.\end{array}$(19) Equation (19(19) $\begin{array}{r}{\dot{x}}_3={\mathrm{\Pi }}_3+\frac{{\mathrm{\Pi }}_3{c}_{10}{x}_{10}}{x_{10}+f_{10}}-{\mu }_3{x}_3.\end{array}$(19) ) describes the dynamics of HHV-8 specific effector cells, $x_3$. The first term represents constant replenishment, ${\mathrm{\Pi }}_3$, of these cells from precursors. The second term represents the proliferation of a proportion of circulating $x_3$ cells due to the presence of the HHV-8 virions, at a constant rate, $c_{10}$ with the half saturation constant, $f_{10}$. The last term represents natural death of HHV-8 specific effector cells at the constant death rate, ${\mu }_3$. (20) $\begin{array}{r}{\dot{x}}_4={\mathrm{\Pi }}_4+\frac{{\mathrm{\Pi }}_4{c}_9{x}_9}{x_9+f_9}-{\mu }_4{x}_4.\end{array}$(20) Similar to above, Equation (20(20) $\begin{array}{r}{\dot{x}}_4={\mathrm{\Pi }}_4+\frac{{\mathrm{\Pi }}_4{c}_9{x}_9}{x_9+f_9}-{\mu }_4{x}_4.\end{array}$(20) ) describes the dynamics of HIV-1 specific effector cells, $x_4$. (21) $\begin{array}{r}{\dot{x}}_5=\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}+{\text{β}}_1{x}_1{x}_9-m_4{x}_4{x}_5-{\mu }_5{x}_5.\end{array}$(21) Equation (21(21) $\begin{array}{r}{\dot{x}}_5=\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}+{\text{β}}_1{x}_1{x}_9-m_4{x}_4{x}_5-{\mu }_5{x}_5.\end{array}$(21) ) represents a class of infected CD4 T cells, $x_5$. The first term represents the proliferation of a proportion of circulating,$x_5$, cells due to the presence of HIV-1 virions and the second term represents the proliferation of a proportion of circulating,$x_5$, cells due to the presence of HHV-8 virions. The third term is the gain from infection of T cells by HIV-1. The fourth term is the lysing of infected CD4 T cells by HIV-1 specific effector cells, at a constant term, $m_4$, and the last term represents natural death of these cells at a constant rate, ${\mu }_5$. (22) $\begin{array}{r}{\dot{x}}_6=\frac{{\mathrm{\Pi }}_2{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_2{\alpha }_6{x}_{10}}{x_{10}+S_{10}}+{\text{β}}_2{x}_2{x}_{10}-m_3{x}_3{x}_6-{\mu }_6{x}_6.\end{array}$(22) Similar to above, Equation (22(22) $\begin{array}{r}{\dot{x}}_6=\frac{{\mathrm{\Pi }}_2{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_2{\alpha }_6{x}_{10}}{x_{10}+S_{10}}+{\text{β}}_2{x}_2{x}_{10}-m_3{x}_3{x}_6-{\mu }_6{x}_6.\end{array}$(22) ) represents a class of infected B cells, $x_6$. The first two terms represent proliferation terms, the third is a gain from the infection of B cells and the fourth term represents the killing of these cells by specific effector cells. (23) $\begin{array}{r}{\dot{x}}_7=r_7{x}_{10}\left(1-\frac{x_{10}}{x_{{10}_{max}}}\right)-{\mu }_7{x}_7-d_3{x}_3{x}_7.\end{array}$(23) Equation (23(23) $\begin{array}{r}{\dot{x}}_7=r_7{x}_{10}\left(1-\frac{x_{10}}{x_{{10}_{max}}}\right)-{\mu }_7{x}_7-d_3{x}_3{x}_7.\end{array}$(23) ) represents the dynamics of infected progenitor cells. The first term accounts for the logistic growth rate of these cells that is assumed to depend on HHV-8. The second term is the progression of these cells to KS at a constant rate ${\mu }_7$ [13]. The third term represents the killing of these cells by HHV-8 specific effector cells at a constant rate $d_3$. (24) $\begin{array}{r}{\dot{x}}_8={\mu }_7{x}_7-{\mu }_8{x}_8.\end{array}$(24) Equation (24(24) $\begin{array}{r}{\dot{x}}_8={\mu }_7{x}_7-{\mu }_8{x}_8.\end{array}$(24) ) represents the dynamics of KS. The first term represents the source from infection of progenitor cells. The second term represents natural loss of KS cells. (25) $\begin{array}{r}{\dot{x}}_9=N_9{\mu }_5{x}_5-{\mu }_9{x}_9.\end{array}$(25) Equation (25(25) $\begin{array}{r}{\dot{x}}_9=N_9{\mu }_5{x}_5-{\mu }_9{x}_9.\end{array}$(25) ) represents the dynamics of HIV-1. The first term represents the production of virions from the bursting of the infected CD4 T cells. The parameter $N_9$ represents the maximum carrying capacity of infected CD4 T cells. The second term is the clearance rate of HIV-1. (26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) Equation (26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ) represents the rate of change of HHV-8. The first term represents the rate at which HHV-8 is produced from bursting of the infected B cells. The parameter $N_{10}$ represents the maximum carrying capacity of infected B cells. The last term accounts for the clearance rate of HHV-8.

4.2. Analysis of the model

4.2.1. Positivity and boundedness of solutions

From Equations (25(25) $\begin{array}{r}{\dot{x}}_9=N_9{\mu }_5{x}_5-{\mu }_9{x}_9.\end{array}$(25) )–(26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ), we have $x_i\left(t\right)\ge x_i\left(0\right)exp(-{\mu }_{i}t)\ge 0$, for all $t\ge 0$ and i = 9, 10. To prove the positivity and ensure that the model (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) )–(26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ) is well-posed, we use the following two conditions, (C1) and (C2): $\begin{array}{r}\left(C1\right):f(x_9,x_{10})=\frac{x_9}{x_9+S_9}+\frac{x_{10}}{x_{10}+S_{10}}\ge 0,\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}(x_9,x_{10})\in {\mathbb{R}}_{+}^2.\\ \left(C2\right):g\left(x\right)=rx\left(1-\frac{x}{x_{max}}\right)\ge 0,\mathrm{f}\mathrm{o}\mathrm{r}0\le x\le \frac{x_{max}}{2}.\end{array}$

Lemma 4.1

Consider the system (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) )–(26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ) and assume that $\left(C1\right)$ and $\left(C2\right)$ hold.

1. If $x_i\left(0\right)\ge 0,i=1,2,3,\dots ,10,$ then the solution $x_i\left(t\right)\ge 0,\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}t>0$.

2. Moreover, $x_i\left(t\right)<\mathrm{\infty },i=1,2,\dots ,10,\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}t\ge 0$.

First we will rearrange the system (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) )–(26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ) into a subsystem of uninfected (S1) states (Equations 17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) –20(20) $\begin{array}{r}{\dot{x}}_4={\mathrm{\Pi }}_4+\frac{{\mathrm{\Pi }}_4{c}_9{x}_9}{x_9+f_9}-{\mu }_4{x}_4.\end{array}$(20) ) and infected (S2) states (Equations 21(21) $\begin{array}{r}{\dot{x}}_5=\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}+{\text{β}}_1{x}_1{x}_9-m_4{x}_4{x}_5-{\mu }_5{x}_5.\end{array}$(21) –26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ). It will be shown that if the noninfected states in (S1) are non-negative for all $t\ge 0$, then the infected states in (S2) are non-negative for all $t\ge 0$.

Proof

Proof of Lemma 4.1(a)

The subsystem of uninfected states (S1) can be written as a system of differential inequalities (27) $\begin{array}{r}\frac{\mathrm{d}{x}_i}{\mathrm{d}t}\ge \left(A_i-\sum _{j=1}^6{B}_{ij}{x}_j\right)x_i+{\hat{\mathrm{\Pi }}}_i\end{array}$(27) where $B_{ij}\ge 0,{\hat{\mathrm{\Pi }}}_i=L_i(x_9,x_{10}),$ for $i=1,2,3,4$. The last component is defined as $\hat{\mathrm{\Pi }}(x_9,x_{10})=(L_1,L_2,L_3,L_4{)}^T$, $\begin{array}{rl}{L}_1& ={\mathrm{\Pi }}_1\left(1+\frac{{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\alpha }_6{x}_{10}}{x_{10}+S_{10}}\right),L_3={\mathrm{\Pi }}_3\left(1+\frac{c_{10}{x}_{10}}{x_{10}+f_{10}}\right)\\ L_2& ={\mathrm{\Pi }}_2\left(1+\frac{{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\alpha }_6{x}_{10}}{x_{10}+S_{10}}\right),L_4={\mathrm{\Pi }}_4\left(1+\frac{c_9{x}_9}{x_9+f_9}\right).\end{array}$Clearly, $\hat{\mathrm{\Pi }}(0,0)>(0,0,0,0{)}^T$ by virtue of (A1(A1) $\begin{array}{rl}{\mathcal{R}}_{V_1}& =\frac{N_9{\mathrm{\Pi }}_1\left({\alpha }_5{\mu }_1+{\text{β}}_1{S}_9\right)}{{\mu }_1{\mu }_9{S}_9}.\frac{{\mu }_4{\mu }_5}{({\mu }_4{\mu }_5+m_4{\mathrm{\Pi }}_4)}<\frac{N_9{\mathrm{\Pi }}_1\left({\alpha }_5{\mu }_1+{\text{β}}_1{S}_9\right)}{{\mu }_1{\mu }_9{S}_9}={\mathcal{R}}_1.\end{array}$(A1) ). Suppose the assertion $x_i\left(t\right)\ge 0$ for $i=1,2,3,4$ is not true. Then there exists a smallest number $t_0,$ such that $\begin{array}{rl}{x}_i\left(t\right)& <0\mathrm{f}\mathrm{o}\mathrm{r}1\le i\le 4,0\le t\le t_0\\ x_i\left(t_0\right)& =0\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{o}\mathrm{n}\mathrm{e}\mathrm{i},\mathrm{s}\mathrm{a}\mathrm{y}{i}_0.\end{array}$Then, $x_{i_0}$ is a decreasing function and we would have $\frac{\mathrm{d}{x}_{i_0}\left(t_0\right)}{\mathrm{d}t}\le 0.$However, from the differential inequality (27(27) $\begin{array}{r}\frac{\mathrm{d}{x}_i}{\mathrm{d}t}\ge \left(A_i-\sum _{j=1}^6{B}_{ij}{x}_j\right)x_i+{\hat{\mathrm{\Pi }}}_i\end{array}$(27) ) for $x_{i_0}\left(t\right)$ we get $\frac{\mathrm{d}{x}_{i_0}\left(t_0\right)}{\mathrm{d}t}\ge {\hat{\pi }}_i>0$which is a contradiction. Hence, if $x_i\left(0\right)\ge 0,i=1,2,3,4$, then $x_i\left(t\right)\ge 0$ for all $t>0,i=1,2,3,4$.

The subsystem of infected states (S2) can be written in the matrix form $\dot{Y}\left(t\right)=MY,$ where $Y=[x_5{x}_6{x}_7{x}_8{x}_9{x}_{10}{]}^T$, and $M=\left(\begin{array}{cccccc}-M_{11}& 0& 0& 0& M_{15}& M_{16}\\ 0& -M_{22}& 0& 0& M_{25}& M_{26}\\ 0& 0& -M_{33}& 0& 0& M_{36}\\ 0& 0& {\mu }_7& -{\mu }_8& 0& 0\\ N_9{\mu }_5& 0& 0& 0& -{\mu }_9& 0\\ 0& N_{10}{\mu }_6& 0& 0& 0& -{\mu }_6\end{array}\right)$with entries $\begin{array}{rl}{M}_{11}& ={\mu }_5+m_4{x}_4,M_{22}={\mu }_6+m_3{x}_3,\\ M_{33}& ={\mu }_7+d_3{x}_3,M_{15}=\frac{{\alpha }_5{\mathrm{\Pi }}_1{S}_9}{(x_9+S_9{)}^2}+{\text{β}}_1{x}_1,\\ M_{16}& =\frac{{\alpha }_6{\mathrm{\Pi }}_1{S}_{10}}{(x_{10}+S_{10}{)}^2},M_{25}=\frac{{\alpha }_5{\mathrm{\Pi }}_2{S}_9}{(x_9+S_9{)}^2},\\ M_{26}& =\frac{{\alpha }_6{\mathrm{\Pi }}_2{S}_{10}}{(x_{10}+S_{10}{)}^2}+{\text{β}}_2{x}_2,M_{36}=r_7\left(1-\frac{2x_{10}}{x_{{10}_{max}}}\right).\end{array}$By virtue of (A1(A1) $\begin{array}{rl}{\mathcal{R}}_{V_1}& =\frac{N_9{\mathrm{\Pi }}_1\left({\alpha }_5{\mu }_1+{\text{β}}_1{S}_9\right)}{{\mu }_1{\mu }_9{S}_9}.\frac{{\mu }_4{\mu }_5}{({\mu }_4{\mu }_5+m_4{\mathrm{\Pi }}_4)}<\frac{N_9{\mathrm{\Pi }}_1\left({\alpha }_5{\mu }_1+{\text{β}}_1{S}_9\right)}{{\mu }_1{\mu }_9{S}_9}={\mathcal{R}}_1.\end{array}$(A1) ) and (A2(A2) $\begin{array}{rl}{\mathcal{R}}_{V_8}& =\frac{N_{10}{\mathrm{\Pi }}_2\left({\alpha }_6{\mu }_2+{\text{β}}_2{S}_{10}\right)}{{\mu }_2{\mu }_{10}{S}_{10}}.\frac{{\mu }_3{\mu }_6}{({\mu }_3{\mu }_6+m_3{\mathrm{\Pi }}_3)}<\frac{N_{10}{\mathrm{\Pi }}_2({\alpha }_6{\mu }_2+{\text{β}}_2{S}_{10})}{{\mu }_2{\mu }_{10}{S}_{10}}={\mathcal{R}}_8.\end{array}$(A2) ), M is a Metzler matrix. Hence, the infected states $x_i\left(t\right)\ge 0$ for all t>0, i = 5, 6, 7, 8, 9, 10.

The proof of Lemma 4.1(b) is given in the Appendix.

4.3. Virus free equilibrium and the basic reproduction number

The system (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) )–(26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ) has a virus free equilibrium, ${ϵ}^0$, given by (28) $\begin{array}{r}{ϵ}^0=\left(\frac{{\mathrm{\Pi }}_1}{{\mu }_1},\frac{{\mathrm{\Pi }}_2}{{\mu }_2},\frac{{\mathrm{\Pi }}_3}{{\mu }_3},\frac{{\mathrm{\Pi }}_4}{{\mu }_4},0,0,0,0,0,0\right).\end{array}$(28) Applying the next-generation matrix approach [32], the basic reproduction number for model (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) )–(26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ) reads as follows: (29) $\begin{array}{r}{\mathcal{R}}_0=\frac{1}{2}\left[{\mathcal{R}}_{V_1}+{\mathcal{R}}_{V_8}+\sqrt{{\left({\mathcal{R}}_{V_1}+{\mathcal{R}}_{V_8}\right)}^2+4{\mathcal{R}}_{V_1}{\mathcal{R}}_{V_8}\left(\mathrm{\Phi }-1\right)}\right],\end{array}$(29) where $\mathrm{\Phi }=\frac{{\mu }_1{\mu }_2{\alpha }_5{\alpha }_6}{K_1{K}_2}=\frac{{\mu }_1{\mu }_2{\alpha }_5{\alpha }_6}{{\mu }_1{\mu }_2{\alpha }_5{\alpha }_6+{\mu }_1{\alpha }_5{\text{β}}_2{S}_{10}+{\mu }_2{\alpha }_6{\text{β}}_1{S}_9+{\text{β}}_1{\text{β}}_2{S}_9{S}_{10}},0\le \mathrm{\Phi }<1,$and (30) $\begin{array}{r}{\mathcal{R}}_{V_1}=\frac{{\mu }_4{\mu }_5{N}_9{\mathrm{\Pi }}_1({\alpha }_5{\mu }_1+{\text{β}}_1{S}_9)}{{\mu }_1{\mu }_9{S}_9({\mu }_4{\mu }_5+m_4{\mathrm{\Pi }}_4)},{\mathcal{R}}_{V_8}=\frac{{\mu }_3{\mu }_6{N}_{10}{\mathrm{\Pi }}_2({\alpha }_6{\mu }_2+{\text{β}}_2{S}_{10})}{{\mu }_2{\mu }_{10}{S}_{10}({\mu }_3{\mu }_6+m_3{\mathrm{\Pi }}_3)},\end{array}$(30) are the reproduction numbers attributed to HIV-1 and HHV-8 infections, respectively. Note that the effect of virus-specific effector cells in (30(30) $\begin{array}{r}{\mathcal{R}}_{V_1}=\frac{{\mu }_4{\mu }_5{N}_9{\mathrm{\Pi }}_1({\alpha }_5{\mu }_1+{\text{β}}_1{S}_9)}{{\mu }_1{\mu }_9{S}_9({\mu }_4{\mu }_5+m_4{\mathrm{\Pi }}_4)},{\mathcal{R}}_{V_8}=\frac{{\mu }_3{\mu }_6{N}_{10}{\mathrm{\Pi }}_2({\alpha }_6{\mu }_2+{\text{β}}_2{S}_{10})}{{\mu }_2{\mu }_{10}{S}_{10}({\mu }_3{\mu }_6+m_3{\mathrm{\Pi }}_3)},\end{array}$(30) ) varies from weak to perfect as $m_j$ changes. It readily follows that (31) $\begin{array}{r}{\mathcal{R}}_{V_1}<\frac{N_9{\mathrm{\Pi }}_1({\alpha }_5{\mu }_1+{\text{β}}_1{S}_9)}{{\mu }_1{\mu }_9{S}_9}=:{\mathcal{R}}_1,{\mathcal{R}}_{V_8}<\frac{N_{10}{\mathrm{\Pi }}_2({\alpha }_6{\mu }_2+{\text{β}}_2{S}_{10})}{{\mu }_2{\mu }_{10}{S}_{10}}=:{\mathcal{R}}_8,\end{array}$(31) where ${\mathcal{R}}_1$, ${\mathcal{R}}_8$ are the reproduction numbers when the virus-specific effector cells are dysfunctional.

Given (29(29) $\begin{array}{r}{\mathcal{R}}_0=\frac{1}{2}\left[{\mathcal{R}}_{V_1}+{\mathcal{R}}_{V_8}+\sqrt{{\left({\mathcal{R}}_{V_1}+{\mathcal{R}}_{V_8}\right)}^2+4{\mathcal{R}}_{V_1}{\mathcal{R}}_{V_8}\left(\mathrm{\Phi }-1\right)}\right],\end{array}$(29) )–(30(30) $\begin{array}{r}{\mathcal{R}}_{V_1}=\frac{{\mu }_4{\mu }_5{N}_9{\mathrm{\Pi }}_1({\alpha }_5{\mu }_1+{\text{β}}_1{S}_9)}{{\mu }_1{\mu }_9{S}_9({\mu }_4{\mu }_5+m_4{\mathrm{\Pi }}_4)},{\mathcal{R}}_{V_8}=\frac{{\mu }_3{\mu }_6{N}_{10}{\mathrm{\Pi }}_2({\alpha }_6{\mu }_2+{\text{β}}_2{S}_{10})}{{\mu }_2{\mu }_{10}{S}_{10}({\mu }_3{\mu }_6+m_3{\mathrm{\Pi }}_3)},\end{array}$(30) ) we make the following observations:

Observation 4.1

If ${\alpha }_5=0$ or ${\alpha }_6=0$, then $\mathrm{\Phi }=0$ and the model reproduction number is given by ${\mathcal{R}}_0=max\{{\mathcal{R}}_{V_1},{\mathcal{R}}_{V_8}\}.$

Observation 4.2

Decreasing/increasing ${\alpha }_5$ or ${\alpha }_6$ decreases/increases the reproduction number. Since ${\alpha }_5$ and ${\alpha }_6$ are the proliferation terms of the uninfected/infected B- and T-cell populations, an optimal control approach will be necessary that will balance the level of proliferation during a potential chemo- or immunotherapy.

Observation 4.3

When $\mathrm{\Phi }=1$, the reproduction number reduces to (32) $\begin{array}{r}{\mathcal{R}}_0=\frac{1}{2}\left[{\mathcal{R}}_{V_1}+{\mathcal{R}}_{V_8}+\sqrt{({\mathcal{R}}_{V_1}+{\mathcal{R}}_{V_8}{)}^2}\right]={\mathcal{R}}_{V_1}+{\mathcal{R}}_{V_8}.\end{array}$(32) It is possible in this case for the infection to persist if ${\mathcal{R}}_{V_1}+{\mathcal{R}}_{V_8}>1$, even if both ${\mathcal{R}}_{V_1}<1$ and ${\mathcal{R}}_{V_8}<1\cdot$

Details on the computation of the reproduction number, ${\mathcal{R}}_0,$ and the interpretation of ${\mathcal{R}}_i$ and ${\mathcal{R}}_{V_i}$ for $i=1,8,$ the reader is directed to the appendix.

4.3.1. Global stability for the virus free equilibrium, ${ϵ}^0$

Theorem 4.1

Decompose the system (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) )–(26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ) as in Lemma A.1 in the Appendix. Then the steady state $U_0=(X^{\ast },\mathbf{0})$ of the system (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) )–(26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ) is globally asymptotically stable for ${\alpha }_5={\alpha }_6=0\mathrm{a}\mathrm{n}\mathrm{d}{\mathcal{R}}_0<1$.

Proof.

Denote $X^{\ast }=\left(\frac{{\mathrm{\Pi }}_1}{{\mu }_1},\frac{{\mathrm{\Pi }}_2}{{\mu }_2},\frac{{\mathrm{\Pi }}_3}{{\mu }_3},\frac{{\mathrm{\Pi }}_4}{{\mu }_4}\right)$. Then following [9], we set $X=(x_1,x_2,x_3,x_4)$, $Y=(x_5,x_6,x_7,x_8,x_9,x_{10})$ and define $\begin{array}{rl}F(X,0)& =\left(\begin{array}{c}{\mathrm{\Pi }}_1-{\mu }_1\\ {\mathrm{\Pi }}_2-{\mu }_2\\ {\mathrm{\Pi }}_3-{\mu }_3\\ {\mathrm{\Pi }}_4-{\mu }_4\end{array}\right),\\ G(X,Y)& =\left(\begin{array}{c}{m}_4(x_4^{\ast }-x_4^0)x_5+{\text{β}}_1(x_1^0-x_1^{\ast })x_9+{\mathrm{\Pi }}_1{\mathrm{\Sigma }}_{j=5}^6{\delta }_j{x}_{j+4}\\ m_4(x_3^{\ast }-x_3^0)x_6+{\text{β}}_2(x_2^0-x_2^{\ast })x_{10}+{\mathrm{\Pi }}_2{\mathrm{\Sigma }}_{j=5}^6{\delta }_j{x}_{j+4}\\ d_3(x_3^{\ast }-x_3^0)x_7+r_7\frac{x_{10}^2}{x_{{10}_{max}}}\\ 0\\ 0\\ 0\end{array}\right),\end{array}$where ${\delta }_j={\alpha }_j\left(\frac{1}{S_{j+4}}-\frac{1}{x_{j+4}+S_{j+4}}\right),j=5,6.$The global stability of the system (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) )–(26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ) at ${ϵ}^0$ requires that $\hat{G}(X,Y)\ge \mathbf{0}\left[13\right].\mathrm{M}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r},x_i^{\ast }>x_i^0,\mathrm{f}\mathrm{o}\mathrm{r}i=3,4$ and $x_j^0>x_j^{\ast },\mathrm{f}\mathrm{o}\mathrm{r}j=1,2.$ Then $X^{\ast }$ is a globally asymptotically stable solution of the system $\frac{\mathrm{d}X}{\mathrm{d}t}=F(X,0)$ since $F(X,0)$ is the limiting function of $\frac{\mathrm{d}X}{\mathrm{d}t}=F\left(X\right(t),Y(t\left)\right)$, that is, $\underset{t\to \mathrm{\infty }}{lim}X\left(t\right)=X^{\ast }$. It follows that $\hat{G}(X,Y)\ge \mathbf{0}$ and so ${ϵ}^0$ is globally asymptotically stable.

4.4. AIDS-KS-present equilibrium, $E^{\ast }$

Theorem 4.2

Consider the system (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) )–(26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ). The KS present equilibrium, $E^{\ast }$, exists if $0<x_{10}^{\ast }<x_{10max}$ and $x_9^{\ast }>0$

where (33) $\begin{array}{r}\left\{\begin{array}{l}{x}_1^{\ast }=\frac{{\mathrm{\Pi }}_1(S_9+x_9^{\ast }(1+{\alpha }_5+2{\alpha }_6))x_{10}^{\ast }+{\mathrm{\Pi }}_1{S}_{10}(S_9+x_9^{\ast }(1+{\alpha }_5))}{(x_9^{\ast }+S_9)(x_{10}^{\ast }+S_{10})({\mu }_1+{\text{β}}_1{x}_9^{\ast })},\\ x_2^{\ast }=\frac{{\mathrm{\Pi }}_2(S_9(1+{\alpha }_6)+x_9^{\ast }(1+{\alpha }_5+{\alpha }_6))x_{10}^{\ast }+{\mathrm{\Pi }}_2{S}_{10}(S_9+x_9^{\ast }(1+{\alpha }_5))}{(x_9^{\ast }+S_9)(x_{10}^{\ast }+S_{10})({\mu }_2+{\text{β}}_2{x}_{10}^{\ast })}\\ x_3^{\ast }=\frac{{\mathrm{\Pi }}_3(f_{10}+x_{10}^{\ast }(1+c_{10}))}{{\mu }_3(x_{10}^{\ast }+f_{10})},\\ x_4^{\ast }=\frac{{\mathrm{\Pi }}_4(f_9+x_9^{\ast }(1+c_9))}{{\mu }_4(x_9^{\ast }+f_9)},\\ x_5^{\ast }=\frac{{\mathrm{\Pi }}_1{\zeta }_1(x_9^{\ast },x_{10}^{\ast })+{\text{β}}_1{x}_1^{\ast }{x}_9^{\ast }}{{\mu }_5+m_4{x}_4^{\ast }},x_6^{\ast }=\frac{{\mathrm{\Pi }}_1{\zeta }_1(x_9^{\ast },x_{10}^{\ast })+{\text{β}}_2{x}_2^{\ast }{x}_{10}^{\ast }}{{\mu }_6+m_3{x}_3^{\ast }},\\ x_7^{\ast }=\frac{{\mu }_3{r}_7{x}_{10}^{\ast }}{({\mu }_3{\mu }_7+d_3{\mathrm{\Pi }}_3(1+{\zeta }_2\left(x_{10}^{\ast }\right)))}(1-\frac{x_{10}^{\ast }}{x_{10max}}),\\ x_8^{\ast }=\frac{{\mu }_7}{{\mu }_8}{x}_7^{\ast },{\zeta }_1(x_9^{\ast },x_{10}^{\ast })=\frac{{\alpha }_5{x}_9^{\ast }}{x_9^{\ast }+S_9}+\frac{{\alpha }_6{x}_{10}^{\ast }}{x_{10}^{\ast }},{\zeta }_2(x_{10}^{\ast })=\frac{c_{10}{x}_{10}^{\ast }}{x_{10}^{\ast }+f_{10}}.\end{array}\right.\end{array}$(33)

where $x_9^{\ast }$ and $x_{10}^{\ast }$ are proved to be positive solutions of the fourth degree polynomials $Q_i,i=1,2$. The reader is directed to the appendix for details.

Remark 4.1

Observe that when $x_9^{\ast }=0$ and $x_{10}^{\ast }=0$, we obtain the virus free equilibrium, ${ϵ}^0$ in Equation (28(28) $\begin{array}{r}{ϵ}^0=\left(\frac{{\mathrm{\Pi }}_1}{{\mu }_1},\frac{{\mathrm{\Pi }}_2}{{\mu }_2},\frac{{\mathrm{\Pi }}_3}{{\mu }_3},\frac{{\mathrm{\Pi }}_4}{{\mu }_4},0,0,0,0,0,0\right).\end{array}$(28) ).

4.5. Non-AIDS-KS-present equilibrium, $E^{\ast \ast }$

Theorem 4.3

Consider the system (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) )–(26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ). The non-AIDS-KS-present equilibrium, $E^{\ast \ast }$, exists if $0<x_{10}^{\ast \ast }<x_{10max}$ and ${\mathcal{R}}_{V_8}>1$

where (34) $\begin{array}{r}\left\{\begin{array}{l}{x}_1^{\ast \ast }=\frac{{\mathrm{\Pi }}_1}{{\mu }_1}(1+\frac{{\alpha }_6{x}_{10}^{\ast \ast }}{x_{10}^{\ast \ast }+S_{10}}),x_2^{\ast \ast }=\frac{{\mathrm{\Pi }}_2}{({\mu }_2+{\text{β}}_2{x}_{10}^{\ast \ast })}(1+\frac{{\alpha }_6{x}_{10}^{\ast \ast }}{x_{10}^{\ast \ast }+S_{10}}),x_3^{\ast \ast }=\frac{{\mathrm{\Pi }}_3}{{\mu }_3}(1+\frac{c_{10}{x}_{10}^{\ast \ast }}{x_{10}^{\ast \ast }+S_{10}}),\\ x_4^{\ast \ast }=\frac{{\mathrm{\Pi }}_4}{{\mu }_4},x_5^{\ast \ast }=0,x_6^{\ast \ast }=\frac{1}{({\mu }_6+m_3{x}_3^{\ast \ast })}({\text{β}}_2{x}_2^{\ast \ast }{x}_{10}^{\ast \ast }+\frac{{\mathrm{\Pi }}_2{\alpha }_6{x}_{10}^{\ast \ast }}{x_{10}^{\ast \ast }+S_{10}})\\ x_7^{\ast \ast }=\frac{r_7{x}_{10}^{\ast \ast }}{({\mu }_7+d_3{x}_3^{\ast \ast })}(1-\frac{x_{10}^{\ast \ast }}{x_{10max}}),x_8^{\ast }=\frac{{\mu }_7{r}_7{x}_{10}^{\ast \ast }}{{\mu }_8({\mu }_7+d_3{x}_3^{\ast \ast })}(1-\frac{x_{10}^{\ast \ast }}{x_{10max}}),x_9^{\ast \ast }=0.\end{array}\right.\end{array}$(34) where $x_{10}^{\ast \ast }$ is a positive zero to (35) $\begin{array}{rl}& G(x_{10})=D_3{x}_{10}^3+D_2{x}_{10}^2+D_1{x}_{10}+D_0,\end{array}$(35) (36) $\begin{array}{rl}& \left\{\begin{array}{l}{D}_3=-{\text{β}}_2{\mu }_{10}({\mu }_3{\mu }_6+m_3{\mathrm{\Pi }}_3(1+c_{10}))<0,\\ D_2=N_{10}{\mathrm{\Pi }}_2{\text{β}}_2{\mu }_3{\mu }_6-{\mathrm{\Pi }}_3{m}_3{\mu }_2{\mu }_{10}-{\mu }_2{\mu }_3{\mu }_6{\mu }_{10}-{\mathrm{\Pi }}_3{S}_{10}{\text{β}}_2{m}_3{\mu }_{10}\\ -{\mathrm{\Pi }}_3{\text{β}}_2{f}_{10}{m}_3{\mu }_{10}-{\mathrm{\Pi }}_3{c}_{10}{m}_3{\mu }_2{\mu }_{10}-S_{10}{\text{β}}_2{\mu }_3{\mu }_6{\mu }_{10}\\ -{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{\mu }_{10}+2N_{10}{\mathrm{\Pi }}_2{\alpha }_6{\text{β}}_2{\mu }_3{\mu }_6-{\mathrm{\Pi }}_3{S}_{10}{\text{β}}_2{c}_{10}{m}_3{\mu }_{10},\\ D_1=N_{10}{\mathrm{\Pi }}_2{S}_{10}{\text{β}}_2{\mu }_3{\mu }_6-{\mathrm{\Pi }}_3{f}_{10}{m}_3{\mu }_2{\mu }_{10}-S_{10}{\mu }_2{\mu }_3{\mu }_6{\mu }_{10}\\ -f_{10}{\mu }_2{\mu }_3{\mu }_6{\mu }_{10}-{\mathrm{\Pi }}_3{S}_{10}{m}_3{\mu }_2{\mu }_{10}\\ +N_{10}{\mathrm{\Pi }}_2{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6-{\mathrm{\Pi }}_3{S}_{10}{\text{β}}_2{f}_{10}{m}_3{\mu }_{10}+N_{10}{\mathrm{\Pi }}_2{\mu }_2{\mu }_3{\mu }_6-{\mathrm{\Pi }}_3{S}_{10}{c}_{10}{m}_3{\mu }_2{\mu }_{10}\\ -S_{10}{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{\mu }_{10}+2N_{10}{\mathrm{\Pi }}_2{\alpha }_6{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6.\\ D_0=S_{10}{f}_{10}{\mu }_2{\mu }_{10}({\mu }_3{\mu }_6+m_3{\mathrm{\Pi }}_3)({\mathcal{R}}_{V_8}-1)>0,\text{if and only if }{\mathcal{R}}_{V_8}>1.\end{array}\right.\end{array}$(36)

Note that $G\left(0\right)=D_0>0$ if and only if ${\mathcal{R}}_{V_8}>1$. Using the continuity property of cubic polynomials, we have $\underset{x_{10}\to \mathrm{\infty }}{lim}G\left(x_{10}\right)=-\mathrm{\infty }.$ Hence, there must be a positive value of $x_{10},$ namely, $x_{10}^{\ast \ast }\in (0,\mathrm{\infty })$ such that $G\left(x_{10}^{\ast \ast }\right)=0.$

5. Numerical simulations of the MAM

Figure 5 shows the sensitivity analysis of parameters from Table A2 on ${\mathcal{R}}_0$. The maximum carrying capacity of infected B cells is strongly positively correlated with ${\mathcal{R}}_0$, whereas the HHV-8 clearance rate has the opposite effect. The logarithm of the reproduction number as a function of these two parameters is also shown. For small values of $N_{10}$, the increase is fast, implying that one can overestimate the severity of the infection. For large values of $N_{10}$, the increase is slow, implying that disease progression is stable and the conclusions are not adversely affected. The relationship between the $\log \left({\mathcal{R}}_0\right)$ and the clearance rate of HHV-8 virions is almost linear, suggesting a moderate negative correlation with the progression of AIDS-KS.

Figure 4. Schematic diagram for the MAM.

Figure 5. PRCCs for parameters of the MAM and log(${\mathcal{R}}_0$) as a function of the two most sensitive parameters, $N_{10}$ and ${\text{β}}_2$: (a) scatter plot for ${\mathcal{R}}_0$, (b) scatter plot for ${\mathcal{R}}_0$ and (c) PRCCs for the model.

The parameter values used in Figures 6–8 are given in Table 2. Assuming an advanced HIV-1 and HHV-8 co-infection stage, the initial condition of the system (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) )–(26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ) was chosen to be the equilibrium point of the system (1(1) $\begin{array}{r}{\dot{X}}_1\left(t\right)=K_2\left(1-\frac{X_5}{{\theta }_{x_5}+X_5}\right)X_1\left(1-\frac{X_1}{x_{1_{max}}}\right)-{\mu }_{x_1}{X}_1.\end{array}$(1) )–(5(5) $\begin{array}{r}{\dot{X}}_5\left(t\right)=K_3\frac{X_1}{{\theta }_{x_1}+X_1}-K_3{X}_5\end{array}$(5) ). Figure 6 shows the dynamics of the uninfected/infected B cells and CD4 T cells, infected progenitor cells, KS cells and virions for a period of 300 days. The population of infected CD4 T cells reaches a peak at 100 days by then the population of infected B cells have already reached a steady state. These results, in line with experimental evidence, support the fact that KS may accelerate the clinical course of HIV-1 infection.

Figure 6. Dynamics of the individual components of the MAM for ${\mathcal{R}}_0=1.7902$: (a) uninfected CD4 T cells and infected CD4 T-cells; (b) uninfected B cells and infected B cells; (c) infected progenitor cells and KS cells, and (d) HIV-1 and HHV-8.

Figure 7. Comparison of the long-term dynamics of the uninfected CD4 T cells and B cells with regard to HIV-1 and HHV-8 load: (a) HIV-1 and uninfected CD4 T cells and (b) HHV-8 and uninfected B cells.

Figure 8. (a) Same as Figure 7(a) but for the first 300 days. (b) The long-term dynamics of infected B cells and KS cells with regards to HHV-8 load: (a) HIV-1 and uninfected CD4 T cells and (b) KS cells, HHV-8 and infected B-cell dynamics.

Table 2. Parameters of the MAM and their definitions.

Variable	Definitions	Initial value	Reference
$x_1$	Uninfected CD4 T cells	${10}^5$ cell ml${}^{-1}$	[25]
$x_2$	Uninfected B cells	${10}^4$ cell ml${}^{-1}$	Estimated
$x_3$	HHV-8 effector cells	cell ml${}^{-1}$	MIM
$x_4$	HIV-1 effector cells	${10}^4$ cell ml${}^{-1}$	Estimated
$x_5$	Infected CD4 T cells	300 cell ml${}^{-1}$	Estimated
$x_6$	Infected B cells	32.4 cell ml${}^{-1}$	MIM
$x_7$	Infected progenitor cells	61.3 cell ml${}^{-1}$	MIM
$x_8$	KS cells	30.1 cell ml${}^{-1}$	MIM
$x_9$	HHV-1 virions	100 virions ml${}^{-1}$	Estimated
$x_{10}$	HHV-8 virions	83.6 cell ml${}^{-1}$	MIM
Parameter	Definitions	Value	Reference
${\alpha }_5$	HIV-1-dependent enhancement of CD4 T-cell proliferation	0.013	Estimated
${\alpha }_6$	HHV-8-dependent enhancement of CD4 T-cell proliferation	0.045	Estimated
${\text{β}}_1$	Infection rate of CD4 T cells by HIV-1	$2.4\times {10}^{-8}$ ml virion${}^{-1}$ day${}^{-1}$	[31]
${\text{β}}_2$	Infection rate of B cells by HHV-8	$1.5\times {10}^{-7}$ ml virion${}^{-1}$ day${}^{-1}$	[31]
${\mathrm{\Pi }}_1$	Source of new CD4 T cells from the thymus	${10}^4$ cell ml${}^{-1}$ day${}^{-1}$	[31]
${\mathrm{\Pi }}_2$	Source of new B cells from the bone marrow	48000 cell ml${}^{-1}$ day${}^{-1}$	[31]
${\mathrm{\Pi }}_3$	Constant source of HHV-8 specific effector cells	$5\times {10}^3$ cell ml${}^{-1}$ day${}^{-1}$	Estimated
${\mathrm{\Pi }}_4$	Constant source of HIV-1 specific effector cells	$2\times {10}^4$ cell ml${}^{-1}$ day${}^{-1}$	[31]
$S_9$	Half saturation for HIV-1-stimulated CD4 T-cell proliferation	$3\times {10}^5$ virions ml${}^{-1}$	Estimated
$S_{10}$	Half saturation for HHV-8-stimulated CD4-T cell proliferation	$2\times {10}^5$ virions ml${}^{-1}$	Estimated
${\mu }_1$	Death rate of CD4 T cells	0.01 day${}^{-1}$	[31]
${\mu }_2$	Death rate of B cells	0.24 day${}^{-1}$	[31]
${\mu }_3$	HHV-8 specific effector cell death rate	0.1 day${}^{-1}$	[31]
${\mu }_4$	HIV-1 specific effector cell death rate	0.1 day${}^{-1}$	[31]
${\mu }_5$	Lytic death rate of infected CD4 T cells	0.24 day${}^{-1}$	[31]
${\mu }_6$	Death rate of free HHV-8 virus	0.33 day${}^{-1}$	[31]
${\mu }_7$	KS progression rate	0.1 day${}^{-1}$	[31]
${\mu }_8$	Natural death rate of KS cells	0.21 day${}^{-1}$	[31]
${\mu }_9$	Clearance rate of HIV-1 virions	3 day${}^{-1}$	[31]
${\mu }_{10}$	clearance rate of HHV-8 virions	0.57 day${}^{-1}$	[31]
$N_9$	Maximum carrying capacity of infected CD4 T cells	1000 virions cell${}^{-1}$	[31]
$N_{10}$	Maximum carrying capacity of infected B cells	700 virions cell${}^{-1}$	[31]
$m_3$	Killing rate of HHV-8 infected B cells	$1.08\times {10}^{-4}$ ml cell${}^{-1}$ day${}^{-1}$	Estimated
$m_4$	Killing rate of HIV-1 infected CD4 T cells	$2\times {10}^{-7}$ ml cell${}^{-1}$ day${}^{-1}$	Estimated
$c_9$	Proliferation rate of HIV-1 specific CTLs due	0.03	[31]
$c_{10}$	Proliferation rate of HHV-8 specific CTLs	0.047	[31]
$d_3$	Killing rate of infected progenitor cells	$5\times {10}^{-4}$ ml cell${}^{-1}$ day${}^{-1}$	Estimated
$r_7$	Maximum proliferation rate of infected progenitor cells	0.33 cells virion${}^{-1}$ day${}^{-1}$	[31]
$f_9$	Half saturation for proliferation of HIV-1 specific CTLs	$3\times {10}^5$ virions ml${}^{-1}$	Estimated
$f_{10}$	Half saturation for proliferation of HHV-8 specific CTLs	$2\times {10}^5$ virions ml${}^{-1}$	Estimated
$x_{{10}_{max}}$	Maximum HHV-8 load	$8\times {10}^9$ virions ml${}^{-1}$	Estimated

Figure 7 compares the long-term dynamics of HIV-1 and HHV-8 populations with uninfected CD4 T and B cells. The uninfected CD4 T-cell population starts declining rapidly after about 300 days and the HIV-1 population switches from a stable to an exponential growth after about 6 years. For the infected B-cell population, this switch occurs after about 8 years. The delay in the switching time from a stable to exponential growth between the HIV-1 and HHV-8 viral load indicates that HHV-1 supports the clinical course of KS. Figure 8 shows that the uninfected CD4 T-cell population peaks at about 60 days, 40 days before the HIV-1 peak and that the KS cell population will have a steeper rise than the relatively stable infected B-cell population. This suggests that at later stages of AIDS-KS, the reservoirs of HHV-8 will be predominantly the infected progenitor cells.

6. Optimal control applied to the MAM

We now apply an optimal control approach to the system (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) )–(26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ). In order to determine the optimal strategy for controlling AIDS-KS with cART, we introduce three time-dependent controls: $u_1\left(t\right),u_2\left(t\right)$ and $u_3\left(t\right)$. The control $u_1\left(t\right)$ is the efficacy of HAART for preventing the infection of CD4 T cells by HIV-1 and $u_2\left(t\right)$ is the efficacy of HAART in preventing the infection of B cells by HHV-8. The third control $u_3\left(t\right)$ represents the efficacy of an anti-KS therapy by enhancing the proliferation of CD4 T and B cells. (37) $\begin{array}{rl}\dot{x_1}& ={\mathrm{\Pi }}_1+\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)x_9}{x_9+S_9}{\mathrm{\Pi }}_1+{\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x}}_{10}\\ & +S_{10}{\mathrm{\Pi }}_1-{\mu }_1{x}_1-\left(1-u_1\left(t\right)\right){\text{β}}_1{x}_1{x}_9.\end{array}$(37) Equation (37(37) $\begin{array}{rl}\dot{x_1}& ={\mathrm{\Pi }}_1+\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)x_9}{x_9+S_9}{\mathrm{\Pi }}_1+{\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x}}_{10}\\ & +S_{10}{\mathrm{\Pi }}_1-{\mu }_1{x}_1-\left(1-u_1\left(t\right)\right){\text{β}}_1{x}_1{x}_9.\end{array}$(37) ) describes the dynamics of the uninfected CD4 T cells. The meanings of the various terms are given in Equation (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) ). However, the proliferation constants, ${\alpha }_i,i=5,6$ in (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) ) are now replaced by $\left({\alpha }_5\right(1-u_3\left)\right)$ and $\left(ϵ{\alpha }_5\right(1-u_3\left)\right)$, respectively. Note that ${\alpha }_6$ is here assumed to be a multiple of ${\alpha }_5$, with ϵ being the constant of proportionality, which is purely done for technical reasons. The infection coefficient, ${\text{β}}_1$, is replaced by $(1-u_1(t\left)\right){\text{β}}_1,$ where $u_1\left(t\right)$ is the efficacy of HAART in preventing the infection of healthy CD4 T cells and hence, $u_1\left(t\right)$ could represent treatment with either the fusion inhibitor or the reverse transcriptase inhibitor. (38) $\begin{array}{rl}\dot{x_2}& ={\mathrm{\Pi }}_2+\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)5x_9}{x_9+S_9}{\mathrm{\Pi }}_2\\ & +\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x_{10}+S_{10}}{\mathrm{\Pi }}_2-{\mu }_2{x}_2-\left(1-u_1\left(t\right)\right){\text{β}}_2{x}_2{x}_{10}.\end{array}$(38) Equation (38(38) $\begin{array}{rl}\dot{x_2}& ={\mathrm{\Pi }}_2+\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)5x_9}{x_9+S_9}{\mathrm{\Pi }}_2\\ & +\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x_{10}+S_{10}}{\mathrm{\Pi }}_2-{\mu }_2{x}_2-\left(1-u_1\left(t\right)\right){\text{β}}_2{x}_2{x}_{10}.\end{array}$(38) ) describes the dynamics of the susceptible B cells. The meanings of the various terms are given in Equation (18(18) $\begin{array}{r}\dot{x_2}={\mathrm{\Pi }}_2+\frac{{\mathrm{\Pi }}_2{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_2{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_2{x}_2-{\text{β}}_2{x}_2{x}_{10}.\end{array}$(18) ). The second and third terms are explained in E (37(37) $\begin{array}{rl}\dot{x_1}& ={\mathrm{\Pi }}_1+\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)x_9}{x_9+S_9}{\mathrm{\Pi }}_1+{\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x}}_{10}\\ & +S_{10}{\mathrm{\Pi }}_1-{\mu }_1{x}_1-\left(1-u_1\left(t\right)\right){\text{β}}_1{x}_1{x}_9.\end{array}$(37) ). The first and fourth terms are explained in Equation (18(18) $\begin{array}{r}\dot{x_2}={\mathrm{\Pi }}_2+\frac{{\mathrm{\Pi }}_2{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_2{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_2{x}_2-{\text{β}}_2{x}_2{x}_{10}.\end{array}$(18) ). The infection coefficient, ${\text{β}}_2$, is replaced by $(1-u_1(t\left)\right){\text{β}}_2,$ where $u_1\left(t\right)$ is the efficacy of HAART in preventing the infection of healthy B cells. It is well documented that AIDS-KS patients undergoing HAART treatment undergo remission and hence, we assume that HAART treatment reduces the infection rate of B cells by HHV-8. (39) $\begin{array}{r}{\dot{x}}_3={\mathrm{\Pi }}_3+\frac{c_{10}{\mathrm{\Pi }}_3{x}_{10}}{x_{10}+f_{10}}-{\mu }_3{x}_3.\end{array}$(39) Equation (39(39) $\begin{array}{r}{\dot{x}}_3={\mathrm{\Pi }}_3+\frac{c_{10}{\mathrm{\Pi }}_3{x}_{10}}{x_{10}+f_{10}}-{\mu }_3{x}_3.\end{array}$(39) ) describes the dynamics of HHV-8 specific effector cells, $x_3$. The terms are as explained in Equation (19(19) $\begin{array}{r}{\dot{x}}_3={\mathrm{\Pi }}_3+\frac{{\mathrm{\Pi }}_3{c}_{10}{x}_{10}}{x_{10}+f_{10}}-{\mu }_3{x}_3.\end{array}$(19) ). (40) $\begin{array}{r}{\dot{x}}_4={\mathrm{\Pi }}_4+\frac{c_9{\mathrm{\Pi }}_4{x}_9}{x_9+f_9}-{\mu }_4{x}_4.\end{array}$(40) Equation (40(40) $\begin{array}{r}{\dot{x}}_4={\mathrm{\Pi }}_4+\frac{c_9{\mathrm{\Pi }}_4{x}_9}{x_9+f_9}-{\mu }_4{x}_4.\end{array}$(40) ) describes the dynamics of HIV-1 specific effector cells, $x_4$.The terms are as explained in Equation (20(20) $\begin{array}{r}{\dot{x}}_4={\mathrm{\Pi }}_4+\frac{{\mathrm{\Pi }}_4{c}_9{x}_9}{x_9+f_9}-{\mu }_4{x}_4.\end{array}$(20) ). (41) $\begin{array}{rl}{\dot{x}}_5& =\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)x_9}{x_9+S_9}{\mathrm{\Pi }}_1+\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x_{10}+S_{10}}{\mathrm{\Pi }}_1\\ & +\left(1-u_1\left(t\right)\right){\text{β}}_1{x}_1{x}_9-m_4{x}_4{x}_5-{\mu }_5{x}_5.\end{array}$(41) Equation (41(41) $\begin{array}{rl}{\dot{x}}_5& =\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)x_9}{x_9+S_9}{\mathrm{\Pi }}_1+\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x_{10}+S_{10}}{\mathrm{\Pi }}_1\\ & +\left(1-u_1\left(t\right)\right){\text{β}}_1{x}_1{x}_9-m_4{x}_4{x}_5-{\mu }_5{x}_5.\end{array}$(41) ) represents a class of infected CD4 T cells, $x_5$. The first two terms are as explained in (37(37) $\begin{array}{rl}\dot{x_1}& ={\mathrm{\Pi }}_1+\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)x_9}{x_9+S_9}{\mathrm{\Pi }}_1+{\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x}}_{10}\\ & +S_{10}{\mathrm{\Pi }}_1-{\mu }_1{x}_1-\left(1-u_1\left(t\right)\right){\text{β}}_1{x}_1{x}_9.\end{array}$(37) ) and (38(38) $\begin{array}{rl}\dot{x_2}& ={\mathrm{\Pi }}_2+\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)5x_9}{x_9+S_9}{\mathrm{\Pi }}_2\\ & +\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x_{10}+S_{10}}{\mathrm{\Pi }}_2-{\mu }_2{x}_2-\left(1-u_1\left(t\right)\right){\text{β}}_2{x}_2{x}_{10}.\end{array}$(38) ) above. The third term is the gain from Equation (37(37) $\begin{array}{rl}\dot{x_1}& ={\mathrm{\Pi }}_1+\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)x_9}{x_9+S_9}{\mathrm{\Pi }}_1+{\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x}}_{10}\\ & +S_{10}{\mathrm{\Pi }}_1-{\mu }_1{x}_1-\left(1-u_1\left(t\right)\right){\text{β}}_1{x}_1{x}_9.\end{array}$(37) ) and the other terms are as explained in Equation (21(21) $\begin{array}{r}{\dot{x}}_5=\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}+{\text{β}}_1{x}_1{x}_9-m_4{x}_4{x}_5-{\mu }_5{x}_5.\end{array}$(21) ). (42) $\begin{array}{rl}{\dot{x}}_6& =\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)x_9}{x_9+S_9}{\mathrm{\Pi }}_2+\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x_{10}+S_{10}}{\mathrm{\Pi }}_2\\ & +\left(1-u_1\left(t\right)\right){\text{β}}_2{x}_2{x}_{10}-m_3{x}_3{x}_6-{\mu }_6{x}_6.\end{array}$(42) Equation (42(42) $\begin{array}{rl}{\dot{x}}_6& =\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)x_9}{x_9+S_9}{\mathrm{\Pi }}_2+\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x_{10}+S_{10}}{\mathrm{\Pi }}_2\\ & +\left(1-u_1\left(t\right)\right){\text{β}}_2{x}_2{x}_{10}-m_3{x}_3{x}_6-{\mu }_6{x}_6.\end{array}$(42) ) represents a class of infected B cells, $x_6$. The first two terms are as explained in (41(41) $\begin{array}{rl}{\dot{x}}_5& =\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)x_9}{x_9+S_9}{\mathrm{\Pi }}_1+\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x_{10}+S_{10}}{\mathrm{\Pi }}_1\\ & +\left(1-u_1\left(t\right)\right){\text{β}}_1{x}_1{x}_9-m_4{x}_4{x}_5-{\mu }_5{x}_5.\end{array}$(41) ) above. The third term is the gain from Equation (38(38) $\begin{array}{rl}\dot{x_2}& ={\mathrm{\Pi }}_2+\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)5x_9}{x_9+S_9}{\mathrm{\Pi }}_2\\ & +\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x_{10}+S_{10}}{\mathrm{\Pi }}_2-{\mu }_2{x}_2-\left(1-u_1\left(t\right)\right){\text{β}}_2{x}_2{x}_{10}.\end{array}$(38) ) and the other terms are as explained in Equation (22(22) $\begin{array}{r}{\dot{x}}_6=\frac{{\mathrm{\Pi }}_2{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_2{\alpha }_6{x}_{10}}{x_{10}+S_{10}}+{\text{β}}_2{x}_2{x}_{10}-m_3{x}_3{x}_6-{\mu }_6{x}_6.\end{array}$(22) ). (43) $\begin{array}{r}{\dot{x}}_7=\left(1-u_1\left(t\right)\right)r_7{x}_{10}\left(1-\frac{x_{10}}{x_{{10}_{max}}}\right)-{\mu }_7{x}_7-d_3{x}_3{x}_7.\end{array}$(43) Equation (43(43) $\begin{array}{r}{\dot{x}}_7=\left(1-u_1\left(t\right)\right)r_7{x}_{10}\left(1-\frac{x_{10}}{x_{{10}_{max}}}\right)-{\mu }_7{x}_7-d_3{x}_3{x}_7.\end{array}$(43) ) represents the dynamics of infected progenitor cells. In the first term, the logistic growth rate $r_7$ is reduced by a factor $(1-u_1(t\left)\right)$ due to the action of HAART in blocking the infection of progenitor cells by HHV-8. The other terms are as explained before in Equation (23(23) $\begin{array}{r}{\dot{x}}_7=r_7{x}_{10}\left(1-\frac{x_{10}}{x_{{10}_{max}}}\right)-{\mu }_7{x}_7-d_3{x}_3{x}_7.\end{array}$(23) ). (44) $\begin{array}{r}{\dot{x}}_8={\mu }_7{x}_7-{\mu }_8{x}_8.\end{array}$(44) Equation (44(44) $\begin{array}{r}{\dot{x}}_8={\mu }_7{x}_7-{\mu }_8{x}_8.\end{array}$(44) ) represents the dynamics of KS. The terms are already explained in Equation (24(24) $\begin{array}{r}{\dot{x}}_8={\mu }_7{x}_7-{\mu }_8{x}_8.\end{array}$(24) ). (45) $\begin{array}{r}{\dot{x}}_9=N_9{\mu }_5\left(1-u_2\left(t\right)\right)x_5-{\mu }_9{x}_9.\end{array}$(45) Equation (45(45) $\begin{array}{r}{\dot{x}}_9=N_9{\mu }_5\left(1-u_2\left(t\right)\right)x_5-{\mu }_9{x}_9.\end{array}$(45) ) represents the HIV-1 dynamics. The first term is decreased by a factor $1-u_2\left(t\right)$ to reflect the action of HAART in blocking the production of infectious and mature HIV-1 virions. The other term is explained in Equation (25(25) $\begin{array}{r}{\dot{x}}_9=N_9{\mu }_5{x}_5-{\mu }_9{x}_9.\end{array}$(25) ). (46) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6\left(1-u_2\left(t\right)\right)x_6-{\mu }_{10}{x}_{10}.\end{array}$(46) Similar to Equation (45(45) $\begin{array}{r}{\dot{x}}_9=N_9{\mu }_5\left(1-u_2\left(t\right)\right)x_5-{\mu }_9{x}_9.\end{array}$(45) ), Equation (46(46) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6\left(1-u_2\left(t\right)\right)x_6-{\mu }_{10}{x}_{10}.\end{array}$(46) ) represents the HHV-8 dynamics.

To this end, we consider the objective (or cost) functional (47) $\begin{array}{r}J\left(u_1,u_2,u_3\right)={\int }_0^T\left[A_1{x}_5+A_2{x}_6+A_3{x}_7+\frac{1}{2}{B}_1{u}_1^2+\frac{1}{2}{B}_2{u}_2^2+\frac{1}{2}{B}_3{u}_3^2\right]\mathrm{d}t,\end{array}$(47) where the control functions $u_1\left(t\right),u_2\left(t\right)$ and $u_3\left(t\right)$ are bounded, Lebesgue integrable functions on $[0,T_f]$ and $A_i,B_i,i=1,2,3$ are positive constants.

6.1. Existence of optimal control

Our control problem is formulated by minimizing the functional J subject to the system (37(37) $\begin{array}{rl}\dot{x_1}& ={\mathrm{\Pi }}_1+\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)x_9}{x_9+S_9}{\mathrm{\Pi }}_1+{\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x}}_{10}\\ & +S_{10}{\mathrm{\Pi }}_1-{\mu }_1{x}_1-\left(1-u_1\left(t\right)\right){\text{β}}_1{x}_1{x}_9.\end{array}$(37) )–(46(46) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6\left(1-u_2\left(t\right)\right)x_6-{\mu }_{10}{x}_{10}.\end{array}$(46) ). That is, we seek to find optimal controls $u_1^{\ast },u_2^{\ast }$ and $u_3^{\ast }$ such that (48) $\begin{array}{r}J\left(u_1^{\ast },u_2^{\ast },u_3^{\ast }\right)=min\left\{J\left(u_1,u_2,u_3\right)|u_1,u_2,u_3\in U\right\}\end{array}$(48) where (49) $\begin{array}{r}U=\left\{(u_1,u_2,u_3)\mathrm{s}\mathrm{u}\mathrm{c}\mathrm{h}\mathrm{t}\mathrm{h}\mathrm{a}\mathrm{t}{u}_1,u_2,u_3\mathrm{a}\mathrm{r}\mathrm{e}\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}0\le u_i\le 1,\mathrm{f}\mathrm{o}\mathrm{r}t\in [0,T_f]\right\}\end{array}$(49) is the control set.

The necessary conditions that an optimal solution must satisfy come from the Pontryagin et al. [26] maximum principle. This principle converts the system (37(37) $\begin{array}{rl}\dot{x_1}& ={\mathrm{\Pi }}_1+\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)x_9}{x_9+S_9}{\mathrm{\Pi }}_1+{\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x}}_{10}\\ & +S_{10}{\mathrm{\Pi }}_1-{\mu }_1{x}_1-\left(1-u_1\left(t\right)\right){\text{β}}_1{x}_1{x}_9.\end{array}$(37) )–(46(46) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6\left(1-u_2\left(t\right)\right)x_6-{\mu }_{10}{x}_{10}.\end{array}$(46) ) with Equation (47(47) $\begin{array}{r}J\left(u_1,u_2,u_3\right)={\int }_0^T\left[A_1{x}_5+A_2{x}_6+A_3{x}_7+\frac{1}{2}{B}_1{u}_1^2+\frac{1}{2}{B}_2{u}_2^2+\frac{1}{2}{B}_3{u}_3^2\right]\mathrm{d}t,\end{array}$(47) ) into a problem of minimizing pointwise a Hamiltonian H, with respect to $u_1,u_2$ and $u_3.$ The Hamiltonian function of the optimal problem is given by (50) $\begin{array}{rl}& H(\mathbf{x},\mathbf{u},{\lambda }_{\mathbf{x}},t)\\ & =A_1{x}_5+A_2{x}_6+A_3{x}_7+\frac{1}{2}{B}_1{u}_1^2+\frac{1}{2}{B}_2{u}_2^2+\frac{1}{2}{B}_3{u}_3^2\\ & +{\lambda }_{x_1}\left({\mathrm{\Pi }}_1+{\alpha }_5(1-u_3){\mathrm{\Pi }}_1\left(\frac{x_9}{x_9+S_9}+\frac{ϵx_{10}}{x_{10}+S_{10}}\right)-{\mu }_1{x}_1-(1-u_1){\text{β}}_1{x}_1{x}_9\right)\\ & +{\lambda }_{x_2}\left({\mathrm{\Pi }}_2+{\alpha }_5(1-u_3){\mathrm{\Pi }}_2\left(\frac{x_9}{x_9+S_9}+\frac{ϵx_{10}}{x_{10}+S_{10}}\right)-{\mu }_2{x}_2-(1-u_1){\text{β}}_2{x}_2{x}_{10}\right)\\ & +{\lambda }_{x_3}\left({\mathrm{\Pi }}_3+\frac{c_{10}{\mathrm{\Pi }}_3{x}_{10}}{x_{10}+f_{10}}-{\mu }_3{x}_3\right)+{\lambda }_{x_4}\left({\mathrm{\Pi }}_4+\frac{c_9{\mathrm{\Pi }}_4{x}_9}{x_9+f_9}-{\mu }_4{x}_4\right)\\ & +{\lambda }_{x_5}\left({\alpha }_5(1-u_3){\mathrm{\Pi }}_1\left(\frac{x_9}{x_9+S_9}+\frac{ϵx_{10}}{x_{10}+S_{10}}\right)+(1-u_1){\text{β}}_1{x}_1{x}_9-m_4{x}_4{x}_5-{\mu }_5{x}_5\right)\\ & +{\lambda }_{x_6}\left({\alpha }_5(1-u_3){\mathrm{\Pi }}_2\left(\frac{x_9}{x_9+S_9}+\frac{ϵx_{10}}{x_{10}+S_{10}}\right)+(1-u_1){\text{β}}_2{x}_2{x}_{10}-m_3{x}_3{x}_6-{\mu }_6{x}_6\right)\\ & +{\lambda }_{x_7}\left((1-u_1)r_7{x}_{10}\left(1-\frac{x_{10}}{x_{{10}_{max}}}\right)-{\mu }_7{x}_7-d_3{x}_3{x}_7\right)+{\lambda }_{x_8}\left({\mu }_7{x}_7-{\mu }_8{x}_8\right)\\ & +{\lambda }_{x_9}\left(N_9{\mu }_5(1-u_2)x_5-{\mu }_9{x}_9\right)+{\lambda }_{x_{10}}\left(N_{10}{\mu }_6(1-u_2)x_6-{\mu }_{10}{x}_{10}\right)\end{array}$(50) where $\mathbf{x}=(x_1,x_2,\dots ,x_{10}),\mathbf{u}=(u_1,u_2,u_3),{\lambda }_x=({\lambda }_{x_1},{\lambda }_{x_2},\dots ,{\lambda }_{x_{10}})$ and ${\lambda }_{x_i},i=1,2,\dots ,10$ are the adjoint or co-state variable corresponding to the state variable $x_i.$ The system of equations is found by taking appropriate partial derivatives of the Hamiltonian with respect to the associated state variable, $x_i.$

Theorem 6.1

Given optimal controls $u_1^{\ast },u_2^{\ast },u_3^{\ast }$ that minimize $J(u_1,u_2,u_3)$ over U, and the solutions $x_1,x_2,\dots x_{10}$ of the corresponding state system (37(37) $\begin{array}{rl}\dot{x_1}& ={\mathrm{\Pi }}_1+\frac{\left({\alpha }_5(1-u_3(t\left)\right)\right)x_9}{x_9+S_9}{\mathrm{\Pi }}_1+{\frac{\left(ϵ{\alpha }_5(1-u_3(t\left)\right)\right)x_{10}}{x}}_{10}\\ & +S_{10}{\mathrm{\Pi }}_1-{\mu }_1{x}_1-\left(1-u_1\left(t\right)\right){\text{β}}_1{x}_1{x}_9.\end{array}$(37) )–(46(46) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6\left(1-u_2\left(t\right)\right)x_6-{\mu }_{10}{x}_{10}.\end{array}$(46) ) with Equation (47(47) $\begin{array}{r}J\left(u_1,u_2,u_3\right)={\int }_0^T\left[A_1{x}_5+A_2{x}_6+A_3{x}_7+\frac{1}{2}{B}_1{u}_1^2+\frac{1}{2}{B}_2{u}_2^2+\frac{1}{2}{B}_3{u}_3^2\right]\mathrm{d}t,\end{array}$(47) ), then there exist adjoint variables ${\lambda }_{x_i}$ satisfying $\frac{\mathrm{\partial }H}{\mathrm{\partial }{x}_i}=-\frac{d{\lambda }_{x_i}}{dt},{\lambda }_{x_i}\left(T_f\right)=0,i=1,2,\dots ,10$and $u_i^{\ast }=min\left\{1,max(0,{\hat{u}}_i)\right\},i=1,2u_3^{\ast }=min\left\{0.5,max(0,{\hat{u}}_3)\right\}$where $u_i^{\ast },i=1,2,3$ are given by $\begin{array}{rl}{\hat{u}}_1& =\frac{x_1{\lambda }_{x_1}+{\text{β}}_1{x}_1{x}_9({\lambda }_{x_5}-{\lambda }_{x_1})+{\text{β}}_2{x}_2{x}_{10}({\lambda }_{x_6}-{\lambda }_{x_2})+r_7{x}_{10}\left(1-\frac{x_{10}}{x_{{10}_{max}}}\right){\lambda }_{x_7}}{B_1},\\ {\hat{u}}_2& =\frac{N_9{\mu }_5{x}_5{\lambda }_{x_9}+N_{10}{\mu }_6{x}_6{\lambda }_{x_{10}}}{B_2},\\ {\hat{u}}_3& =\left(\frac{{\alpha }_5{x}_9}{x_9+S_9}+\frac{ϵ{\alpha }_5{x}_{10}}{x_{10}+S_{10}}\right)\frac{({\mathrm{\Pi }}_1({\lambda }_{x_1}+{\lambda }_{x_5})+{\mathrm{\Pi }}_2({\lambda }_{x_2}+{\lambda }_{x_6}))}{B_3}.\end{array}$

Proof.

Corollary 4.1 of Fleming and Rishel [12] gives the existence of an optimal control due to the convexity of the integrand of J with respect to $u_1,u_2$ and $u_3$. In addition, it also guarantees a-priori boundedness of the state solutions and the Lipschitz property of the state system with respect to the state variables. The differential equations governing the adjoint variables are obtained by differentiation of the Hamiltonian function, evaluated at the optimal control. Then the adjoint differential equations can be written as $\begin{array}{rl}\frac{\mathrm{d}{\lambda }_{x_1}}{\mathrm{d}t}& ={\mu }_1{\lambda }_{x_1}-(1-u_1){\text{β}}_1{x}_9({\lambda }_{x_5}-{\lambda }_{x_1}),\\ \frac{\mathrm{d}{\lambda }_{x_2}}{\mathrm{d}t}& ={\mu }_2{\lambda }_{x_2}-(1-u_1){\text{β}}_2{x}_{10}({\lambda }_{x_6}-{\lambda }_{x_2}),\\ \frac{\mathrm{d}{\lambda }_{x_3}}{\mathrm{d}t}& ={\mu }_3{\lambda }_{x_3}+m_3{x}_6{\lambda }_{x_6}+d_3{x}_7{\lambda }_{x_7},\\ \frac{\mathrm{d}{\lambda }_{x_4}}{\mathrm{d}t}& =\mu {\lambda }_{x_4}+m_4{x}_5{\lambda }_{x_5},\\ \frac{\mathrm{d}{\lambda }_{x_5}}{\mathrm{d}t}& =-A_1+({\mu }_5+m_4{x}_4){\lambda }_{x_5}-N_9{\mu }_5(1-u_2){\lambda }_{x_9},\\ \frac{\mathrm{d}{\lambda }_{x_6}}{\mathrm{d}t}& =-A_2+({\mu }_6+m_3{x}_3){\lambda }_{x_6}-N_{10}{\mu }_6(1-u_2){\lambda }_{x_{10}},\\ \frac{\mathrm{d}{\lambda }_{x_7}}{\mathrm{d}t}& =-A_3+({\mu }_7+d_3{x}_3){\lambda }_{x_7}-{\mu }_7{\lambda }_{x_8},\\ \frac{\mathrm{d}{\lambda }_{x_8}}{\mathrm{d}t}& ={\mu }_8{\lambda }_{x_8},\\ \frac{\mathrm{d}{\lambda }_{x_9}}{\mathrm{d}t}& =-\frac{(1-u_3){\alpha }_5{S}_9}{(x_9+S_9{)}^2}({\mathrm{\Pi }}_1({\lambda }_{x_1}+{\lambda }_{x_5})+{\mathrm{\Pi }}_2({\lambda }_{x_2}+{\lambda }_{x_6}))\\ & +(1-u_1){\text{β}}_1{x}_1({\lambda }_{x_1}-{\lambda }_{x_5})-\frac{c_9{\mathrm{\Pi }}_4{f}_9}{(x_9+f_9{)}^2}{\lambda }_{x_4}+{\mu }_9{\lambda }_{x_9},\\ \frac{\mathrm{d}{\lambda }_{x_{10}}}{\mathrm{d}t}& =-\frac{(1-u_3)ϵ{\alpha }_5{S}_{10}}{(x_{10}+S_{10}{)}^2}({\mathrm{\Pi }}_1({\lambda }_{x_1}+{\lambda }_{x_5})+{\mathrm{\Pi }}_2({\lambda }_{x_2}+{\lambda }_{x_6}))\\ & +(1-u_1){\text{β}}_2{x}_2({\lambda }_{x_2}-{\lambda }_{x_6})-\frac{c_{10}{\mathrm{\Pi }}_3{f}_{10}}{(x_{10}+f_{10}{)}^2}{\lambda }_{x_3}\\ & +{\mu }_{10}{\lambda }_{x_{10}}-\frac{(1-u_1)r_7(x_{{10}_{max}-2x_{10}})}{x_{{10}_{max}}}\end{array}$In what follows, we differentiate the Hamiltonian H with respect to $u_i\left(t\right),i=1,2,3$ to obtain $u_i^{\ast }=min\{1,max(0,{\hat{u}}_1)\},i=1,2.u_3^{\ast }=min\{0.5,max(0,{\hat{u}}_3)\}.$By standard control arguments involving the bounds on the controls, we conclude $u_i^{\ast }=\left\{\begin{array}{ll}{\hat{u}}_i,& \mathrm{i}\mathrm{f}0<{\hat{u}}_i<1;\\ 0,& \mathrm{i}\mathrm{f}{\hat{u}}_i\le 0;\\ 1,& \mathrm{i}\mathrm{f}{\hat{u}}_i\ge 1.\end{array}\right.$which can be compactly written as $u_i^{\ast }=min\{1,max(0,{\hat{u}}_i)\},i=1,2\text{and}{u}_3^{\ast }=\left\{\begin{array}{ll}{\hat{u}}_3,& \mathrm{i}\mathrm{f}0<{\hat{u}}_3<0.5;\\ 0,& \mathrm{i}\mathrm{f}{\hat{u}}_3\le 0;\\ 1,& \mathrm{i}\mathrm{f}{\hat{u}}_3\ge \mathrm{0.5.}\end{array}\right.$Analogously, this can be written as $u_3^{\ast }=min\{0.5,max(0,{\hat{u}}_3)\}\text{■}$

7. Numerical simulations of optimal treatment regiments

Figure 9 shows the profiles of the optimal variables $u_1^{\ast },u_2^{\ast }$ and $u_3^{\ast }$ in the first 20 days. The initial values $u_i^{\ast }\left(20\right),i=1,2,3$ are used in an iterative procedure for periods much longer than 20 days. Figure 10 depicts the dynamics of the infected cell populations and viral loads for fixed values of $u_2$ and $u_3$ during 300 days for different values of the HAART efficacy parameter $u_1$. For optimal efficacy $u_1^{\ast }\approx 0.79$, the HIV-1 and HHV-8 populations drop below the level of detection. This parameter combination, however, is not unique as illustrated in Figure 11, which shows that reduced infected CD4 T-cell-specific HAART efficacy (from $u_1^{\ast }\approx 0.79$ to $u_1^{\ast }\approx 0.67$) and increased infected B-cell-specific HAART efficacy (from $u_2^{\ast }\approx 0.48$ to $u_2^{\ast }\approx 0.68$) in combination with reduced KS therapy-specific efficacy (from $u_3^{\ast }\approx 0.25$ to $u_3^{\ast }\approx 0.15$) is also able to control HIV-1 and HHV-8 co-infection. It is possible to find the triple $(u_1^{\ast },u_2^{\ast },u_3^{\ast })$ but this is not unique. This conclusion makes it easier to develop cheaper combinations that would suit the budgets of developing countries.

Figure 9. Optimal controls in the first 20 days. (a) Control profiles of optimal controls, $u_1^{\ast },u_2^{\ast },u_3^{\ast }$.

Figure 10. Population dynamics of infected cells and viruses during 300 days for varying values of the infected CD4 T-cell-specific HAART efficacy $u_1$: (a) infected T cells with $u_2=0.48,u_3=0.25$; (b) infected B cells with $u_2=0.48,u_3=0.25$; (c) infected progenitor cells with $u_2=0.48,u_3=0.25$; (d) KS cells with $u_2=0.48,u_3=0.25$; (e) HIV-1 with $u_2=0.48,u_3=0.25$; (f) HHV-8 with $u_2=0.48,u_3=0.25$.

Figure 11. Population dynamics of infected cells and viruses during 300 days for a new combination of optimal controls: (a) Infected T cells with $u_2=0.68,u_3=0.15$; (b) infected B cells with $u_2=0.68,u_3=0.15$; (c) infected progenitor cells with $u_2=0.68,u_3=0.15$; (d) KS cells with $u_2=0.68,u_3=0.15$; (e) HIV-1 with $u_2=0.68,u_3=0.15$; (f) HHV-8 with $u_2=0.68,u_3=0.15$.

8. Discussion

Currently there is no treatment available to eradicate HHV-8 infection and the purpose of anti-KS therapies is directed at slowing disease progression. KS is the most common neoplasm associated with AIDS and hence, therapies centre on the use of HAART in AIDS-KS patients. In this work, we formulate two models with innate and adaptive mechanism, in order to study the dynamics of non-AIDS KS and AIDS-KS and to make predictions about the efficacy of anti-KS therapies.

Some evidence suggests that immune activation is a requisite for the development of Classic KS [11,22]. Hence, we included the effect of the immune response in the MIM and showed that in silico treatment of NAKS can significantly reduce HHV-8 infection. We determined a critical value of the efficacy threshold for the innate immune response below which KS burden is diminished even if ${\mathcal{R}}_0>1$. This is completely novel to the best of our knowledge and has important implications, since therapies involving cytokines such as IL-2 [28] may have adverse side effects. Because it is not known how long the innate response for HIV-1 and HHV-8 co-infection lasts, this study is not able to suggest the dosage and frequency of treatment regimen. It is hoped a clinical study can explore the potential revealed by this study.

About 30% of Classic KS patients develop a second malignancy like non-Hodgkin lymphoma [10], hence, an early diagnosis that may prevent second malignancies is essential. The MIM demonstrated the potential for controlling HHV-8 infection and consequently, also HIV-1 and HHV-8 co-infection. Currently, there are early HIV-1 tests capable of revealing that an individual has been in contact with someone infected with HIV-1 (http://hivinsite.ucsf.edu/insite?page=basics-01-01). The strategy described in this study can be used to prevent HIV-1 infections before they develop for example in rape victims who are treated with antiretroviral drugs for 28 days immediately after the crime.

In the absence of any treatment the long-term prognosis of AIDS-KS is poor. This is shown by Figures 7 and 8 of the MAM. However, if cART (e.g. HAART plus KS therapy) is administered at optimal levels, both HIV-1 and HHV-8 infection can decline to undetectable levels (Figures 10 and 11). The simulations indicate that HAART at suboptimal levels cannot eradicate the viral load, which can possibly lead to the emergence of drug-resistance mutations [2]. For cART, when the efficacy of a potential fusion inhibitor or reverse transcriptase inhibitor is kept at optimal level, a very low level efficacy for KS therapy is sufficient to control the HHV-8 infection.

The implications of our study are three-fold. First, early intervention in the form of KS-therapy can control HHV-8 infection from developing into KS and possibly also slow HIV-1 infection to develop into AIDS. Second, optimal control of HIV-1 infection using HAART is an integral part of a successful AIDS-KS therapy and should be used in combination with low-level KS therapy. These recommendations have the potential to impact the therapeutic goal in KS, by focusing on short term control and thus, aiding long-term remission. Third, it is hoped that the existence of an infectivity threshold, which is critical for the progression of KS from asymptomatic to symptomatic stage, can help to study for example, the HHV-8 infectivity of transfusions [1].

The MIM revealed that externally administered cytokine drugs can reduce the chances of KS developing. This is supported in part by a study by [29,33] involving IL-2 and IL-12 which showed that administration of cytokines can prevent or delay the development of infection. However, to the best of our knowledge there are no clinical studies that have determined the safe levels of these cytokine drugs that can be taken. Because it is not known how long the innate response for HIV-1/HHV-8 co-infection lasts, this study is not able to suggest the dosage and frequency of treatment regimen. We believe a clinical study can explore the potential revealed by our study.

Abbreviations

KS, Kaposi's Sarcoma; AIDS, Acquired immune deficiency syndrome; HHV-8, Human Herpes virus-8; HAART, Highly Active Antiretroviral Therapy

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The first author gratefully acknowledges the funding received from Simons Foundation based at Botswana International University of Science and Technology (BIUST).

Appendix

(i) Positivity and boundedness of solution

Lemma

b

Let $y_1\left(t\right)=x_1\left(t\right)+x_5\left(t\right)$ and $y_2\left(t\right)=x_2\left(t\right)+x_6\left(t\right)$ denote the total sub-populations of the $CD{4}^{+}$ T cells and B cells at time t, respectively, where $x_i\left(t\right)\ge 0$ by Lemma $1\left(a\right)$.

Table A1. Range of parameters for sensitivity analysis for the MIM.

Parameter	Range	Units
$K_1$	0.11–0.39	day${}^{-1}$
$K_2$	0.0154–0.0193	day${}^{-1}$
$K_3$	0.013–0.054	day${}^{-1}$
${\theta }_{x_1}$	100–500	cell mm${}^{-3}$
${\theta }_{x_5}$	0.002–0.069
$x_{1_{max}}$	200–650	cell mm${}^{-3}$
$x_{2_{max}}$	${10}^4-5.1\times {10}^5$	virions mm${}^{-3}$
${\mu }_{x_3}$	0.1–0.3	day${}^{-1}$
${\mu }_{x_4}$	0.001–0.067	ml day${}^{-1}$
${\mu }_{x_1}$	0.2–0.5	day${}^{-1}$
${\mu }_{x_2}$	0.45–0.65	day${}^{-1}$
$N_{x_2}$	100–1000	virions cell${}^{-1}$

Table A2. Range of parameters for sensitivity analysis for the MAM.

Parameter	Range	Units
$N_9$	100–1000	virions cell${}^{-1}$
$N_{10}$	100–1000	virions cell${}^{-1}$
${\mathrm{\Pi }}_1$	${10}^3-{10}^4$	cell ml${}^{-1}$ day${}^{-1}$
${\mathrm{\Pi }}_2$	$70\times {10}^3-4.275\times {10}^5$	cell ml${}^{-1}$ day${}^{-1}$
${\mathrm{\Pi }}_3$	${10}^2-{10}^3$	cell ml${}^{-1}$ day${}^{-1}$
${\mathrm{\Pi }}_4$	${10}^4-3\times {10}^4$	cell ml${}^{-1}$ day${}^{-1}$
$r_7$	$0.045-0.055$	cells virion${}^{-1}$ day${}^{-1}$
$d_3$	$5.47\times {10}^{-5}-1.09\times {10}^{-4}$	ml cell${}^{-1}$ day${}^{-1}$
$m_3$	$5.47\times {10}^{-5}-1.09\times {10}^{-4}$	ml cell${}^{-1}$ day${}^{-1}$
$m_4$	${10}^{-7}-{10}^{-5}$	ml cell${}^{-1}$ day${}^{-1}$
${\alpha }_5$	${10}^{-2}-1.5\times {10}^{-4}$
${\alpha }_6$	$1.4\times {10}^{-2}-5.5\times {10}^{-2}$
$S_9$	$2\times {10}^5-5\times {10}^5$	virions ml${}^{-1}$
$S_{10}$	${10}^5-3\times {10}^5$	virions ml${}^{-1}$
${\text{β}}_1$	$2.4\times {10}^{-8}-2.5\times {10}^{-7}$	ml virion${}^{-1}$ day${}^{-1}$
${\text{β}}_2$	$4.75\times {10}^{-9}-4.75\times {10}^{-7}$	ml virion${}^{-1}$ day${}^{-1}$
${\mu }_9$	2–5	day${}^{-1}$
${\mu }_{10}$	0.45–0.65	day${}^{-1}$

Adding Equations (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) ) and (21(21) $\begin{array}{r}{\dot{x}}_5=\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}+{\text{β}}_1{x}_1{x}_9-m_4{x}_4{x}_5-{\mu }_5{x}_5.\end{array}$(21) ), and then (18(18) $\begin{array}{r}\dot{x_2}={\mathrm{\Pi }}_2+\frac{{\mathrm{\Pi }}_2{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_2{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_2{x}_2-{\text{β}}_2{x}_2{x}_{10}.\end{array}$(18) ) and (22(22) $\begin{array}{r}{\dot{x}}_6=\frac{{\mathrm{\Pi }}_2{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_2{\alpha }_6{x}_{10}}{x_{10}+S_{10}}+{\text{β}}_2{x}_2{x}_{10}-m_3{x}_3{x}_6-{\mu }_6{x}_6.\end{array}$(22) ) we get $\begin{array}{rl}\frac{\mathrm{d}{y}_j\left(t\right)}{\mathrm{d}t}& \le {\mathrm{\Pi }}_j\left(1+2{\alpha }_5\frac{x_9}{x_9+S_9}+2{\alpha }_6\frac{x_{10}}{x_{10}+S_{10}}\right)-{\mu }_j{x}_j-{\mu }_{j+4}{x}_{j+4},j=1,2.\\ & \le \underset{M_j}{\underbrace{{\mathrm{\Pi }}_j(1+2({\alpha }_5+{\alpha }_6))}}-{\hat{\mu }}_j{y}_j\left(t\right),\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\hat{\mu }}_j=min\{{\mu }_j,{\mu }_{j+4}\}\end{array}$It is easy to show that as $t\to \mathrm{\infty },y_j\left(t\right)\le \frac{M_j}{{\hat{\mu }}_j},j=1,2.$ Equations (25(25) $\begin{array}{r}{\dot{x}}_9=N_9{\mu }_5{x}_5-{\mu }_9{x}_9.\end{array}$(25) ) and (26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ) can be represented by $\frac{\mathrm{d}{x}_j\left(t\right)}{\mathrm{d}t}=N_j{\mu }_{j-4}{x}_{j-4}-{\mu }_j{x}_j\le N_j{\mu }_{j-4}\frac{M_{j-8}}{{\hat{\mu }}_{j-8}}-{\mu }_j{x}_{j}j=9,10.$It is easy to show that $\begin{array}{rl}\underset{t\to \mathrm{\infty }}{lim sup}{x}_j\left(t\right)& \le \frac{N_j{\mu }_{j-4}{M}_{j-8}}{{\mu }_j{\hat{\mu }}_{j-8}},j=9,10,\\ \underset{t\to \mathrm{\infty }}{lim sup}Z\left(t\right)& \le \frac{{\mathrm{\Pi }}_3(1+c_{10})+{\mathrm{\Pi }}_4(1+c_9)}{\xi },\\ \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\xi & =min\{{\mu }_3,{\mu }_4\}\mathrm{a}\mathrm{n}\mathrm{d}Z\left(t\right)=x_3\left(t\right)+x_4\left(t\right),\}\\ \underset{t\to \mathrm{\infty }}{lim sup}{x}_7\left(t\right)& \le \frac{{\mu }_6{r}_7{N}_{10}{M}_2}{{\hat{\mu }}_2{\mu }_7{\mu }_{10}}\left(1-\frac{{\mu }_6{N}_{10}{M}_2}{{\hat{\mu }}_2{\mu }_{10}{x}_{{10}_{max}}}\right)\\ \underset{t\to \mathrm{\infty }}{lim sup}{x}_8\left(t\right)& \le \frac{{\mu }_6{\mu }_7{r}_7{N}_{10}{M}_2}{{\hat{\mu }}_2{\mu }_8{\mu }_9{\mu }_{10}}\left(1-\frac{{\mu }_6{N}_{10}{M}_2}{{\hat{\mu }}_2{\mu }_{10}{x}_{{10}_{max}}}\right)\end{array}$

(ii) Virus Free Equilibrium and the Basic Reproduction Number

The following refers to the reproduction numbers given in Section 4.3 (A1) $\begin{array}{rl}{\mathcal{R}}_{V_1}& =\frac{N_9{\mathrm{\Pi }}_1\left({\alpha }_5{\mu }_1+{\text{β}}_1{S}_9\right)}{{\mu }_1{\mu }_9{S}_9}.\frac{{\mu }_4{\mu }_5}{({\mu }_4{\mu }_5+m_4{\mathrm{\Pi }}_4)}<\frac{N_9{\mathrm{\Pi }}_1\left({\alpha }_5{\mu }_1+{\text{β}}_1{S}_9\right)}{{\mu }_1{\mu }_9{S}_9}={\mathcal{R}}_1.\end{array}$(A1) (A2) $\begin{array}{rl}{\mathcal{R}}_{V_8}& =\frac{N_{10}{\mathrm{\Pi }}_2\left({\alpha }_6{\mu }_2+{\text{β}}_2{S}_{10}\right)}{{\mu }_2{\mu }_{10}{S}_{10}}.\frac{{\mu }_3{\mu }_6}{({\mu }_3{\mu }_6+m_3{\mathrm{\Pi }}_3)}<\frac{N_{10}{\mathrm{\Pi }}_2({\alpha }_6{\mu }_2+{\text{β}}_2{S}_{10})}{{\mu }_2{\mu }_{10}{S}_{10}}={\mathcal{R}}_8.\end{array}$(A2) From (A2(A2) $\begin{array}{rl}{\mathcal{R}}_{V_8}& =\frac{N_{10}{\mathrm{\Pi }}_2\left({\alpha }_6{\mu }_2+{\text{β}}_2{S}_{10}\right)}{{\mu }_2{\mu }_{10}{S}_{10}}.\frac{{\mu }_3{\mu }_6}{({\mu }_3{\mu }_6+m_3{\mathrm{\Pi }}_3)}<\frac{N_{10}{\mathrm{\Pi }}_2({\alpha }_6{\mu }_2+{\text{β}}_2{S}_{10})}{{\mu }_2{\mu }_{10}{S}_{10}}={\mathcal{R}}_8.\end{array}$(A2) ) we can deduce that

• ${\mathcal{R}}_{V_i},i=1,8$ are the reproduction numbers when virus i specific effector cells are active but their effect varies from being weak to perfect as $m_j,j=3,4,$ increases · 

• ${\mathcal{R}}_i,i=1,8$ are the reproduction numbers when virus i specific effector cells are dysfunctional.

1. Let $x=(x_5,x_6,x_7,x_8,x_9,x_{10}{)}^T,y=(x_1,x_2,x_3,x_4{)}^T$. The system (17(17) $\begin{array}{r}\dot{x_1}={\mathrm{\Pi }}_1+\frac{{\mathrm{\Pi }}_1{\alpha }_5{x}_9}{x_9+S_9}+\frac{{\mathrm{\Pi }}_1{\alpha }_6{x}_{10}}{x_{10}+S_{10}}-{\mu }_1{x}_1-{\text{β}}_1{x}_1{x}_9.\end{array}$(17) )– 26(26) $\begin{array}{r}{\dot{x}}_{10}=N_{10}{\mu }_6{x}_6-{\mu }_{10}{x}_{10}.\end{array}$(26) ) can be written as $\begin{array}{rl}\frac{\mathrm{d}x}{\mathrm{d}t}& =\mathcal{F}(x,y)-\mathcal{V}(x,y),\\ \frac{\mathrm{d}y}{\mathrm{d}t}& =g_j(x,y)\end{array}$The matrices for new infections and other class transition, respectively given by F and V, are (A3) $\begin{array}{r}\begin{array}{rl}F& =\left(\begin{array}{cccccc}0& 0& 0& 0& \frac{{\mathrm{\Pi }}_1{K}_1}{{\mu }_1{S}_9}& \frac{{\alpha }_6{\mathrm{\Pi }}_1}{S_{10}}\\ 0& 0& 0& 0& \frac{{\alpha }_5{\mathrm{\Pi }}_2}{S_9}& \frac{{\mathrm{\Pi }}_2{K}_2}{{\mu }_2{S}_{10}}\\ 0& 0& 0& 0& 0& r_7\\ 0& 0& {\mu }_7& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right),\\ V& =\left(\begin{array}{cccccc}\frac{a_{11}}{{\mu }_4}& 0& 0& 0& 0& 0\\ 0& \frac{a_{22}}{{\mu }_3}& 0& 0& 0& 0\\ 0& 0& \frac{a_{33}}{{\mu }_3}& 0& 0& 0\\ 0& 0& 0& {\mu }_8& 0& 0\\ -N_9{\mu }_5& 0& 0& 0& {\mu }_9& 0\\ 0& -N_{10}{\mu }_6& 0& 0& 0& {\mu }_{10}\end{array}\right)\\ F{V}^{-1}& =\left(\begin{array}{cccccc}\frac{{\mu }_4{\mu }_5{N}_9{\mathrm{\Pi }}_1{K}_1}{{\mu }_1{\mu }_9{a}_{11}{S}_9}& \frac{{\mu }_3{\mu }_6{\alpha }_6{N}_{10}{\mathrm{\Pi }}_1}{{\mu }_{10}{a}_{22}{S}_{10}}& 0& 0& \frac{{\mathrm{\Pi }}_1{K}_1}{{\mu }_1{\mu }_9{S}_9}& \frac{{\alpha }_6{\mathrm{\Pi }}_1}{{\mu }_{10}{S}_{10}}\\ \frac{{\mu }_4{\mu }_5{\alpha }_5{N}_9{\mathrm{\Pi }}_2}{{\mu }_9{a}_{11}{S}_9}& \frac{{\mu }_3{\mu }_6{N}_{10}{\mathrm{\Pi }}_2{K}_2}{{\mu }_2{\mu }_{10}{a}_{22}{S}_{10}}& 0& 0& \frac{{\alpha }_5{\mathrm{\Pi }}_2}{{\mu }_9{S}_9}& \frac{{\mathrm{\Pi }}_2{K}_2}{{\mu }_2{\mu }_{10}{S}_{10}}\\ 0& \frac{{\mu }_3{\mu }_6{r}_7{N}_{10}}{{\mu }_{10}{a}_{22}}& 0& 0& 0& \frac{r_7}{{\mu }_{10}}\\ 0& 0& \frac{{\mu }_3{\mu }_7}{a_{33}}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)\\ K_1& ={\alpha }_5{\mu }_1+{\text{β}}_1{S}_9,K_2={\alpha }_6{\mu }_2+{\text{β}}_2{S}_{10},a_{11}={\mu }_4{\mu }_5+m_4{\mathrm{\Pi }}_4,\\ a_{22}& ={\mu }_3{\mu }_6+m_3{\mathrm{\Pi }}_3,a_{33}={\mu }_3{\mu }_7+d_3{\mathrm{\Pi }}_3\end{array}\end{array}$(A3) (iii) Global stability for the virus free equilibrium, ${ϵ}^0$

Lemma A.1

The lemma is based on the work of Castillo-Chavez et al. [9]. Consider the system (A4) $\begin{array}{r}\begin{array}{rl}\frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}& =F(\mathbf{(}\mathbf{X}\mathbf{,}\mathbf{I}\mathbf{)}),\\ \frac{\mathrm{d}\mathbf{I}}{\mathrm{d}t}& =G(\mathbf{(}\mathbf{X}\mathbf{,}\mathbf{I}\mathbf{)}),G(\mathbf{X}\mathbf{,}\mathbf{0})=0,\end{array}\end{array}$(A4) where $\mathbf{X}\in {\mathbb{R}}^m$ denotes the components of the uninfected states, $\mathbf{I}\in {\mathbb{R}}^n$ denotes the components of the infected states and ${\mathbf{U}}_{\mathbf{0}}=({\mathbf{X}}^{\ast },\mathbf{0})$ denotes the disease-free equilibrium of (A4(A4) $\begin{array}{r}\begin{array}{rl}\frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}& =F(\mathbf{(}\mathbf{X}\mathbf{,}\mathbf{I}\mathbf{)}),\\ \frac{\mathrm{d}\mathbf{I}}{\mathrm{d}t}& =G(\mathbf{(}\mathbf{X}\mathbf{,}\mathbf{I}\mathbf{)}),G(\mathbf{X}\mathbf{,}\mathbf{0})=0,\end{array}\end{array}$(A4) ). Assume the conditions (H1) and (H2) below are satisfied $\begin{array}{rl}\left(\mathrm{H}1\right)\mathrm{F}\mathrm{o}\mathrm{r}\frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}& =F(\mathbf{X}\mathbf{,}\mathbf{0}),{\mathbf{X}}^{\mathbf{\ast }}\mathrm{i}\mathrm{s}\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{y}\mathrm{a}\mathrm{s}\mathrm{y}\mathrm{m}\mathrm{p}\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{y}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}(\mathrm{g}.\mathrm{a}.\mathrm{s}),\\ \left(\mathrm{H}2\right)G(\mathbf{X}\mathbf{,}\mathbf{I})& =\mathbf{A}\mathbf{I}-\hat{G}(\mathbf{X}\mathbf{,}\mathbf{I}),\hat{G}(\mathbf{X}\mathbf{,}\mathbf{I})\ge \mathbf{0}\mathrm{f}\mathrm{o}\mathrm{r}(\mathbf{X}\mathbf{,}\mathbf{I})\in \mathrm{\Omega },\end{array}$where $\mathbf{A}=D_{\mathbf{I}}G\left({\mathbf{X}}^{\mathbf{\ast }}\mathbf{,}\mathbf{0}\right)$ is an M-matrix (the off diagonal elements of A are nonnegative) and Ω is the region where the model makes biological sense. Then ${\mathbf{U}}_{\mathbf{0}}$ is globally asymptotically stable.

(iv) AIDS-KS present equilibrium, $E^{\ast }$

The coordinates $x_9^{\ast }$ and $x_{10}^{\ast }$ in Theorem 4.2 are positive solutions to $Q_i,i=1,2$ where (A5) $\begin{array}{r}{Q}_1\left(x_{10}\right)=B_4{x}_{10}^4+B_3{x}_{10}^3+B_2{x}_{10}^2+B_1{x}_{10}+B_0=0,\end{array}$(A5) with (A6) $\begin{array}{r}\left\{\begin{array}{l}{B}_4=-{\text{β}}_2{\mu }_{10}(x_9+S_9)({\mu }_3{\mu }_6+(1+c_{10})m_3{\mathrm{\Pi }}_3)<0,\\ B_3=N_{10}{\mathrm{\Pi }}_2{S}_9{\text{β}}_2{\mu }_3{\mu }_6-S_9{\mu }_2{\mu }_3{\mu }_6{\mu }_{10}-{\mathrm{\Pi }}_3{m}_3{\mu }_2{\mu }_{10}{x}_9-{\mu }_2{\mu }_3{\mu }_6{\mu }_{10}{x}_9-{\mathrm{\Pi }}_3{S}_9{m}_3{\mu }_2{\mu }_{10}\\ -{\mathrm{\Pi }}_3{S}_9{S}_{10}{\text{β}}_2{m}_3{\mu }_{10}-{\mathrm{\Pi }}_3{S}_9{\text{β}}_2{f}_{10}{m}_3{\mu }_{10}-{\mathrm{\Pi }}_3{S}_9{c}_{10}{m}_3{\mu }_2{\mu }_{10}\\ -S_9{S}_{10}{\text{β}}_2{\mu }_3{\mu }_6{\mu }_{10}+N_{10}{\mathrm{\Pi }}_2{\text{β}}_2{\mu }_3{\mu }_6{x}_9\\ -{\mathrm{\Pi }}_3{S}_{10}{\text{β}}_2{m}_3{\mu }_{10}{x}_9-S_9{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{\mu }_{10}\\ -{\mathrm{\Pi }}_3{\text{β}}_2{f}_{10}{m}_3{\mu }_{10}{x}_9-{\mathrm{\Pi }}_3{c}_{10}{m}_3{\mu }_2{\mu }_{10}{x}_9-S_{10}{\text{β}}_2{\mu }_3{\mu }_6{\mu }_{10}{x}_9\\ -{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{\mu }_{10}{x}_9+N_{10}{\mathrm{\Pi }}_1{\alpha }_5{\text{β}}_2{\mu }_3{\mu }_6{x}_9+N_{10}{\mathrm{\Pi }}_1{\alpha }_6{\text{β}}_2{\mu }_3{\mu }_6{x}_9+\\ N_{10}{\mathrm{\Pi }}_2{\alpha }_5{\text{β}}_2{\mu }_3{\mu }_6{x}_9+N_{10}{\mathrm{\Pi }}_2{\alpha }_6{\text{β}}_2{\mu }_3{\mu }_6{x}_9\\ -{\mathrm{\Pi }}_3{S}_{10}{\text{β}}_2{c}_{10}{m}_3{\mu }_{10}{x}_9+N_{10}{\mathrm{\Pi }}_1{S}_9{\alpha }_6{\text{β}}_2{\mu }_3{\mu }_6+N_{10}{\mathrm{\Pi }}_2{S}_9{\alpha }_6{\text{β}}_2{\mu }_3{\mu }_6-{\mathrm{\Pi }}_3{S}_9{S}_{10}{\text{β}}_2{c}_{10}{m}_3{\mu }_{10}\\ B_2=N_{10}{\mathrm{\Pi }}_2{S}_{10}{\text{β}}_2{\mu }_3{\mu }_6{x}_9-{\mathrm{\Pi }}_3{S}_9{f}_{10}{m}_3{\mu }_2{\mu }_{10}-S_9{S}_{10}{\mu }_2{\mu }_3{\mu }_6{\mu }_{10}\\ -{\mathrm{\Pi }}_3{S}_{10}{m}_3{\mu }_2{\mu }_{10}{x}_9-S_9{f}_{10}{\mu }_2{\mu }_3{\mu }_6{\mu }_{10}\\ -{\mathrm{\Pi }}_3{f}_{10}{m}_3{\mu }_2{\mu }_{10}{x}_9-S_{10}{\mu }_2{\mu }_3{\mu }_6{\mu }_{10}{x}_9-f_{10}{\mu }_2{\mu }_3{\mu }_6{\mu }_{10}{x}_9\\ -{\mathrm{\Pi }}_3{S}_9{S}_{10}{m}_3{\mu }_2{\mu }_{10}-{\mathrm{\Pi }}_3{S}_9{S}_{10}{m}_3{\mu }_2{\mu }_{10}\\ -S_9{S}_{10}{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{\mu }_{10}+N_{10}{\mathrm{\Pi }}_2{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{x}_9-{\mathrm{\Pi }}_3{S}_{10}{\text{β}}_2{f}_{10}{m}_3{\mu }_{10}{x}_9+N_{10}{\mathrm{\Pi }}_1{\alpha }_5{\mu }_2{\mu }_3{\mu }_6{x}_9\\ +N_{10}{\mathrm{\Pi }}_1{\alpha }_6{\mu }_2{\mu }_3{\mu }_6{x}_9\\ -{\mathrm{\Pi }}_3{S}_{10}{c}_{10}{m}_3{\mu }_2{\mu }_{10}{x}_9-S_{10}{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{\mu }_{10}{x}_9+N_{10}{\mathrm{\Pi }}_2{S}_9{S}_{10}{\text{β}}_2{\mu }_3{\mu }_6\\ +N_{10}{\mathrm{\Pi }}_2{S}_9{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6-{\mathrm{\Pi }}_3{S}_9{S}_{10}{\text{β}}_2{f}_{10}{m}_3{\mu }_{10}\\ +N_{10}{\mathrm{\Pi }}_1{S}_9{\alpha }_6{\mu }_2{\mu }_3{\mu }_6+N_{10}{\mathrm{\Pi }}_1{S}_9{\alpha }_6{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6+N_{10}{\mathrm{\Pi }}_2{S}_9{\alpha }_6{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6\\ +N_{10}{\mathrm{\Pi }}_1{S}_{10}{\alpha }_5{\text{β}}_2{\mu }_3{\mu }_6{x}_9+N_{10}{\mathrm{\Pi }}_2{S}_{10}{\alpha }_5{\text{β}}_2{\mu }_3{\mu }_6{x}_9\\ +N_{10}{\mathrm{\Pi }}_1{\alpha }_5{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{x}_9+N_{10}{\mathrm{\Pi }}_1{\alpha }_6{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{x}_9+N_{10}{\mathrm{\Pi }}_2{\alpha }_5{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{x}_9\\ +N_{10}{\mathrm{\Pi }}_2{\alpha }_6{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{x}_9\\ B_1=N_{10}{\mathrm{\Pi }}_2{S}_9{S}_{10}{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6-S_9{S}_{10}{f}_{10}{\mu }_2{\mu }_3{\mu }_6{\mu }_{10}-{\mathrm{\Pi }}_3{S}_{10}{f}_{10}{m}_3{\mu }_2{\mu }_{10}{x}_9\\ -S_{10}{f}_{10}{\mu }_2{\mu }_3{\mu }_6{\mu }_{10}{x}_9-{\mathrm{\Pi }}_3{S}_9{S}_{10}{f}_{10}{m}_3{\mu }_2{\mu }_{10}\\ +N_{10}{\mathrm{\Pi }}_1{S}_9{\alpha }_6{f}_{10}{\mu }_2{\mu }_3{\mu }_6+N_{10}{\mathrm{\Pi }}_2{S}_{10}{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{x}_9+N_{10}{\mathrm{\Pi }}_1{S}_{10}{\alpha }_5{\mu }_2{\mu }_3{\mu }_6{x}_9\\ +N_{10}{\mathrm{\Pi }}_1{\alpha }_5{f}_{10}{\mu }_2{\mu }_3{\mu }_6{x}_9+N_{10}{\mathrm{\Pi }}_1{\alpha }_6{f}_{10}{\mu }_2{\mu }_3{\mu }_6{x}_9\\ +N_{10}{\mathrm{\Pi }}_1{S}_{10}{\alpha }_5{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{x}_9+N_{10}{\mathrm{\Pi }}_2{S}_{10}{\alpha }_5{\text{β}}_2{f}_{10}{\mu }_3{\mu }_6{x}_9\\ B_0=N_{10}{\mathrm{\Pi }}_1{S}_{10}{\alpha }_5{f}_{10}{\mu }_2{\mu }_3{\mu }_6{x}_9>0,\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}{x}_9>0.\end{array}\right.\end{array}$(A6) The existence of a positive value, $x_{10}^{\ast }$, for $Q_1\left(x_{10}\right)$ can be justified as follows:

Note that $Q_1\left(0\right)=B_0>0.$ Due to the continuity of $Q_1\left(x_{10}\right)$, we have $\underset{x_{10}\to \mathrm{\infty }}{lim}{Q}_1\left(x_{10}\right)=-\mathrm{\infty }$ since $B_4<0.$ Hence, there exists a positive value for $x_{10},$ say, $x_{10}^{\ast }$ such that $Q_1\left(x_{10}^{\ast }\right)=0.$

Using Equation (25(25) $\begin{array}{r}{\dot{x}}_9=N_9{\mu }_5{x}_5-{\mu }_9{x}_9.\end{array}$(25) ) along with $x_{10}^{\ast }$ just obtained above, we get (A7) $\begin{array}{r}{Q}_2\left(x_9\right)=E_4{x}_9^4+E_3{x}_9^3+E_2{x}_9^2+E_1{x}_9+E_0=0,\end{array}$(A7) with (A8) $\begin{array}{r}\left\{\begin{array}{l}{E}_4=-{\text{β}}_1{\mu }_9(x_{10}^{\ast }+S_{10})({\mu }_4{\mu }_5+({\mu }_9+c_9)m_4{\mathrm{\Pi }}_4)<0,\\ E_3=N_9{\mathrm{\Pi }}_1{S}_{10}{\text{β}}_1{\mu }_4{\mu }_5-S_{10}{\mu }_1{\mu }_4{\mu }_5{\mu }_9-{\mathrm{\Pi }}_4{m}_4{\mu }_1{\mu }_9{x}_{10}^{\ast }-{\mu }_1{\mu }_4{\mu }_5{\mu }_9{x}_{10}^{\ast }-{\mathrm{\Pi }}_4{S}_{10}{m}_4{\mu }_1{\mu }_9\\ -{\mathrm{\Pi }}_4{S}_9{S}_{10}{\text{β}}_1{m}_4{\mu }_9-{\mathrm{\Pi }}_4{S}_{10}{\text{β}}_1{f}_9{\mu }_4{\mu }_9-{\mathrm{\Pi }}_4{S}_{10}{c}_9{m}_4{\mu }_1{\mu }_9\\ -S_9{S}_{10}{\text{β}}_1{\mu }_4{\mu }_5{\mu }_9+N_9{\mathrm{\Pi }}_1{\text{β}}_1{\mu }_4{\mu }_5{x}_{10}^{\ast }\\ -{\mathrm{\Pi }}_4{S}_9{\text{β}}_1{m}_4{\mu }_9{x}_{10}^{\ast }-S_{10}{\text{β}}_1{f}_9{\mu }_4{\mu }_5{\mu }_9\\ -{\mathrm{\Pi }}_4{\text{β}}_1{f}_9{m}_4{\mu }_9{x}_{10}^{\ast }-{\mathrm{\Pi }}_4{c}_9{m}_4{\mu }_1{\mu }_9{x}_{10}^{\ast }-S_9{\text{β}}_1{\mu }_4{\mu }_5{\mu }_9{x}_{10}^{\ast }\\ +2N_9{\mathrm{\Pi }}_1{\alpha }_5{\text{β}}_1{\mu }_4{\mu }_5{x}_{10}^{\ast }+2N_9{\mathrm{\Pi }}_1{\alpha }_6{\text{β}}_1{\mu }_4{\mu }_5{x}_{10}^{\ast }-{\mathrm{\Pi }}_4{S}_9{\text{β}}_1{c}_9{m}_4{\mu }_9{x}_{10}^{\ast }\\ +2N_9{\mathrm{\Pi }}_1{S}_{10}{\alpha }_5{\text{β}}_1{\mu }_4{\mu }_5-{\mathrm{\Pi }}_4{S}_9{S}_{10}{\text{β}}_1{c}_9{m}_4{\mu }_9\\ E_2=N_9{\mathrm{\Pi }}_1{S}_9{\text{β}}_1{\mu }_4{\mu }_5{x}_{10}^{\ast }-{\mathrm{\Pi }}_4{S}_{10}{f}_9{m}_4{\mu }_1{\mu }_9-S_9{S}_{10}{\mu }_1{\mu }_4{\mu }_5{\mu }_9\\ -{\mathrm{\Pi }}_4{S}_9{m}_4{\mu }_1{\mu }_9{x}_{10}^{\ast }-S_{10}{f}_9{\mu }_1{\mu }_4{\mu }_5{\mu }_9\\ -{\mathrm{\Pi }}_4{f}_9{m}_4{\mu }_1{\mu }_9{x}_{10}^{\ast }-S_9{\mu }_1{\mu }_4{\mu }_5{\mu }_9{x}_{10}^{\ast }-f_9{\mu }_1{\mu }_4{\mu }_5{\mu }_9{x}_{10}^{\ast }\\ -{\mathrm{\Pi }}_4{S}_9{S}_{10}{c}_9{m}_4{\mu }_1{\mu }_9-{\mathrm{\Pi }}_4{S}_9{S}_{10}{m}_4{\mu }_1{\mu }_9\\ -S_9{S}_{10}{\text{β}}_1{f}_9{\mu }_4{\mu }_5{\mu }_9+N_9{\mathrm{\Pi }}_1{\text{β}}_1{f}_9{\mu }_4{\mu }_5{x}_{10}^{\ast }-{\mathrm{\Pi }}_4{S}_9{\text{β}}_1{f}_9{m}_4{\mu }_9{x}_{10}^{\ast }\\ +N_9{\mathrm{\Pi }}_1{\alpha }_5{\mu }_1{\mu }_4{\mu }_5{x}_{10}^{\ast }+N_9{\mathrm{\Pi }}_1{\alpha }_6{\mu }_1{\mu }_4{\mu }_5{x}_{10}^{\ast }\\ -{\mathrm{\Pi }}_4{S}_9{c}_9{m}_4{\mu }_1{\mu }_9{x}_{10}^{\ast }-S_9{\text{β}}_1{f}_9{\mu }_4{\mu }_5{\mu }_9{x}_{10}^{\ast }+N_9{\mathrm{\Pi }}_1{S}_9{S}_{10}{\text{β}}_1{\mu }_4{\mu }_5\\ +N_9{\mathrm{\Pi }}_1{S}_{10}{\text{β}}_1{f}_9{\mu }_4{\mu }_5-{\mathrm{\Pi }}_4{S}_9{S}_{10}{\text{β}}_1{f}_9{m}_4{\mu }_9\\ +N_9{\mathrm{\Pi }}_1{S}_{10}{\alpha }_5{\mu }_1{\mu }_4{\mu }_5+2N_9{\mathrm{\Pi }}_1{S}_{10}{\alpha }_5{\text{β}}_1{f}_9{\mu }_4{\mu }_5+2N_9{\mathrm{\Pi }}_1{S}_9{\alpha }_6{\text{β}}_1{\mu }_4{\mu }_5{x}_{10}^{\ast }\\ +2N_9{\mathrm{\Pi }}_1{\alpha }_5{\text{β}}_1{f}_9{\mu }_4{\mu }_5{x}_{10}^{\ast }+2N_9{\mathrm{\Pi }}_1{\alpha }_6{\text{β}}_1{f}_9{\mu }_4{\mu }_5{x}_{10}^{\ast }\\ E_1=N_9{\mathrm{\Pi }}_1{S}_9{S}_{10}{\text{β}}_1{f}_9{\mu }_4{\mu }_5-S_9{S}_{10}{f}_9{\mu }_{1}m{u}_4{\mu }_5{\mu }_9-{\mathrm{\Pi }}_4{S}_9{f}_9{m}_4{\mu }_1{\mu }_9{x}_{10}^{\ast }\\ -S_9{f}_9{\mu }_1{\mu }_4{\mu }_5{\mu }_9{x}_{10}^{\ast }-{\mathrm{\Pi }}_4{S}_9{S}_{10}{f}_9{m}_4{\mu }_1{\mu }_9\\ +N_9{\mathrm{\Pi }}_1{S}_{10}{\alpha }_5{f}_9{\mu }_1{\mu }_4{\mu }_5+N_9{\mathrm{\Pi }}_1{S}_9{\text{β}}_1{f}_9{\mu }_4{\mu }_5{x}_{10}^{\ast }+N_9{\mathrm{\Pi }}_1{S}_9{\alpha }_6{\mu }_1{\mu }_4{\mu }_5{x}_{10}^{\ast }\\ +N_9{\mathrm{\Pi }}_1{\alpha }_5{f}_9{\mu }_1{\mu }_4{\mu }_5{x}_{10}^{\ast }+N_9{\mathrm{\Pi }}_1{\alpha }_6{f}_9{\mu }_1{\mu }_4{\mu }_5{x}_{10}^{\ast }+2N_9{\mathrm{\Pi }}_1{S}_9{\alpha }_6{\text{β}}_1{f}_9{\mu }_4{\mu }_5{x}_{10}^{\ast }\\ E_0=N_9{\mathrm{\Pi }}_1{S}_9{\alpha }_6{f}_9{\mu }_1{\mu }_4{\mu }_5{x}_{10}^{\ast }>0.\end{array}\right.\end{array}$(A8) Analogously, arguing as before, we can establish the existence of a positive value of $x_9,$ say, $x_9^{\ast }$ satisfying $Q_2\left(x_9^{\ast }\right)=0.$