Bistable Dynamics and Hopf Bifurcation in a Refined Model of Early Stage HIV Infection

Recent clinical studies have shown that HIV disease pathogenesis can depend strongly on many factors at the time of transmission, including the strength of the initial viral load and the local availability of CD4+ T-cells. In this article, a new within-host model of HIV infection that incorporates the homeostatic proliferation of T-cells is formulated and analyzed. Due to the effects of this biological process, the influence of initial conditions on the proliferation of HIV infection is further elucidated. The identifiability of parameters within the model is investigated and a local stability analysis, which displays additional complexity in comparison to previous models, is conducted. The current study extends previous theoretical and computational work on the early stages of the disease and leads to interesting nonlinear dynamics, including a parameter region featuring bistability of infectious and viral clearance equilibria and the appearance of a Hopf bifurcation within biologically relevant parameter regimes.


1.
Introduction. Mathematical modeling of the in-host behavior of viral infections has become an indispensable tool to biological researchers in recent decades. New models have been used to describe the dynamical behavior of various infectious diseases such as HIV, HBV, and influenza, among others. Within this field, testing specific hypotheses based on clinical data is often difficult since samples cannot be taken frequently from patients, and viral load detection techniques may lack a necessary level of precision. Thus, new predictive models play a central role and are continually needed to further our understanding of disease dynamics. One such mathematical model that has been quite useful, known as the standard model of viral dynamics [23,26,27], describes the early stage in-host behavior of HIV infection. In general, the time course of this disease typically consists of three distinct phases.
The first, known as the acute stage, is characterized by a rapid fluctuation in both the healthy T-cell and virion populations, usually lasting from 2 to 10 weeks [13]. Symptoms during this stage, often described as "flu-like", include fever, swollen glands, sore throat, rash, and fatigue. During the second stage, known as chronic infection, the size of the uninfected T-cell population and viral load maintain relatively constant states, with the latter known as the viral set point. Without the aid of antiretroviral treatment, this period can persist for up to 10 years but can vary greatly among individual patients [9,24,29]. Finally, within the third stage the viral load experiences exponential growth with a correspondingly rapid decrease in the healthy T-cell population. This leads to the onset of AIDS, defined clinically as a T-cell count of an HIV-positive patient measured below 200 cells/mm 3 . While the standard model of viral dynamics has been extremely successful in reproducing the acute and chronic stages, it has been shown that both the infected and infection-free equilibrium states of healthy T-cells, infected T-cells, and virions induced by this model are globally asymptotically stable. This property was first investigated analytically in [5], and later proved using a Lyapunov function in [18]. The global stability of these states implies that the long time asymptotic behavior of the system depends only upon parameter values in the model. As a direct consequence, both the equilibrium values and stability properties of these equilibria are independent of initial conditions. However, a number of recent clinical studies [7,14,15,19] have shown that additional factors beyond these parameters may have a significant impact on the development or clearance of a persistent infection. Such factors include the initial viral load and the availability of target CD4 T-cells at the time of transmission. Because early events during infection may determine both the pathogenic consequences of the virus and its sensitivity to interventions or treatment strategies to combat the disease, Igarashi et al. [14] evaluated the effects of inoculum size on the development of the disease in macaques infected with Simian/Human Immunodeficiency Virus (SHIV). In particular, the results of this study showed that macaques who were administered large intravenous SHIV inocula experienced irreversible CD4+ T lymphocyte depletion and developed clinical disease. In contrast, rhesus monkeys receiving 50% tissue culture infective doses (or less) of virus survived the acute stage with reduced but stable levels of CD4+ T lymphocytes and produced antibodies capable of neutralizing SHIV. In short, although SHIV induced an extremely rapid and profound depletion of T-cells in all infected rhesus monkeys, the loss of this T-cell subset was not irreversible in animals inoculated with small amounts of virus. A similar investigation has since been conducted in [19] yielding analogous results.
A different study [15] has highlighted the importance of the T-cell count at the time of primary viral infection. In particular, during further studies of seventeen rhesus macaques, SHIV infection was found to emerge only in a single monkey whose T-cell count had been markedly depleted by monoclonal antibody (mAb) treatment at the time of primary viral infection, while none of the remaining sixteen monkeys inoculated with SHIV, but not treated with the mAb, developed immunodeficiency. A similar outcome was observed in [7]. Hence, the availability of target T-cells at the time of viral transmission can also play a large role in the establishment of a persistent infection, as differing strengths of the susceptible T-cell population may promote or inhibit viral replication.
Since the equilibria of the standard model of viral dynamics are globally stable in mutually exclusive parameter regimes, these empirical results suggest that, in order to appropriately describe early HIV (or SHIV) dynamics, additional factors must be considered in the model development. Other authors [1,3,6,10] have further posited such considerations, asserting the need for a variety of secondary biological characteristics including variability of host susceptibility to initial infection, withinhost competition between different viruses for target cells at the initial site of virus replication, and the effects of the innate immune response. These ingredients should play a realistic role in disease pathogenesis and long time dynamics. In general, previous in-host models do not account for effects arising from the strength of the initial viral load or variations in the T-cell count at the time of transmission, as they describe the tendency to viral infection or clearance based solely upon parameters and not on initial conditions. Another element of disease pathogenesis overlooked within the standard model is the homeostatic mechanism that regulates the peripheral T-cell pool. Recent clinical studies have displayed the importance that homeostasis of the susceptible T-cell population may play during infection [4,22], as the replenishment of target cells provides additional opportunities for HIV infection by freely moving virions. In the current study, we focus on the homeostatic proliferation of T-cells, i.e. the process by which T-cells in a lymphopenic host divide in the absence of cognate antigen to reconstitute the peripheral lymphoid compartment, which is believed to be driven by the presence of foreign antigens [20,31]. A few long-term models of HIV infection [12,13,25] have incorporated the homeostatic proliferation of the T-cell population within their formulation, but the dynamical effects of this biological mechanism are not well-understood. Other authors [8] have considered an acute stage model incorporating a logistic growth term, depending only upon healthy T-cells, to represent the body's propensity to regulate the T-cell population. However, in the setting of HIV-induced lymphopenia, it was determined [4] that the homeostatic proliferation of CD4+ T-cells is driven primarily in response to the viral load, while naive CD4+ T-cells are also recruited into the proliferating pool due to CD4+ T-cell depletion. Therefore, in the presence of HIV, such a regulatory mechanism should depend on the strength of the viral load in addition to the size of the T-cell pool.
Based on the aforementioned experimental findings concerning the influence of initial conditions and T-cell homeostasis, we explore a refined model of early stage infection dynamics that incorporates the ability of the immune system to maintain the T-cell count even when the number of such cells is depleted by the presence of the virus. In accounting for such effects, it will be shown that this model will accurately portray the dependence of equilibria on initial conditions by producing a biologically relevant parameter regime featuring bistability of the infected and uninfected equilibrium states. Hence, the model proposed herein will account for both of the aforementioned processes. In the next section, the new model of early HIV infection is discussed, and a study of parameter identifiability is conducted. In Section 3, we prove that exactly three states exist -one uninfected and two infected equilibria. In Section 4, the local stability properties of equilibria are characterized in terms of parameter values. In particular, we identify a biologically-important region of the parameter space within which both the relevant infected equilibrium and the viral clearance state are locally stable. This illustrates that the development of a persistent infection will depend crucially on initial conditions, and we further explore the basins of attraction generated by these equilibria. Finally, we show that the system experiences a Hopf bifurcation that gives rise to oscillatory behavior within a certain parameter regime. To conclude the paper, appendices containing proofs of the aforementioned results are provided.
2. Model and Parameters. The proposed dynamical model couples a nonlinear system of three ordinary differential equations given by Here, T (t) denotes the population of healthy T-cells, I(t) the population of these cells which have been infected, and V (t) the size of the virion population. The parameter λ represents the source of new cells arising from general production, while the healthy cell death rate is denoted by d T . The interaction, or mass action, term kT V , where k is the infection rate, represents the infection of healthy Tcells and the subsequent conversion of these cells to infected lymphocytes, with corresponding death rate d I . The parameter p is the rate at which new virions are created by the infected cell population, and the clearance rate of free virus particles is given by d V . See Table 1 for a complete list of parameters and variables with representative initial values. Additionally, in this model we do not consider distinct compartments within the host since the dynamics of interest take place over many weeks, while transfer between these compartments occurs on the time scale of hours. The term ρ C+V T V describes the homeostatic production of T-cells due to the presence of the virus and subsequent decline in healthy T-cells, both of which may vary over the course of infection. Here, ρ is the maximum growth rate and C is the half-velocity constant of growth. Note that the behavior of this term is limited by the growth and decay of the virus population. In particular, the function no virions are present in the system, this so-called Michaelis-Menten term vanishes and the basic dynamics are the same as the standard virus model. This is consistent with the actual immune response as the body need not further augment the T-cell population in the absence of virions. Contrastingly, as the virus population grows large, the infected host's immune system replenishes the T-cell population so as to balance the effects resulting from its depletion, and this occurs at a growing rate whose maximal impact is ρT . Regardless of the limited rate of growth within this term, the inclusion of homeostatic proliferation, as we will show, has a profound affect on the dynamics of the system. While the new model (3CM) can be derived from a bottom-up approach merely by adding the homeostatic proliferation term to the standard model of viral dynamics, it also stems directly from a top-down approach. More specifically, (3CM) can be fully derived from a reduced description of long-term models that were proposed in [12,13] to accurately represent all three stages of HIV infection within a host. In particular, a dynamic active subspace decomposition of the twenty-seven dimensional parameter space within the three-stage model of [12], which features seven different in-host populations, was performed in [25]. This decomposition produces a global sensitivity analysis of the parameter space and indicates exactly which parameters are important to the evolution of the model during each of the three distinct phases of disease progression. Upon eliminating those parameters (of which there  [13] Dimensionless Populations ( * omitted in exposition)  [12,13] Dimensionless Parameters Table 1. Variables and Parameters were nineteen) that are found to be negligible throughout the acute stage, a total of four populations -namely the influence of latently-infected T-cells, macrophages, infected macrophages, and the cytotoxic lymphocyte response -completely decouple from the model. Hence, the reduced system (3CM) results, providing a more precise description of the early stage behavior of the disease than the standard viral dynamics model. Figure 1 contains a representative simulation of (3CM) and includes a comparison to the early stage behavior of the long-term model of [12]. We note that other models of HIV infection [2,13,30] have also incorporated such a Michaelis-Menten term to describe homeostasis, though the current article will contain the first dynamical analysis of such a model. 2.1. Dimensionless system. To reduce the size of the parameter space, the original model (3CM) is recast in dimensionless form. The resulting system, in which dimensionless populations have been renamed T * , I * , and V * is  Table 1, and a comparison with the full three-stage model in [12] -T-cell count (left) and viral load (right). where The values of dimensionless parameters are summarized within Table 1, and in the future we will remove the * notation and deal solely with the dimensionless system. Notice that each new parameter is positive since the original variables are positive. The complete derivation of (3CM*) from (3CM) can be found in Appendix A. The model (3CM*) contains only five parameters, and each may play a role in the dynamics of the system. However, we will typically fix the values of α 1 , α 2 , and β while considering variations in R 0 and R m , which represent the usual basic reproduction number (as in the standard viral model) and a new reproduction number generated by the addition of the Michaelis-Menten term, respectively.

