EXAMINING HIV PROGRESSION MECHANISMS VIA MATHEMATICAL APPROACHES

The progression of HIV infection to AIDS is unclear and under-examined. Many mechanisms have been proposed, including a decline in immune response, increase in replication rate, involution of the thymus, syncytium inducing capacity, activation of the latently infected cell pool, chronic activation of the immune system, and the ability of the virus to infect other immune cells. The significance of each mechanism in combination has not been studied. We develop a simple HIV viral dynamics model incorporating proposed mechanisms as parameters that are allowed to vary. In the entire parameter space, we derive two formulae for the basic reproduction number (R0) by considering the infection starting with a single infected CD4 T cell and a single virion, respectively. We then show that both formulae are equivalent. Analytically, we derive conditions for the occurrence of backward and forward bifurcations. Numerically, we perform uncertainty and sensitivity analysis to identify model parameters that significantly affect disease progression, which include the infection rate, infected cell death and proliferation rates, and viral production and clearance rates. Focusing on these identified significant parameters, a series of numerical bifurcation analyses demonstrate various HIV/AIDS progression dynamics through one or two slowly changing parameters.


Introduction
HIV can infect all cells in the immune system and the central nervous system which have a CD4 receptor on the cell surface, including T helper cells, monocytes, macrophages, and dendritic cells. However, the main target of the HIV virus is the CD4 T-helper lymphocyte, the main driver of the immune response. A large reduction in the number of CD4 T-helper cells seriously weakens the immune system, affecting the ability to fight opportunistic diseases. When the CD4 T-cell count reaches a measurement below 200 cells per microlitre of blood, a patient is diagnosed with AIDS. Treatment with antiviral therapy can delay the onset of AIDS, but this depends on drug adherence and the evolution of drug resistance.
The progression of HIV to AIDS is marked by a decrease in the CD4 T-cell count, and an increase in the viral load [74]. The mechanism by which HIV infection transforms into AIDS is unclear [15]. Several factors such as a decline in immune response [23,22], increase in replication rate [8,41,60], involution of the thymus [75,7], syncytium inducing capacity [23,40], activation of the latently infected cell pool [27,20,58,14,59,10], chronic activation of the immune system [51,38,61,35,30], and the ability of the virus to infect other immune system cells [50,70,67,19] have been associated with HIV progression to AIDS. However, it is not known which factors most drive disease progression. Mathematical modelling is well suited to provide insight into this complex process [4,47,49].
Mathematical modelling studies of HIV infection in-host are abundant in the literature. The vast majority of these studies, however, focus on the acute and latent stages of infection. Mathematical models of HIV progression to AIDS are presented by [73], [1], [24], [74], [17], and [34]. While these models have been successful in presenting T-cell count and viral load curves which demonstrate HIV progression to AIDS, these studies have fallen short of determining what characteristics of the underlying biology most drive the decline in CD4 T-cell population and increase in viral load. This may be due to the fact that many of these studies include models composed of a large number of equations, making model analysis very challenging. The majority of these studies also include a large number of unknown parameters, therefore identifying model parameters linked to specific biological processes which drive HIV progression can be difficult.
A recent study [74], which focused on the effects of low level homeostatic proliferation and/or immune activation on T-cell decline, concluded that a single mechanism cannot explain the slow decline in T-cell population size. This motivates our current study.
We have developed a mathematical model of HIV infection in-host progression to AIDS consisting of three ordinary differential equations. The model includes the effects of thymic involution, density dependent proliferation of CD4 T-cell population, and a growth term in the productively infected cell pool signifying the addition of infected cells from other infected cell pools including the latently infected cell pool. The model is presented in Section 2. Most of the analytical results are presented in Section 3. The basic reproduction number R 0 is derived using two different approaches which are shown to be equivalent. Conditions for the occurrence of backward and forward bifurcations are derived. Sensitivity analyses on the endemic equilibrium are presented in Section 4 to rank the influence of all model parameters, and identify those that most drive HIV progression to AIDS. In Section 5 we further study the identified parameters through numerical bifurcation analyses and simulations. Finally, we draw our conclusions in Section 6.