2.2.
Parameter Identifiability. With the dimensionless system determined, we study parameter identifiability in (3CM*) as this model can provide useful simulations only if the parameters involved can be discerned from data. In particular, we first conduct a test of (3CM*) developed for differing models in [21], [32], and [33] to understand the structural identifiability of parameters. Structural identifiability is used to characterize the one-to-one property of the map that takes the parameter space to the set of system outputs (i.e., the information encapsulated by collected data). In order to evaluate this property for (3CM*) we use the so-called Multiple Time Points (MTP) method developed in [32], which entails the construction of an invertible identification function Φ(θ), from the parameter space to the set of observable outputs, that preserves the structure of the differential equations model. In particular, because invertibility of such a mapping is required, we wish to ultimately conclude that ∂Φ ∂θ has full rank. To begin, we first describe the space of output values. Because healthy and infected T-cell counts can be both difficult to measure and unreliable, data is most easily gathered from an individual's viral load. Hence, model outputs in this context will be regarded as values of the viral load and its derivatives, as the latter are needed to compensate for the lack of T-cell data but can be generated from values of V . Thus, we begin to construct Φ by first eliminating the populations T and I within (3CM*) in favor of derivatives of V . This procedure involves merely taking derivatives in (3CM*) and representing T and I in terms ofV ,V , and ... V and yields a single equation to represent the original three-dimensional system of ODEs, namely .. .
where f is given by (3) below and θ = (R 0 , R m , α 1 , α 2 , β) T is the vector of parameters. Hence, any solution of (3CM*) can be characterized by satisfying this relationship at time t.
In order to identify the five distinct parameters in the model, five identification equations are needed, and this requires us to satisfy the above ODE at five different time points, say t k , for k = 1, ..., 5. Given this, we denote V k = V (t k ), with the same notation for derivatives (e.g.,V k =V (t k )), and construct the identification function Φ : By construction, the model (3CM*) is trivially satisfied at fitted parameter values θ * given in Table 1 as Φ(θ * ) = 0, and we are thus interested in whether Φ is invertible for values of θ ≈ θ * . Now, computation of Φ requires knowledge of V and its derivatives at five different time points, but values for these derivatives are typically unavailable either as collected data or via direct simulation of the model. Thus, we must require additional values of the viral load in order to compute them. In particular, eight values of V are needed to numerically approximate ... V k for k = 1, ..., 5 by using a suitable finite difference approximation. Therefore, we choose three additional time values t 6 , t 7 , and t 8 at which V must be known, and note that the values ofV k ,V k , and ... V k for k = 1, ..., 5 are merely determined by values of the viral load at multiple time points; for example, ... V k depends upon V 1 , ..., V 8 for every k = 1, ..., 5. In order to conclude that the model is locally structurally identifiable, we must show that the corresponding Jacobian matrix ∂Φ ∂θ is invertible near the fitted values θ = θ * . Since this is nearly impossible to perform analytically, we instead take a computational approach. First, we symbolically represent the matrix ∂Φ ∂θ using (2) and (3). Then, fixing a specific vector of parameter values θ, we simulate the output variable V (t; θ) at a chosen sequence of times using (3CM*) and Matlab's ode15s solver, and then numerically approximate its derivatives. Finally, we use these simulated values of V k ,V k ,V k , and ... V k to compute the resulting rank of the Jacobian for these particular parameter values. Repeating this calculation over a grid of parameter values within the biologically reasonable ranges [0.5θ * k , 1.5θ * k ] for k = 1, ..., 5, we find rank ∂Φ ∂θ = 5 for every simulation. Thus, we conclude that ∂Φ ∂θ 8 STEPHEN PANKAVICH, NATHAN NERI, AND DEBORAH SHUTT is of full rank and the associated parameters are structurally identifiable, at least locally, in the range of parameter values used for these simulations.
While this analysis provides a theoretical assurance that parameter values can be identified from exactly observed viral load data, clinical measurements will always contain some level of error. Even a model such as (3CM*), in which parameters can be uniquely identified locally within the parameter space, may yield unreliable parameter estimates due to noisy fluctuations or measurement error in the data. Hence, we also study the practical identifiability of parameters, namely the relative proximity of fit parameters to their true values given uncertainty within obtained data, by using a Monte Carlo method outlined within [21] and [32].
To estimate the differences in parameter fits generated from noisy data, we will use a metric known as the average relative estimation error (ARE). Prior to precisely defining this quantity, we outline the algorithm for generating such values. Beginning with the previously fit vector of parameters θ * , we first use the numerical ODE solver to generate a time course of baseline viral load values V * j = V (t j ; θ * ) for a chosen set of times t j , j = 1, ..., M . Next, we choose a sensitivity threshold δ > 0 and a number of Monte Carlo trials N ∈ N, then define the perturbed viral load dataV where ǫ n j ∼ N (0, δ) represents an unbiased normally-distributed measurement error for each fixed j and n. With random error introduced within the simulated data, we perform N parameter fits of this data to generate N new vectors of parameter valuesθ n , for n = 1, ..., N . Finally, to evaluate the variations in these fits generated by the noise, we define the ARE for each parameter by whereθ n k is the estimate of the kth parameter ofθ n arising from the nth perturbation with variance δ. Hence, this Monte Carlo method simulates the introduction of Gaussian measurement noise within the viral load data (4) based on the output model (3CM*) and computes the expected response in parameter values from these variations.
The ARE algorithm was applied to simulations of N = 1000 distinct simulated noisy measurement sets for each δ noise level of 5, 10, 15, ..., 30 percent of the fit value θ * k , at time points t = 0, 1, 2, ..., 90. The results are summarized in Table 2. Hence, we find a collection of small relative errors, with only the error in β rising above δ% of the true value. That being said, β also displays relatively minor fluctuations throughout the simulations, even for noisy data, and other parameters possess even less deviation from their fit values. Thus, (3CM*) appears to be quite robust with respect to variations in measurement data for the purposes of parameter fitting.
3. Steady States. We begin an analysis of the dynamics of (3CM*) by first determining all steady states and investigating their regions of biological relevance within the parameter space. This information will be used extensively in the next section in which the local dynamics of solutions is characterized. We first compute the associated steady states, which are given in the form of an ordered triple (T, I, V ) T .  In particular, we find exactly three solutions to the algebraic system guaranteed by and they are summarized within the following theorem.
Theorem 3.1. The only time-independent solutions of (3CM*) are where a, b, and c are defined by and the dimensionless parameters are defined by (1).
Here, E u is the uninfected steady state while E + i and E − i represent states of persistent infection. In the future, when referring to components of equilibria, we will use a bar to distinguish between the components of these states (e.g., V ) and timedependent solutions (e.g., V (t)). Additionally, we will distinguish amongst the same components of different equilibria using subscripts (i.e., V u , V + , and V − ). Since all components represent scaled population sizes, we impose restrictions on the values for which these steady states are biologically reasonable. All three populations of E u will remain nonnegative for all times. However, for both E + i and E − i , we must require that the infected T-cell population, I, and the virus population, V , be real and positive. To ensure real valued populations, we impose the restriction b 2 − 4ac ≥ 0, which is equivalent to the condition The requirement that all populations of the infected states be positive forces other restrictions. For positivity of the E + i state, we must impose either or the condition For all populations within the E − i state to remain positive, we must impose the condition These constraints are justified within Appendix B. Additionally, Figure 2 provides a graphical summary of the restrictions on parameters necessary to guarantee positivity of corresponding equilibria. We will often refer to such a region as the "region of existence" of an equilibrium state, and in studying equilibria, we will always assume that parameters are within the region of existence of the state under consideration. Clearly, these restrictions depend upon only three parametersβ, R 0 , and R m . However, because β does not greatly affect the qualitative structure of the system, we will fix this parameter and focus on the behavior of the system depending only upon R 0 and R m . A similar approach will be taken in the investigation of stability properties of equilibria, which may further depend upon α 1 and α 2 , but these two additional parameters are ratios of death and clearance rates, which are fairly well-known. Thus, we will later fix these parameters as well, and again focus on the relationship between R 0 and R m .  Table 1.
Theorem 4.1. If R 0 < 1 then the infection-free equilibrium E u is locally asymptotically stable, whereas if R 0 > 1 then it is unstable. Additionally, for all parameter values that guarantee the positivity of the components of E − i , this equilibrium is unstable. Finally, let V + denote the value of the viral load for E + i , which can be expressed in terms of β, R 0 , and R m by Theorem 3.1. Then, the equilibrium E + i is locally asymptotically stable if the condition is satisfied, and otherwise unstable. These parameter regimes are summarized by Figure 3.
In short, only the uninfected steady state, E u , and the infected steady state, E + i , are locally asymptotically stable within their respective biologically relevant regions, as described by Figure 3. The (light) green region denotes the portion of the (R m , R 0 ) plane in which E u is locally asymptotically stable, while the area that is shaded (dark) red denotes the corresponding local stability region for the E + i state. Interestingly, there is a small overlap of these two regions in which both steady states are locally stable, namely the striped triangular region. An illustration of the change in the regions of Figure 3 generated by differing values of β is provided in Figure 4. Finally, we note that with the original fitted parameter values, the reproduction numbers are R 0 = 1.53 and R m = 0.923, respectively, which corresponds to the development of a persistent viral infection.

4.1.
Basins of Attraction in the Bistable Region. Since both steady states are locally stable in the overlapping region of Figure 3, we expect that differing longterm behavior, and hence different disease outcomes, may arise from variations in initial data. Indeed, this is the case, and we demonstrate this by considering two different simulations of the model with (R m , R 0 ) values within the bistable region. In order to display population values on a biologically pertinent scale, simulation values for T (t) and V (t) are referenced and displayed in the original, dimensional variables, rather than for the dimensionless system.  Table 1.    To demonstrate this further, we simulate the progression of the model for varying (R m , R 0 ) values with other dimensionless parameters held constant and test whether the large-time behavior of these solutions tends towards the E u or E + i steady state. The results of these simulations over the bistable region is shown in Figure 6.
While these simulations provide information regarding the qualitative difference between solutions in this parameter region, they fail to describe how close initial data must be to equilibrium in order to guarantee their stability, i.e. the basin of attraction. To study these basins of the two stable equilibrium states, we perform a number of perturbative computational studies at differing points in the bistable region. The results, shown in Figures 7 and 8, display the sensitivity to initial conditions at each location. Notice that locations closer to R 0 = 1 and further to the right in the bistable region display a greater basin of attraction for E + i than those to the left of this region in the (R m , R 0 ) plane, as displayed by Figures 8(c), (e), and (f). The appearance of a strip of persistent infection steady states increases while moving from left to right within this region. We see that the initial viral load has a minor influence on the shape of the strip for various locations as the width of each strip is not uniform, while the initial T-cell count has a more pronounced affect on the long term behavior. As these initial values are dimensionless, when T 0 is rescaled to represent an actual T-cell count, the basins of attraction for the infected steady state E + i are increased by a factor of nearly 10 3 . Similarly, the basin of attraction corresponding to perturbations in the viral load are increased by around 200 when V 0 is rescaled to represent a true viral load.
An individual infected with HIV typically possesses baseline parameters (see Table 1) corresponding to a location in the (R m , R 0 ) plane lying above R 0 = 1, and hence within the E + i stability region. However, certain parameter values, including the infection rate k and the rate of viral production p, are known to vary widely amongst individuals, and thus the feasible region of attained values within the (R m , R 0 ) plane is quite vast. In particular, k has been reported to be as small as 10 −6 [13] and as large as 10 −2 [12], and this uncertainty could allow the values of R m and R 0 to vary within the range 10 −2 − 10 2 . Thus, we see that the biologically feasible parameter regime extends even into the unshaded region of Figure 3, which leads us to study the dynamics there, as well. Figure 3, we notice that there is a region of the parameter space within which no equilibrium point is locally stable. Thus, one may expect that a different attracting set inherits this property for such parameter values. Indeed, this is the case, and as we will show using R m as a bifurcation parameter, a Hopf bifurcation occurs at the boundary of this domain. In particular, we will take any R 0 > 1, vary R m to move within this region of the (R m , R 0 ) plane, and investigate the stability properties of E + i as they change along the right boundary of the (dark) red region within Figure 3. Since R m will be used to move through the parameter space, we will alter notation when necessary in order to denote certain quantities that depend on this parameter. As demonstrated within  Appendix C, the Jacobian of system (3CM*) evaluated at E + i is