Model
The model consists of three ordinary differential equations describing the uninfected and productively infected target cell populations, and the infectious HIV viral load. The model represents a simple extension from the basic model of virus dynamics by [49], including density dependent proliferation of the T-cells, and time dependent production rates of the target cells. We let x, y, and v respectively represent the number of uninfected T-cells, infected T-cells, and infectious virus particles at time t. The model is as follows: All parameters are nonnegative. We assume that they are nonzero unless explicitly stated otherwise. We now describe each equation in more detail: Equation for x: CD4 T-cells are produced by the thymus. We assume that this is affected by two mechanisms: natural production that can vary by the age of the individual, and thymic involution [75,7], which depends on the viral load. We assume that T-cell production will decrease as the viral load increases [28,42,48,52,65] given a saturating constant (i.e., HIV induced involution of the thymus will not occur immediately upon infection with the HIV virus [25]). We will examine the influence of the parameter λ. To maintain homeostasis in T-cell count in uninfected hosts, T-cells must proliferate to account for the decreasing thymic production (homeostatic proliferation) [64]. This is modeled by logistic growth. In the absence of population limitation, the average uninfected and infected T-cell proliferation rates are p and r, respectively. In the presence of T-cells, the proliferation rate reduces to zero as the total T-cell load (x + y) approaches the carrying capacity T . This density dependent proliferation term implies that the increase of the infected cell population negatively affects the healthy T-cell population. Similar to the basic model of virus dynamics [49], uninfected cells either die (rate d) or become infected by HIV infection (βxv).
Equation for y: Productively infected CD4 T cells are produced by the infection of uninfected T cells by virus (βxv). We also assume that these cells can proliferate [13], and that latently infected cells can enter the productively infected cells pool. This is all captured by the term ry(1 − (x + y)/T ). Infected cells die at a rate of a. Note that this model is similar to a model presented by [58], where the latently infected cell pool was considered. In our study we consider changes in r, a, and β to reflect changes in the infected cell population, immune system killing of infected cells (immune activation), and changes in infection rate (which reflects the use of drug therapy), respectively. Equation for v: Infectious virus particles are produced from productively infected CD4 T cells and from other types of productively infected cells including monocytes, macrophages and dendritic cells [69,12]. The production of infectious free virus is represented as kq, where k is the number of infectious virions produced by a productively infected CD4 T-cell, and q gives the ratio of the total infected cell pool including CD4 T cells, monocytes, macrophages and dendritic cells to the productively infected CD4 T-cell pool. Thus, when q = 1 only productively infected T-cells produce virus, but as q increases, the production of the virus becomes more dependent on productively infected cells of other cell types [39]. This allows us to account for a shift in infection to monocytes, dendritic cells and macrophages when there is a limited pool of CD4 T cells to infect [26]. Note that we have elected not to explicitly model these cell types so as to make the analysis of the model feasible. Infectious virus is lost due to the infection of an uninfected cell (βxv), or is cleared by the immune system (uv). We note that the loss of virus particles in the infection process is ignored in the vast majority of models of HIV infection in-host. However, it has been found that this term can play a role in the infection dynamics [37,32], T-cell count and viral load [9,29,31]. Therefore, this term may be important in studying the progression to AIDS. In our study we consider changes in k, q, u, and β to reflect changes in virus production (including drug therapy), the ratio of productively infected T cells to other types of productively infected cells, virus clearance (immune activation), and cell infection (drug therapy), respectively.
In summary, the progression of HIV to AIDS may involve a few or many model parameters. For example, a decrease in the production of the CD4 T cells (λ), which represents an aging thymus, or thymic involution, may contribute. Also, it is possible that a change in the death rate of infected cells (a), in the clearance rate of the virus (u), in the production of the infected cells (r), or in the production of infectious virus particles (k) may affect disease progression. A shift from the infection of CD4 T cells to monocytes, macrophages and dendritic cells will also play a role (q). We therefore assume that all model parameters may change slowly with progression of the infection. This assumption allows us to study the key underlying mechanisms related to HIV progression to AIDS. That is, we assume that all model parameters can vary with time, but that this change is gradual with respect to time. We formulate a model of ordinary differential equations with parameter values varying in feasible ranges. In this paper, we provide analytical and numerical analyses that aid in determining key parameters related to this progression process.