Hopf Bifurcation. From
is the E + i viral population given by Theorem 3.1. Recall that this steady state value satisfies the quadratic equation and thus varies with R m when all other parameters are fixed. The characteristic polynomial associated to ∇f (E + i ) is (8), and in particular the relationship Notice that d 1 , d 2 > 0 since all parameters are positive. Additionally, we establish the following result, which will be useful in our study of the roots of (9) and in the proof of Theorem 4.1. The proof of Lemma 4.2 is contained within Appendix B. Next, denote the corresponding roots of the characteristic polynomial (9) by η i (R m ), i = 1, 2, 3. It follows that the mapping R m → η i (R m ) is smooth, as displayed within Figure 9. To begin our study of the eigenvalues of this system, we first show the existence of a negative real root within the local stability region of E + i . Lemma 4.3. For any R m > 0 within the region of existence for E + i , the characteristic polynomial possesses at least one real, negative root. Additionally, all real roots are negative.
Proof. By the Fundamental Theorem of Algebra, the above characteristic polynomial will have exactly three roots. Thus, it will either possess one real root and two complex roots, or three real roots. Using Lemma 4.2, we find d 0 , d 1 , d 2 > 0 in the region of existence for E + i , and thus any real root must be negative. Hence, the  Table 1.
characteristic polynomial has either three negative real roots or one negative real root and two complex conjugate roots. In either case, the conclusions follow.
Next, we define to be the second Hurwitz determinant of the characteristic polynomial (9). Further, for a fixed value of R 0 > 1, let R * m be the value of R m in the region of existence of E + i such that D 2 (R * m ) = 0. We will see that the curve D 2 (R m ) = 0 corresponds to the right edge of the E + i stability region in Figure 3. The following result shows that as the value of R m is increased beyond R * m , so that the point (R m , R 0 ) lies outside of the E + i stability region, the complex eigenvalues become purely imaginary. Lemma 4.4. At the point R m = R * m , two eigenvalues of ∇f (E + i ), denoted η 2 (R * m ) and η 3 (R * m ), are purely imaginary and conjugate, while the third, η 1 (R * m ), is real and negative.
Proof. At the value R m = R * m , we have D 2 (R * m ) = 0 and thus d 0 (R * m ) = d 1 d 2 . Therefore, (9) can be written as which possess the roots Finally, with this understanding of η 1 , η 2 , and η 3 , we can prove the existence of a Hopf bifurcation across the curve D 2 (R m ) = 0.

STEPHEN PANKAVICH, NATHAN NERI, AND DEBORAH SHUTT
Theorem 4.5. For R 0 > 1, a Hopf bifurcation occurs at the critical value R m = R * m . In particular, as the value of R m crosses R * m , the equilibrium point E + i becomes unstable and a stable limit cycle branches from the equilibrium.
Proof. With Lemmas 4.3 and 4.4, we need only show the transversality condition, dη2 dRm (R * m ) = 0 to prove the existence of a Hopf bifurcation. Using the characteristic polynomial (9) evaluated at η 2 (R m ), we take the R m -derivative to find Evaluating the derivative at R m = R * m and substituting the known value for Thus, multiplying by the conjugate and taking the real part, we find Re .
From this, it follows that Re dη2 dRm (R * m ) = 0 if and only if d ′ 0 (R * m ) = 0. Using (10), we compute this term as Now, using (8) we compute dV + dRm so that and after some algebra and use of (8), this implies Since V + > 0 for R 0 > 1, we conclude that this term is strictly positive. Finally, a brief computation using (10) shows that the remaining term in d ′ 0 (R m ) satisfies Since R 0 > 1 and V + > 0, this term is strictly positive, as well. Hence, we find d ′ 0 (R * m ) = 0 and the proof is complete. To supplement these analytical results, we also include representative simulations of the system for specific values of R 0 > 1 near the point R m = R * m within Figures 10 -12. In particular, Figure 10 demonstrates the stability of the infected steady state within the region in which only E i + is locally asymptotically stable. In such a case (R 0 = 1.25, R m = 6.5 here), the complex eigenvalues of the Jacobian matrix evaluated at E i + possess negative real part, and Figure 10 11, the bifurcation parameter R m is adjusted so that R m ≈ R * m and the location in the parameter plane is at the border of the region of E + i stability and the unshaded region seen within Figure 3. The associated complex eigenvalues now have real parts which approach zero, and the emergence of a periodic orbit in Figure 11(a) becomes more visible. Finally, R m is further increased so that the parameter plane location lies within the unshaded region. The complex eigenvalues now have positive real part, and the solution settles into a periodic orbit as seen within Figure 12 (d) Infected T-cell population Figure 11. Simulation of the dimensionless system for Rm = 9 and R 0 = 1.25. The parameter location is at the border of the (dark) red and unshaded regions of Figure 3. Hence, complex eigenvalues of the system linearized about E + i are of the form ±iβ, signaling a transition in the asymptotic behavior of solutions and the emergence of a stable orbit.
time of transmission has the strongest impact on the long term behavior (as shown by Figure 8), the initial viral load also affects the ability of the virus to establish a persistent infection. Hence, the proposed model can account for differing infection dynamics that were displayed by clinical experiments due to variations in the initial size of the viral inoculum or initial strength of the T-cell count. Additionally, the model yields results that highlight the sensitivity of infection dynamics with respect to variations in initial data. In particular, Figure 8 shows that only certain ranges of the T-cell count promote viral infection, though their width may change with parameter values. This contradicts the general idea that an increase in CD4 T-cells, which serve as target cells for HIV, will necessarily give rise to a greater likelihood of infection, since within the parameter region of bistability, increasing T 0 can actually push the system from the infected state to an uninfected regime. This result is consistent with the clinical immunosuppressant studies of [7,15] which found that in some cases, a depleted T-cell count could give rise to a greater possibility of infection. Therefore, differing initial strengths of the susceptible T-cell population may either promote or inhibit viral replication, and a "sweet spot" in the initial state (seen as   Figure 3). Complex eigenvalues of the system linearized about E + i are of the form α ± iβ with α > 0 yielding instability of the equilibrium. With initial conditions taken near the bifurcation, a stable periodic orbit appears. the dark strips within Figure 8) appears to exist within which infection occurs. In general, the refined model (3CM) further highlights the biological implications of T-cell homeostasis, namely that the propensity of the T-cell population to replenish itself, not merely at a constant rate but due to the appearance of a pathogen, can account for differing infection outcomes.
With the emergence of a Hopf bifurcation in (3CM), there also exists a region of the (R m , R 0 ) plane within which solutions display stable oscillations, thereby mimicking observed biological behavior such as viral blips. Such oscillatory blips are transient spikes in the size of the viral load which often occur in a patient undergoing antiretroviral therapy during the chronic stage of infection [28]. Various studies have found that 20% to 60% of patients with viral suppression experience viral load blips (depending on the antiretroviral regimen used and the frequency of viral load testing), and perhaps one-third of these experience repeated blips. Within the proposed model, blips can occur from a sudden, but small, change in parameter values that increase R m or decrease R 0 from the region in which E + i is stable into the region in which a periodic orbit becomes stable. Such changes in parameter values may realistically arise from interruption or sudden alteration of antiretroviral therapy, which strongly influences the values of k and p in the model, and thus the dimensionless parameters R 0 and R m , as well. Hence, (3CM) may be extended to further explain phenomena during chronic infection when the effects of antiretroviral therapy can drastically alter parameter values in the (R m , R 0 ) plane. Finally, the proposed model can also be generalized to describe the entire time course of infection [12,13] by incorporating additional components and biological effects, though a full dynamical analysis would likely be prohibitive in such a case. and compute the gradient of this function To prove the stability result concerning the uninfected steady state, we evaluate the Jacobian at E u to find Clearly, the eigenvalues are η 1 = −1 and η 2 , η 3 given by the two roots of the quadratic η 2 + (α 1 + α 2 )η + α 1 α 2 (1 − R 0 ) = 0.
By the Routh-Hurwitz criterion, the latter eigenvalues both have negative real part if and only if R 0 < 1. Hence, we conclude by the Hartman-Grobman theorem that E u is locally asymptotically stable if R 0 < 1. Contrastingly, if R 0 > 1 then E u is unstable. Next, we establish the stability properties of the infected steady states. Again, we compute the Jacobian, but evaluate it only using the T-cell steady state value, T = 1, as the value of the viral load differs for E + i and E − i . With this, we find where V ± is given by Theorem 3.1 for either steady state and, in both cases, satisfies the quadratic The associated characteristic polynomial is η 3 + d 2 η 2 + d 1 η + d 0 = 0 where d 0 = α 1 α 2 V ± (1 + βV ± ) 2 −R m + (1 + βV ± ) 2 , d 1 = (α 1 + α 2 )R 0 , The Routh-Hurwitz conditions are clearly met for d 1 and d 2 since all parameter values are positive. We only concern ourselves with the sign of d 0 and showing the other remaining condition, namely D 2 = d 1 d 2 − d 0 > 0. After some algebra and use of (12), we can rewrite d 0 as d 0 = α 1 α 2 1 + βV ± βV 2 ± + R 0 − 1 which means that d 0 > 0 if and only if the condition α 1 α 2 1 + βV ± βV 2 ± + R 0 − 1 > 0 is satisfied. First considering the E + i steady state, we find by Lemma 4.2 that the above inequality holds at every point in the region of existence. To study the condition D 2 > 0, we merely note that this is equivalent to α 1 α 2 1 + βV + βV 2 + + R 0 − 1 < (α 1 + α 2 )R 0 (α 1 + α 2 + R 0 ) and results in (7). Thus, the condition for stability of E + i is complete. Finally, to establish the instability of E − i in every parameter regime that guarantees postivity of this equilibrium, we will show that the conditions which ensure V − > 0 in this case violate the stability criteria. It was previously shown that the conditions b < 0, c > 0, and b 2 − 4ac ≥ 0 are needed in order to arrive at a positive root V − for E − i which satisfies (12), in which case Hence, writing the root of interest, namely and multiplying by the conjugate of the numerator, we find as b < 0, c > 0, and b 2 − 4ac ≥ 0. Rewriting this inequality as bV − + c ≥ −c and substituting parameters for b and c yields Using (12), this inequality is exactly −βV 2 − ≥ −(1 − R 0 ). Rearranging finally yields the condition βV 2 − + R 0 − 1 ≤ 0, which, considering the postivity of V − and parameters, violates the stability criterion α 1 α 2 1 + βV ± βV 2 ± + R 0 − 1 > 0, thereby implying instability of E − i .