Basic Properties of Solutions.
3.1.1. Well-posedness of the Solutions of Model (1). We first remark that Equation (1) is reasonable in the sense that no population size can be negative, and no population grows unbounded. For this well-posed problem, we first prove positivity and then further show boundedness of the solutions of model (1). The result is written in Theorem 3.1.
Theorem 3.1. Assume that initial conditions for system (1) are non-negative. Then the solution of system (1) exists and is nonnegative for all t ∈ [0, ∞). Furthermore, the solution is unique and bounded.
The nonnegativity of the solutions follows. Next, we prove that every nonnegative solution of model (1) is attracted to the bounded set Then, we calculate the derivative of (x + y) to obtain It follows that lim sup t→+∞ (x + y) ≤ B. As to the virus population, we obtain that from which it follows that lim sup t→+∞ v(t) ≤ kq(B + 1)/u. The proof is finished.

Steady States.
To get the equilibrium solutions of model (1), we set the right sides of the three equations equal to zero as follows: Firstly, the third equation in Equation (2) gives Substituting the preceding expression into the second equation of Equation (2) yields The first factor from Equation (4) gives a disease-free equilibrium (DFE), E 0 = (x 0 , 0, 0). Here, x 0 is the non-negative root of Since all parameter values are taken positive, there is a unique positive root of Equation (5), namely The second factor from Equation (4) being zero yields v as a function of x, that is Moreover, assuming that y > 0 and v > 0, we can rewrite the second equation of Equation (2) as Substituting the above expression into the first equation of Equation (2) shows that any infected equilibrium (x, y, v) must satisfy λ = H(x), where we define We conclude that the following assertions hold for any given positive parameter values.
(a) Every infected equilibrium is of the form and satisfies the equation λ = H(x).
We remark that the generating rate of naive T cells from thymus λ is the critical measurement monitoring the thymus involution during HIV infection. Moreover, the equation λ = H(x) must be satisfied for any infected equilibrium. This relates a parameter that changes over age to parameters that change in disease progression. Given this important relationship, we further analyze model dynamics with respect to the change of λ through a sensitivity analysis.
3.2. The Basic Reproductive Ratio R 0 . The basic reproductive ratio R 0 is defined to be the number of infected cells (or virions) produced by a single infected cell (or virion) when the infection is introduced into a totally susceptible population of target cells (i.e., the DFE E 0 : (x 0 , 0, 0)). In the present case, two different expressions for R 0 can be derived, depending on whether infected cells or virions are considered; however, the threshold condition "R 0 > 1" for instability of the DFE is valid for both expressions.
3.2.1. Deriving R 0 from a Single Infected Cell Gives R y 0 . First, we consider R y 0 , the number of "next generation" infected cells produced by a single infected cell introduced at the DFE. There are two ways that an infected cell can give rise to new infected cells: (i) by cellular proliferation, and (ii) by producing infectious virions which then infect uninfected cells. For (i), a single infected cell proliferates at rate r(1 − x 0 /T ); since the typical lifetime of an infected cell is 1/a, we see that a single infected cell produces a total of r(1 − x 0 /T )/a new infected cells by proliferation. For (ii), a single infected cell produces a total of kq/a infectious virions during its lifetime, and each of these infectious virions has probability βx 0 /(u + βx 0 ) of infecting a healthy cell before the virion dies. Therefore, we obtain R y 0 > 1 is the threshold of the emergence of the infected equilibrium E 1 . Note that if r = 0, then R y 0 is similar to definitions of R 0 found in previous studies which include the viral loss due to infection in the virus equation in Equation (1) [33, 31, 32].

3.2.2.
Deriving R 0 from a Virion Gives R v 0 . Next we consider R v 0 , the number of "next generation" infectious virions produced by a single infectious virion introduced at the DFE. This is more subtle, because we need to count possibilities, such as, a virion (called v 1 , say) infecting a cell (called c 1 ), which proliferates to create a new infected cell (c 2 ), which in turn proliferates to create a new cell (c 3 ), and this cell c 3 produces a new virion v 2 . Although the line from v 1 to v 2 involves several generations of infected cells, it does not involve any other virions; in that sense, we consider virion v 2 to be part of the "next generation" of virions produced by virion v 1 . Following this understanding, we define the quantity N y to be the total number of infected cells descended from a single infected cell by proliferation only (i.e. over all generations but excluding any new infections by virions), including the initial infected cell. If r = 0, then N y = 1. If r > 0, then N y satisfies the equation because 1 represents the initial cell and r(1−x 0 /T )/a represents the number of first-generation offspring of a single infected cell during its lifetime (as argued for R y 0 above), and each of those first-generation offspring produces a total progeny of N y over all generations by proliferation alone.
Since the number of infected cells is negligible in this scenario, we only include the uninfected cells x 0 in the density dependence term r(1 − x 0 /T )/a.
Observe that if a ≤ r 1 − x0 T , then the only nonnegative solution of Equation (11) is N y = +∞. This makes sense because the inequality a < r 1 − x0 T tells us that when y is small, the death rate of an infected cell is less than its proliferation rate, and so the number of infected cells increases. Thus, the population of infected cells will never die out, even if the immune system is 100% effective at immediately killing all free virions; hence N y is infinite. In contrast, if a > r 1 − x0 T , then we can solve Equation (11) to obtain the finite value Finally, we have where we should interpret the final expression to be +∞ if it is negative or undefined. Equation (13) arises because each infectious virion has probability βx 0 /(u + βx 0 ) of infecting a new cell, which in turn would produce a total of N y infected cells by proliferation alone, and each of these infected cells would produce kq/a infectious virions. Note again that, if r = 0 at the beginning of infection, Equation (13) is similar to definitions of R 0 found in previous studies which include the viral loss due to infection in the virus equation in Equation (1) by [33,31,32].
Simple algebraic manipulations show that R v 0 < 1 if and only if R y 0 < 1. Given the relationship of parameters a and r discussed above, we mark these parameters for further investigation later in this study. (1), the Jacobian matrix associated with its linearized system evaluated at E 0 = (x 0 , 0, 0) takes the form

Linear Analysis for DFE (E 0 ). Regarding model
Its corresponding characteristic polynomial takes the form as The positive values of all parameters and x 0 (that is T ) assuming that a non-existent infected cell population is possible. The first factor of Equation (15) gives a negative real root (or eigenvalue). Roots derived from the second factor of Equation (15) take the form of L 1, 2 = 1 2 −a 01 ± a 2 01 − 4a 02 . The local stability of E 0 is then determined by the sign of the real parts of L 1 and L 2 .
Proof. We rewrite a 02 in (15) as from which the lemma is immediate.
The three roots of the characteristic polynomial P (L) = 0 are as follows: (i) If a 02 > 0: three negative real roots, or one negative real root and a pair of complex conjugate roots with negative real parts; (ii) If a 02 = 0: two negative real roots and one zero root; (iii) If a 02 < 0: one positive and two negative real roots.
In particular, no root can be purely imaginary.
Proof. Since the characteristic equation P (L) = 0 has one negative root and the other parameters satisfy Lemma 3.2, the proof is obvious. On the other hand, the well-posed property of the solutions of model (1) guarantee that the solutions stay positive with positive initial conditions. Since the DFE E 0 is located on the x-axis, it cannot be enclosed by any periodic solution.
, and a 02 = 0 are equivalent thresholds to determine the stability of the DFE.
Proof. The condition for a stable equilibrium E 0 is Similar argument proves that the condition for an unstable equilibrium E 0 is a 02 < 0 ⇔ R v 0 > 1 ⇔ R y 0 > 1 and the static bifurcation condition is equivalent to the two thresholds, that is a 02 = 0 ⇔

3.4.
Nonlinear Analysis for DFE (E 0 ). Theorem 3.4 shows that the DFE has a zero eigenvalue when a 02 = 0 (or R 0 = 1), which is equivalent to The other two eigenvalues, −a 00 and −a 01 shown in Equation (15), are both negative. We choose right and left nullvectors corresponding to the zero eigenvalues as We obtain < w, v >= 1. That is the inner product (< ·, · >) of the column and row vectors (v and u) is one. Then we have B > 0 is satisfied for all positive parameter and x 0 values. Moreover, we denote Here, A and B in Equation (19) and Equation (20) are coefficients of the center manifold of model (1) at the DFE E 0 when R 0 = 1. The center manifold up to third order is written aṡ Applying the results from [68] and [11], we have the following result.

Uncertainty and Sensitivity Analysis
In the previous analysis we identified three parameters that were critical to determining the existence of the infected equilibrium and the stability of the disease-free equilibrium, namely, λ, r, and a. In the following we determine how changes in these parameters affect the infected equilibrium.
4.1. Analytical Sensitivity Analysis. Suppose that we want to say how a stable infected equilibrium changes as one parameter varies, with all other parameters being held constant. For example, Theorem 4.1(i) below says that the x-coordinate of the fixed point is an increasing function of λ. Before we state the results formally, we shall show that the Implicit Function Theorem guarantees that slightly varying a single parameter leads to a differentiable curve of stable fixed points passing through a given stable infected equilibrium.
Suppose (x 0 , y 0 , v 0 ) is a stable infected equilibrium corresponding to a particular choice of strictly positive parameter values (λ 0 , 0 , d 0 , β 0 , . . .). (Observe that x 0 , y 0 , and v 0 are necessarily strictly positive.) We shall focus on λ for concreteness, but the same argument works for any parameter. Define the function G : (0, ∞) 4 → R 3 by Then the fixed point equations used to determine the solutions to Equation (2) are equivalent to writing G = 0. Write the 3 × 4 derivative matrix DG in the block form where J is the 3 × 3 Jacobian matrix (14) (except that now we treat λ as a variable) and w is the column vector By our definition of stability, det(J) < 0 at this equilibrium. Thus the Implicit Function Theorem guarantees the existence of a δ > 0 and a continuously differentiable function Φ : (λ 0 − δ, λ 0 + δ) → (0, ∞) 3 such that Φ(λ 0 ) = (x 0 , y 0 , v 0 ) and the equation G(Φ(λ), λ) = 0 holds for every λ ∈ (λ 0 −δ, λ 0 + δ), where we write Φ(λ) = (x(λ), y(λ), v(λ)); that is, Φ(λ) is a curve of fixed points parametrized by λ. Since the eigenvalues are continuous functions of the coefficients of the characteristic polynomial (e.g. see Theorem 1.4 of [43]), we can choose δ small enough so that the fixed point Φ(λ) is stable for every λ in (λ 0 − δ, λ 0 + δ). In particular, the derivatives in the following theorem all exist, and we can safely interpret dx dλ (λ 0 ), for example, to be the rate of change of x with respect to λ at the fixed point (x 0 , y 0 , v 0 ). Proof. We first consider varying λ. Differentiating the equation G(Φ(λ), λ) = 0 with respect to λ gives J(dΦ/dλ) + w = 0 (where we write dΦ/dλ as a 3 × 1 column matrix), or equivalently Let K ij be the ij entry of J −1 . Then Equations (22)(23) show that (omitting the 0 subscripts on the parameters) dx dλ Writing J ij as the ij entry of J, matrix algebra says that At an equilibrium point, the diagonal elements of J given in Equation (14) can be rewritten as follows. Using the first equation of Equation (2), we obtain From the second equation of Equation (2), we obtain And using the last equation of Equation (2), we obtain Then the first equation in Equation (24) yields We assume that (x, y, v) is a positive stable equilibrium, then obtain det(J) < 0 and dx/dλ > 0. This proves (i). The second equation in Equations (24) yields, (28)).
It follows that d(x + y) dλ = − βvu (v + ) det(J) which is positive whenever (x, y, v) is a stable infected equilibrium. This proves (ii).
The proofs for varying a and r are very similar to the arguments for λ, with the following differences. The vector w is now  The resulting analogues of Equation (24) are We then have from Equation (14) that which has negative numerator and denominator, and, with the help of Equation (26), which has positive numerator and negative denominator. These suffice to explain parts (iii), (iv), (vi), and (vii). Finally, parts (v), and (viii) follow from the relation v/y = kq/(u + βx) (from the fixed point equations of (2)) and parts (i), (iii), and (vi).
The correlation between the infected T-cell proliferation rate r and the uninfected T-cell count (x), infectious free virus (v), and ratio of infectious free virus to infected cell count (v/y) are all determined by the sign of the difference between the total number of uninfected and infected T-cells (x + y) and the carrying capacity (T ). The magnitude of the correlation, which is difficult to ascertain in the analysis, will be demonstrated by the numerical sensitivity analysis results which are presented in Figure 1 in Section 4.2. Moreover, Figure 1 also confirms results (i)-(v) in Theorem 4.1.
The progression of HIV to AIDS is marked by a decrease in the CD4 T-cell count and an increase in the viral load. The results of Theorem 4.1 show that decreases in CD4 T-cell count are feasible when λ is reduced, a is increased, and r is increased (but only if sgn(x + y − T ) < 0). These results also show that the viral load increases when a is decreased and r is increased (but only if sgn(x + y − T ) < 0).

4.2.
Numerical Sensitivity Analysis. We now further explore changes in the infected equilibrium with respect to all parameter values through uncertainty and sensitivity analysis using Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCC) [46,5]. LHS gives adequate quality assurance on model predictions and PRCC determines the main parameters driving infection.
The LHS is performed using the numerical algorithm described in [44]. We adopt uniform distributions for model parameters with parameter value ranges from the modeling and clinical literature (see Table 2). Parameter values are chosen randomly without replacement. We only record parameter sets satisfying four filter conditions to ensure that equilibrium solutions lie within realistic ranges for HIV/AIDS (see Table 1). Following this procedure, a PRCC value is calculated for each model parameter. PRCC values range between -1 and 1 with the sign determining whether a model parameter is negatively (-) or positively (+) correlated with the specified model output. A PRCC value |P RCC| > 0.5 was considered statistically significant based on empirical knowledge. We chose 100, 000 parameter sets and 910 of these satisfied the filter criteria (Table 1). Finally, within the 910 parameter sets, 6 parameter sets produced oscillating solutions for Equation (1). A further study of those oscillating behaviors is carried out in the next section. The occurrence of periodic solutions agree with the result from [78] suggesting that backward bifurcation implies rich dynamical behaviors.
The PRCC values relating the model parameters to the infected equilibrium are shown in Figure 1. It suggests that parameters r, a, u, β, and k q significantly affect population sizes of the uninfected and infected T-cells, and the viral load. Alternatively, parameters d, p, , and λ do not significantly affect the infected equilibrium.
In terms of the CD4 T-cell count, sensitivity analysis results suggest that an increase in the uninfected CD4 T-cell count (x) and a reduction in the infected and total CD4 T-cell count (y and x + y) can occur as the infected cell death rate (a) increases and proliferation and/or addition to the infected cell pool (r) decreases. Moreover, Equations (12) and (13) confirm this PRCC result, that is the relation between a and r is important to determine the total CD4 T-cell count (x + y) and in turn determine the fate of the disease.
The parameter influence ranking on the infectious viral load (v) is shown in Figure 1 (left column, middle row). It shows that the most influential factors are the viral clearance rate (u) and the viral production rate (kq), followed by the rates of infected cell death (a) and proliferation and/or addition to the infected cell pool (r). Here, increases in kq and r, and decreases in u and a increase the viral load.
The ratios of the infectious free virus to the infected cell count (v/y) and the total CD4 T-cell count to the infectious free virus ((x + y)/v) are shown in the bottom row of Figure 1. Here we see that u and kq are the only parameters that significantly affect v/y and (x + y)/v. In fact, we see that a decrease in u and an increase in kq both increase the viral load, and decrease the T-cell count, and therefore significantly affect HIV progression to AIDS.
The influence of the infection rate (β) has direct implications on HIV drug therapy, where the goal is to decrease viral load and increase the CD4 T-cell count. Here, we see that a reduction in β can increase both x and y and decrease v (although this is not in the significant range of the sensitivity analysis for y and v).
In summary, the LHS-PRCC results suggest that an increase in CD4 T-cell count and a decrease in HIV viral load are significantly influenced by the death rate of the infected cells killed by immune response (a) and the proliferation rate and addition of infected macrophages, monocytes and dendritic cells and activation of latently infected cells (r), but that HIV progression to AIDS is mostly determined by viral clearance rate (u) and the viral production rate (kq). The production rate of uninfected CD4 T-cell (λ) is not a significant factor to model dynamics. Note that the above conclusions are drawn from numerical sensitivity analysis results, which are valid in the chosen parameter regions infromed by literature review (see Table 2).

Bifurcation Analysis and Numerical Simulations
Based on the identified significant parameters from the sensitivity analysis results, we now focus our attention on a, r, kq, u and β.

Bifurcation Analysis on
Infected Equilibrium E 1 . In this section, we will further investigate how parameters r, kq, u, β and a affect disease progression through a bifurcation analysis of model (1). Saddle-node, transcritical and Hopf bifurcations are computed through symbolic computation. We first evaluate the Jacobian matrix of model (1) at E 1 , which yields the corresponding characteristic polynomial: where y and v take the form of (3) and (7) as Moreover, x satisfies F (x) = c 10 x 5 + c 11 x 4 + c 12 x 3 + c 13 x 2 + c 14 x + c 15 = 0, where, + ru 2 (dr − ap) − 4βλr 2 u), In the next subsection, we will numerically explore the possible bifurcations at the endemic equilibrium E 1 in terms of the above formulas and the following theorem from [76].

Numerical
Results. We carry out numerical bifurcation analysis based on Theorem 5.1. The results in Figures 2, 3, and 5 are verified by MatCont [21]. We first choose a and r as bifurcation parameters. Other parameter values are taken from Table 2. First, we set a as the control parameter and r as the bifurcation parameter. Model (1) shows a forward bifurcation when a = 2.9333 and backward bifurcations when a = 2.9438, a = 3.11, a = 3.5, and a = 3.8. The corresponding 1-and 2-dimensional bifurcation diagrams are shown in Figure 2(a) and (b). We obtain a limit point curve (plotted in green) consisting of Limit Point (LP) bifurcation points and a Hopf curve (H) connecting both the neutral saddle (in magenta) and Hopf bifurcation (in red) points at Bogdanov-Takens bifurcations (BT). Moreover, BT points are also intersections between the limit point and Hopf curves. Hopf bifurcation serves as an oscillation source, which induces the oscillating viral load. The stability of the bifurcating limit cycles are determined by the sign of the first Lyapunov coefficient of the corresponding Hopf bifurcation. A positive first Lyapunov coefficient indicates a subcritical Hopf bifurcation, which induces unstable limit cycles; while a negative first Lyapunov coefficient implies a supercritical Hopf bifurcation, which bifurcates stable limit cycles. A zero first Lyapunov coefficient represents an occurrence of a generalized Hopf bifurcation (GH), which separates subcritical and supercritical Hopf bifurcations in the Hopf curve. For the parameter values of a and r, indicating the death and proliferation rates of the infected CD4 T-cell, we plot a close-up 2-dimensional bifurcation diagram r vs a in biologically meaningful ranges in Figure 2(b). The blue curve represents a branching point (BP) or transcritical bifurcation, which is plotted according to R 0 (a, r) = 1. Here R 0 is the basic reproduction number. The top Hopf curve represents Hopf bifurcation because the vertical axis range is below the BT point in Figure 2(a), that is a < 3.114226. Therefore, for positive a and r values, two red Hopf curves and a blue transcritical curve enclose a parameter range plotted in yellow, in which model (1) demonstrates oscillating viral loads.
Three parameter values a = 2.9333, a = 2.9428 and a = 3.8, are chosen in the yellow region in the closeup of Figure 2(b). The corresponding 1-dimensional bifurcation diagrams are plotted in the left column of Figure 3. Oscillations occur between two supercritical Hopf bifurcations in Figure 3(a) and (b). The oscillation peak values are shown in terms of the uninfected CD4 T-cell population (x). The corresponding oscillation periods are followed in the right panel. Large viral load oscillations occur around r = 0.8 for both a = 2.9333 and a = 2.9428. For the case when a = 3.8, the transcritical bifurcation occurs at a very large r value. The DFE is locally stable in the biologically meaningful region of r. Numerical simulation in Figure 3(c) shows that the oscillation period obtained from the supercritial Hopf bifurcation approaches infinity. The infinite period oscillation may be due to the homoclinic cycle bifurcating from a local BT bifurcation. Therefore, large viral load oscillation may occur around r = 3 for a = 3.8.
Biologically, if the proliferation rate (r) and death rate (a) for the infected CD4 T-cell population are small, the uninfected CD4 T-cell population can stabilize/oscillate at a high level (see Figure 2(a)). As the infection progresses, and the infected cell proliferation rate (r) increases, the uninfected T-cell  population can then stabilize/oscillate at a low level. It follows that the increase of the infected CD4 T-cell death rate (a) shows a prohibitive effect on disease progression.
The oscillatory behaviors discussed here are of interest as they may indicate the existence of viral blips and provide the oscillatory base over which stochastic effects may generate large periodic viral loads. In Figure 3 we show model outcomes with a regular oscillatory behavior (see sub-figure (a) and (b)), and also cases that more closely reflect the occurrence of viral blips (see sub-figure (c)). We now  Table 2  . consider stochastic variability on the infected cell proliferation rate (r), modeled by a mean-reverting process [2,3]. The corresponding stochastic differential equation is written as where r e is an approximated steady-state infected cell proliferation rate, and α r is related to the changing speed of the proliferation rate. Moreover, lim t→∞ E(r(t)/r(0)) = r e and lim t→∞ V ar(r(t)/r(0)) = σ 2 r r 2 e 2αr−σ 2 r . Figure 4 shows a simulation of the stochastic system considering a moderate variability, i.e. σ 2 r < 2α r . We observe that the viral load oscillates with varying period and amplitude. A further exploration of a full stochastic model based on Model (1) is a course for future work.
Our analysis and simulations show that progression in the infected CD4 T-cell proliferation rate (r) and the death rate (a) can generate various stages of HIV progression, including low and high level infection, oscillating viral loads, and, perhaps, viral blips.
To complete our analysis, we now consider kq, u and β as bifurcation parameters, along with a and r. Figure 5 shows two parameter bifurcation diagrams for every parameter pair. These figures demonstrate the same outcomes as those previously discussed (saddle-node bifurcation, transcritical bifurcation, Bogdanov-Takens bifurcation, supercritical and subcritical Hopf bifurcations, and generalized Hopf bifurcation).

Discussion
Mathematical models describing HIV infection in-host offer a way to understand the dynamics of HIV during different disease stages. The vast majority of mathematical studies in the literature focus on the acute and latent stages of infection, ignoring the progression from HIV to AIDS. We have developed a mathematical model that includes biological mechanisms that have been associated with HIV progression to AIDS. These include: thymic involution, a reduction in T-cell production by age of an individual, density dependent logistic growth in the CD4 T-cell population, and the effects of the latently infected cell pool, immune system exhaustion, and contributions of production of free virus from other cells pools. It is currently unknown how much each of the processes listed above contribute to the health status of the host in the progression of HIV to AIDS. In the absence of prior data on the problem that we are considering, we determine the most significant parameters which can substantially change the model output behaviors through an uncertainty and sensitivity analysis, using a statistical based Latin hypercube sampling (LHS) and partial rank correlation coefficient (PRCC) analysis on all model parameters on clinically feasible parameter intervals. This uncertainty and sensitivity analysis provides the most influential parameters to different model outcomes of interest at the infected equilibrium. We find that the proliferation and death rates of the productively infected CD4 T-cell pool (r and a) most affect the total T-cell count (x + y), while the production rate of free virus from other cell pools and/or a change in production by productively infected CD4 T cells (kq) and the death rate of virus (u) significantly influence the total viral load. We also find that the ratio of total cell count to virus ((x + y)/v) is most affected by kq and u, and that HIV progression to AIDS is most affected by decreases in u (immune exhaustion), and increases in kq (activation and proliferation of latently infected cells and other infected cell types).
To show the fate of the disease, we carry out bifurcation analyses on a, r, kq, u, and β. Corresponding numerical simulations provide parameter regions within HIV progression stages, clearance, relapse, remission, and recurrence, can occur. Our numerical and analytical results not only can be related to specific biological processes, but also give the relative importance of each process.
Monocytes and dendritic cells have been implicated as HIV reservoirs using in vitro and ex vivo models of viral infection [16]. Monocytes and macrophages, after infection with HIV virus are resistant to cytopathic effects and persist throughout the course of infection as long-term stable reservoirs for HIV-1 which produce virus, and can disseminate the virus to tissues [36]. These cells have also been shown to contribute to the pathogenesis of HIV via the impairment of effector functions [36]. Although it is known that HIV interacts with monocytes, macrophages and dendritic cells, some key questions remain to be answered to fully understand the pathogenesis of HIV to AIDS. For instance, the relative contributions of these cell types in the development and persistence of the HIV latently infected cells reservoir, the activation of these cells, and the individual contributions of these cells in both viral and host aspects in the progression to AIDS remain to be elucidated. Mathematical models explicitly including these cells lines can contribute to this area of study. Our results confirm the importance of macrophages and dendritic cells to HIV progression to AIDS.