UNCERTAINTY QUANTIFICATION IN MODELING HIV VIRAL MECHANICS

. We consider an in-host model for HIV-1 infection dynamics developed and validated with patient data in earlier work [7]. We revisit the earlier model in light of progress over the last several years in understanding HIV-1 progression in humans. We then consider statistical models to describe the data and use these with residual plots in generalized least squares problems to develop accurate descriptions of the proper weights for the data. We use recent parameter subset selection techniques [5, 6] to investigate the impact of estimated parameters on the corresponding selection scores. Bootstrapping and asymptotic theory are compared in the context of conﬁdence intervals for the resulting parameter estimates.


1.
Introduction. In 2008, Banks, et al., published a mathematical model [7] that was used to quantify the interaction between the Human Immunodeficiency Virus (HIV-1) and the immune system on an individual patient level. The authors of [7] performed this quantification by estimating parameters for the mathematical model using longitudinal data from a clinical study involving structured treatment interruptions (STIs) in which patients cycled on and off antiretroviral therapy (ART) for some durations of time. The model captured the transient viremia experienced by some patients on therapy with viral load levels suppressed below the detection limit, i.e., the reemergence of virus from latently infected cells even after treatment had successfully reduced the level of free virus. In the present effort we return to the model of [7] for further analysis. A major motivation for revisiting this model is its potential to be readily modified to aid in understanding and predicting outcomes when long term ART patients are removed from treatment.
Individuals infected with (HIV-1) and subsequently treated with ART can have disparate outcomes. A recent study showed that 11 HIV-1 patients maintained low HIV-1 levels without ART after being treated with ART for an average of 36.5 months [26]. However, in the majority of cases, HIV-1 patients must be treated with ART indefinitely, living with toxic side effects [1,27] and with the potential risk of developing Acquired Immunodeficiency Syndrome (AIDS) [23]. Several factors have been elucidated that can help explain how HIV-1 escapes the action of ART and leads to disparities in ART outcomes. One such factor that has been of particular interest recently is the ability of HIV-1 to lay dormant in so called "reservoir cells". In this scenario, a subset of infected cells harbor the virus in a latent state in which the virus is not produced at all or at undetected levels [3,12]. Since ART targets the various stages of viral replication, HIV-1 that is integrated into the DNA of these reservoir cells can evade ART and reseed the HIV-1 infection when patients are taken off ART [7,29,37].
There is evidence that the presence and size of cellular HIV-1 reservoirs are associated with the ability of the patient to control the HIV-1 infection, as well as external factors such as the timing and duration of ART. For example, some patients called "elite controllers" possess the innate ability to suppress the HIV-1 infection without ART. In these patients, the HIV-1 reservoir is significantly lower than in patients that require ART [17,20,25]. Among HIV-1 patients that require ART, those treated during the earliest stages of infection have a smaller reservoir than those treated at later stages [4]. Furthermore, some patients that are able to suppress HIV-1 after the termination of long term ART therapy have low HIV-1 reservoirs in central-memory CD4+ T cells [26]. These clinical findings emphasize the need to quantify the interaction between HIV-1 replication, the host immune system, latent reservoirs, and disease stage in a patient-specific manner. Such quantification could enable more accurate predictions of how patients will respond to removal from ART and lead to improved personalized treatment regimens.
While Banks et al., found that some of their model parameters could be taken as constant across patients, other parameters, including some describing the latent HIV-1 reservoir, varied between patients. Since these parameters were used to predict different HIV-1 patient responses to STIs, it is important to quantify the confidence in their estimation. The assessment of parameter uncertainty is a highly non-linear problem that depends on the mathematical model, the statistical error model for parameter estimation, and the data set. Although parameter uncertainty quantification is widely held as an important step in the modeling process, this is usually done only after an iterative modeling process (as described in [11]) in which conceptual as well as mathematical models have been developed and validated with data. The first stages are often focused on the predictive capability and biological accuracy of the model. Since parameter uncertainty was not assessed for the Banks et al HIV-1 model, we perform this assessment here and discuss the resulting biological implications of our findings. This will be of particular importance to our future efforts using modifications of the model of [7] as a foundation for inclusion of reservoir cell dynamics in long term predictive models.
We examine the uncertainty in parameter estimates for data from several patients that underwent one or more ART interruptions. But first, because the field of HIV-1 biology has developed rapidly since the model of [7] was published, we provide an in depth model description to verify that the Banks et al., HIV-1 model agrees with current biological understanding of HIV progression. We proceed by describing the data and inverse problem methodology. We then analyze the statistical error model and use these results to perform parameter estimation and quantify within-host parameter uncertainty with asymptotic distributional analysis. The results with asymptotic analysis reveal potentially unsettling aspects in regard to certain parameter estimates in that relative large standard errors with corresponding large confidence intervals are obtained. We then turn to a recently developed methodology, parameter subset selection [5,6], in attempts to interpret and understand these findings. Finally we use bootstrapping techniques to refine and improve our uncertainty quantification for the model of [7].
2. Model description. A system of ordinary differential equations (ODE) has been developed and validated with clinical data in [7] to describe the pathogenesis of HIV infection, including the role of CD4+ memory cells that serve as a major reservoir of latently infected cells. This model does not reflect the nature of all immune responses and all host and viral dynamics but instead captures the most prominent biological features thought to be represented in the data. It has also been shown to accurately recapitulate the within host dynamics of CD4+ T cells and virus levels from clinical patient data [7] and to be an effective predictive model.
The biological model is depicted in Figure 1, where T 1 represents uninfected activated (antigen-specific) CD4+ T cells and T * 1 is infected activated CD4+ T cells. Uninfected resting, or not activated, CD4+ T cells are represented by T 2 and infected resting CD4+ T cells are given by T * 2 . Infectious free virus (virus that is capable of infecting other cells in the plasma) is represented by V I and non-infectious free virus (virus that is yielded inactive by protease inhibitors) is represented by V N I . HIV-specific effector CD8+ T cells are given by E 1 and HIV-specific memory CD8+ T cells are represented by E 2 . Compartmental diagram for the mathematical model given in [7]. Solid black arrows indicate death/clearance, solid grey arrows indicate birth/input. PI and RTI denote interactions affected by protease and reverse transcriptase inhibitors, respectively.
The corresponding compartmental ODE model for in-host HIV infection dynamics developed in [7], based on balance laws, is given by HIV-specific memory CD8+ T cells Table 1. Explanation of the states used in the model (1). equations: with initial condition vector [ . A number of terms in this model are based on the law of mass action, where the rate of change in the size of a population is directly proportional to the size of the population at any given time, so the population changes with an exponential rate. Examples of these terms are −d 1 T 1 and γ T T 1 in theṪ 1 equation. Other terms are based on Michaelis-Menten kinetics, where the rate saturates (approaches a maximum). An example of this type of term is a T V I V I +K V in theṪ 1 equation. The functions ξ 1 (t) = 1 u(t) and ξ 2 (t) = 2 u(t) represent the impact of the treatment, where 1 is the effectiveness of the reverse transcriptase inhibitor, 2 is the effectiveness of the parameter explanation δ viral produced lysis rate of T * 1 d 2 T 2 and T * 2 natural death rate δ E2 death rate of E 2 m rate of removal by cell lysis of T * 1 from the system by E 1 γ T rate at which T 1 and T * 1 differentiate into T 2 and T * 2 , respectively c natural clearance rate of V I and V N I δ E1 constant death rate of E 1 λ E source term for E 1 k 2 production rate of T * 2 due to encounters between T 2 and V I that is less than k 1 ρ 1 rate of removal of V I through successful infection of T 1 ρ 2 rate of removal of V I through successful infection of T 2 d 1 natural death rate of T 1 2 relative effectiveness of protease inhibitor (PI) a A activation rate of T 2 and T * 2 by non-HIV antigen 1 relative effectiveness of reverse transcriptase inhibitor (RTI) p T net proliferation of T 1 and T * 1 due to clonal expansion and programmed contraction p E net proliferation of E 1 due to clonal expansion and programmed contraction k 1 production rate of T * 1 from encounters between T 1 and V 1 N T number of RNA copies produced during the process of T * 1 lysis Table 2. Explanation of the parameters used in the model (1).
protease inhibitor, and u(t) is the HAART drug level, which is 0 when the patient is off treatment and 1 when the patient is on treatment. The efficacy of treatment is represented by f , where 0 ≤ f ≤ 1 in theṪ * 2 compartment. The factors 10 3 are included since CD4+ and CD8+ T cells are measured as cells/µL of blood and virus is measured as RNA copies/mL of plasma. Individual parameters are defined in Table 2. For a more in-depth explanation of the model, refer to [7].
Initially, we looked into the accuracy of this model, since [7] was published in 2008 and it is possible that a number of advancements in the understanding of HIV could have occurred in the intervening five years. We attempted to use current knowledge to confirm the accuracy of this model. We investigated six assumptions: (1) whether the proliferation of CD4+ and the activation of CD8+ memory T cells was correctly represented, (2) whether the role CD4+ cell-produced interferons play in CD8+ cell differentiation was properly included, (3) whether it was accurate to exclude a death rate for infected activated CD4+ cells, (4) whether it is correct to assume that one free virus particle is responsible for each new infection, (5) whether the dynamic nature between CD4+ cells and immune effector cells is accurate, and (6) whether the assumptions made about the activation of resting CD4+ T cells are consistent with recent research understanding.
First we confirmed the accuracy of the proliferation of CD4+ and activation of CD8+ memory T cells into effector T cells in our model. A paper from 2011 [15] outlines the proliferation of CD4+ T cells and activation of CD8+ T cells during an HIV infection. Naive CD4+ T cells are recruited in response to CD4+ T cell depletion whereas CD8+ memory T cell activation is driven by HIV RNA levels [15].
We confirmed that these were accurately incorporated into the model. The density dependent regulation of CD4+ cells is incorporated into the model by the terms λ T Ks V I +Ks − d 2 T 2 that describe how the rate of CD4+ production increases as CD4+ levels decrease. Viral load level is also incorporated into these terms and reflect the diminished thymic production of CD4+ cells when viral load is too high [7]. The HIV RNA dependent activation of CD8+ memory cells into CD8+ effector cells is incorporated into the model by the term p E a E V I V I +K V E 2 , where virus level promotes the activation of CD8+ memory cells (E 2 ), adding to the pool of CD8+ effector cells (E 1 ). The parameter p E represents the net proliferation of CD8+ effector cells once activated. Therefore the model is a reasonable approximation to these dynamics.
Second, we investigated the role of interferons produced by CD4+ T cells to confirm their role in the model. Recent research has emerged about the importance of interferons in the relationship between CD4+ and CD8+ cells [15,28,38]. In our model, a dotted line titled CD4+ help connects non-infected activated CD4+ cells (T 1 ) and infected activated CD4+ cells to the differentiation of immune effector cells (E 1 ) into immune memory cells (E 2 ): Figure 1. This is consistent with the role of Type 1 interferons in the response of CD8+ cells to viral infection [15,38]. Type 1 interferons are cytokines that are induced directly by virus infection. IFN directly promotes the proliferation of antigen-specific CD8+ T cells [15,38]. This mechanism is also reasonably incorporated into the model. Third, the model includes a natural death rate for uninfected activated (T 1 ), uninfected resting (T 2 ), and infected resting cells (T * 2 ), but not for infected activated cells (T * 1 ). We investigated the current understandings of these assumptions. There is research that suggests active infected cells do not experience apoptosis. One study [24] analyzed lymph node sections of four perinatally HIV-infected children and identified productively infected cells as well as apoptotic nuclei. They found that though the two identified cell types were often close to one another, they were never overlapping. They took similar samples from an additional child and four adults and observed overlap but only rarely. Also, newly HIV-infected cells produce viral accessory proteins, which likely assist them in avoiding apoptosis [16]. These proteins include Negative Regulatory Factor (Nef), Viral Infector Factor (Vif), and Viral Protein R (Vpr). For example, Nef hinders p53-dependent UVinduced apoptosis, which otherwise occurs in response to viral attachment and entry involving HIV Gp120. Therefore, it is reasonable to exclude a natural death term for infected activated cells (T * 1 ). Fourth, in [7], the simplifying assumption is made that one free virus particle is responsible for each new infection. The free virus compartment V I includes the term 10 3 [(1−ξ 1 (t))ρ 1 k 1 T 1 +(f ξ 1 (t))ρ 2 k 2 T 2 ]V I representing the removal of free virus through the infection of T 1 and T 2 . The parameters ρ 1 and ρ 2 , which represent the number of virus particles lost due to a single successful infection are both set equal to 1 (copies mL-blood/cells mL-plasma). In [7] it is stated that ρ i can be changed to be greater than 1 to account for multiple virus particles being responsible for each new infection. We investigated whether our model should include ρ i = 1 or ρ i > 1. With reference to a paper from 2010 [28] and a paper from 2013 [16], we decided to leave ρ i = 1.
Fifth, we investigated the terms describing the production and removal of CD8+ immune effector T cell (E 1 ). In the E 1 compartment the term λ E represents the constant differentiation of naive CD8+ T cells (assumed readily available) into effector cells (E 1 ) and the term δ E1 E 1 represents exponential natural death of E 1 . The combination of the birth term have on the effector CD8+ cell compartment (Ė 1 ). This interaction is dynamic because the presence of T * 1 both stimulates the proliferation of E 1 and causes the aging?of The presence of T * 1 stimulates the proliferation of E 1 through antibody presentation [39]. At the same time, infected activated CD4+ cells (T * 1 ) contribute to the death of the effector CD8+ cells (E 1 ) due to the mechanism of activation-induced cell death (AICD). CD8+ cells have Fas ligand proteins on their surface, which bind to the Fas proteins of target cells and induce apoptosis in the target cell. However, CD8+ cells also have Fas proteins on their surfaces and while virgin CD8+ T cells are insensitive to the ligation of their own Fas proteins, as CD8+ T cells are repeatedly stimulated by infected active CD4+ cells, they become increasingly sensitive and susceptible to Fas-mediated killing [35]. Therefore the net effect of these terms onĖ 1 is a function of T * 1 . This relationship can be seen in Figure 2.
representing early infection when antigen presentation by emerging T * 1 stimulates the production of effector cells [39].
representing the scenario when death due to repetitive T * 1 -induced stimulation overwhelms T * 1 stimulated birth. Therefore CD8+ immune effector T cell (E 1 ) generation is depicted in this model in a biologically feasible way.
Finally, it was noted that a large section of [7] pertaining to the latent reservoir CD4+ T cell (T * 2 ) activation into infected activated CD4+ T cells (T * 1 ) is based on a paper written in 2002 and makes the assumption that activation of infected HIV-specific resting CD4+ T cells (T * 2 ) depends on the virus concentration with a half saturation constant K V and maximum activation rate a T . These assumptions are also made for activation of uninfected HIV-specific resting CD4+ T cells (T 2 ) into uninfected activated CD4+ T cells (T 1 ). We investigated the accuracy of these assumptions based on recent literature. Recent research [32,34] supports the assumption dendritic cell (DC) mediated activation of HIV latent reservoir cells. Dendritic cells are antigen-presenting cells that activate the adaptive immune system and can stimulate latent cells to become active. Therefore, the amount of virus can directly affect the amount of antigen-presenting cells that activate HIVspecific resting CD4+ T cells (T * 2 ).

Data and methods.
3.1. Data sets. Data used by the authors of [7] include patients from a clinical trial at Massachusetts General Hospital that underwent anti-retroviral therapy (ART) and had at least one ART interruption. The collected data for each patient are total CD4+ T cell count/µL-blood and total RNA copies/mL-plasma. If there are less than 400 copies/mL-plasma for a standard assay or 50 copies/mL-plasma for an ultra-sensitive assay, then the viral load value is censored to be the limit of quantification. For these data points, the expectation maximization (EM) algorithm was applied in [7]. For this work, we chose to model the data for Patients 1, 27, and 33 from [7] as representatives of patients who underwent three, two, or one structured treatment interruptions (STI), respectively. The EM data points for each patient were combined with the uncensored data to form a data set that could 1 +K d , and the difference between the birth and death terms against T * 1 be used for parameter estimation and uncertainty analysis. In total, there were 102 data points for viral load and 84 data points for the CD4+ T cell count over 1527 days for Patient 1, 36 viral load data points and 35 CD4+ T cell counts over 1379 days for Patient 27, and 75 viral load data points and 52 CD4+ T cell counts over 1302 days for Patient 33.

3.2.
Inverse problem methodology. In order to estimate a vector of parameters q, we use the data for CD4+ T cell count, {y i 1 } N1 i=1 , and data for the viral load, {y j 2 } N2 j=1 . The mathematical model divides total CD4+ T cells into uninfected activated (T 1 ), infected activated (T * 1 ), uninfected resting (T 2 ), and infected resting (T * 2 ), so the total CD4+ T cell count is represented in the model by Similarly, since viral load RNA copies are separated into infectious free virus (V I ) and non-infectious free virus (V N I ), total RNA copies is represented in the model by z 2 = V I + V N I . We used a vector observation version of iterative weighted or generalized least squares framework with the following statistical error models for the CD4+ and viral RNA observables, with the corresponding mathematical models z 1 and z 2 , respectively: where .., N 2 , and q 0 are the usual hypothesized "true" parameter values [9,11]. Hereγ = (γ 1 , γ 2 ) are weighting factors that were determined using analysis of residuals [9,11]. In principle it is possible to use a multistage estimation procedure [19] to determine the values of the statistical model parameters γ 1 , γ 2 , but we did not do that here. Instead we followed standard statistical methodology [9,13,14,19,31] of evaluating residuals for randomness patterns for numerous values of γ 1 and γ 2 . As reported in Section 4.1 below we found γ 1 = 0 and γ 2 = 1.2 were appropriate.
The iterative weighted or generalized least squares estimator [9,11,13,18,19, 31] for n = (N 1 , N 2 ) observations is given by (hereafter we take γ 1 = 0, with the variance components given by (this follows from (2) and (3)) We note that equations (4)-(6) involve the unknown (and to be estimated) parameters q 0 and hence one needs to employ a generalized least squares algorithm as explained below to carry out the desired minimization. We also discuss below the need to weight the error for the viral load in proportion to the model solution as well as to determine the best choice of the value for the exponent γ = γ 2 . All inverse problems involving the 15 parameters q andσ 2 1 ,σ 2 2 were performed in Matlab using the f minsearch function and a vector observation version of the iterative generalized least squares estimation procedure as described in Section 3.2.9 of [11], Sections 3.2.5-3.2.6 of [9] and Algorithm 3.1 below. In this algorithm ε is a user defined threshold tolerance used to set a termination criterion and "./ denotes element-wise division.
Let g 1 = T 1 + T * 1 + T 2 + T * 2 and g 2 = V I + V N I . The sensitivities are computed by solving the system of equations d dt where x and q are the state variables and the parameters being estimated, respectively. We define the 2 × p matrices for i = 1, ..., N 1 , and for j = 1, ...., N 2 We also define the 2 × 2 matrices Then the p × p matrices D iT and ∂z 2 ∂q l (t j 2 ; q 0 ), k, l = 1, ..., p, j = 1, ..., N 2 , respectively. We then define the p×p Fisher matrix F (q 0 ) = (F k,l (q 0 )) with elements The approximate Fisher matrix (7), obtained by evaluating (12) atq n ≈ q 0 , is then used to compute the standard errors for each element (k = 1, 2, . . . , p) ofq n , given by The 100(1-ρ)% confidence intervals can be computed based on the confidence level parameters associated with the generalized least squares estimators q n = q n GLS : where ρ is chosen to be small (e.g., ρ=0.05 for 95% confidence intervals) and in the examples here). The corresponding 95% confidence intervals are then given by with ρ = 0.05.

3.4.
Numerical solution of the model and sensitivity equations. The ordinary differential equation model (1) and the sensitivity equations (8) were solved on a log10 scale with the Matlab solver ode15s. As noted in [7], the system (1) is more efficiently solved as a log-transformed system. Consequently, the sensitivity equations (8) are also more efficiently solved when log transformed because their solution depends on the solution of the original system (1). Solutions to (1) were mapped back to the non-log transformed states prior to computing the least squares minimizer in (4). Hence, solving the log-transformed version of (1) does not affect the iterative generalized least squares Algorithm 3.1 or the statistical error model

Parameter estimation and evaluation of the statistical error model.
Prior to uncertainty quantification, we investigated whether the number of structured treatment interruptions can affect parameter estimation. To accomplish this, we calculated standard errors and the corresponding 95% confidence intervals for patients with one, two, and three structured treatment interruptions (STIs). We analyzed data from Patient 1 (3 STIs), Patient 27 (2 STIs), and Patient 33 (1 STI) from [7]. We estimated patient specific parameters using a generalized least squares framework (see Eq. (2) -(6)). We estimated a subset of parameters in the mathematical model, since Banks, et al., found that only 16 out of 39 total parameters in the model had significant variation between patients. We also fixed V I (0) to the initial viral load data point under the assumption that the amount of non-infectious virus is negligible at the time of treatment initiation. We argue that this assumption is reasonable because the primary cause of non-infectious production is ART (see the term ξ 2 (t) in (1)). Thus, we estimated the subset of 15 parameters {λ T , We variedγ = (γ 1 , γ 2 ) in order to test the assumption that the normalized residuals e i 1 , e j 2 in equations (2) and (3) are i.i.d.. The validity of any uncertainty analysis depends on the correctness of these underlying statistical modeling assumptions. We found that the model accurately fit the virus data for Patients 1 ( Figure  3), 27 ( Figure 5), and 33 ( Figure 7) when using γ 1 = 0 and for several different values of γ = γ 2 between 0 and 2, and that the corresponding residuals [11] were found to be approximately randomly distributed when γ 2 = 1.2 (Figures 4, 6, and  8). Thus, we set γ = γ 2 = 1.2 in the uncertainty analyses we present below. We note that in the figures below we plotted and analyzed the normalized residuals for viral load

4.2.
Uncertainty analysis using asymptotic theory. We next used the asymptotic theory techniques described in Section 3.3 to compute uncertainty for parameters in our model of HIV infection and immune system dynamics (1). The resulting parameter estimates and standard errors for Patients 1, 27, and 33 are given in Table 3. We found that the standard errors for several parameters in each patient were larger than the parameter estimates themselves. Together with our suspicion that there are likely several unidentifiable (or correlated in the statistical sense) parameters in our model (1), these results suggest that the application of asymptotic theory for parameter uncertainty quantification will yield unreasonably large standard errors for the set of 15 parameters {λ T , d 1 , ε 1 , k 1 , a T , ε 2 , N T , b E2 , a E , p E , a A , p T , T 1 (0), T * 1 (0), T 2 (0)} regardless of the number of structured treatment interruptions. We next investigated whether there was a subset of parameters in our model which we could estimate and expect reasonable standard errors, i.e.,  such that the 95% confidence interval is narrow relative to the size of the parameter estimate itself.

Parameter subset selection.
To better understand how individual parameters contributed to the overall high asymptotic standard errors we found for our model, we investigated whether there were subsets of the 15 parameter set P S = {λ T , d 1 , ε 1 , k 1 , a T , ε 2 , N T , b E2 , a E , p E , a A , p T , T 1 (0), T * 1 (0), T 2 (0)} for which the application of asymptotic theory would result in reasonable standard errors. We used an algorithm recently developed in [5,6] that selects a parameter subset based  on minimizing the sum of normalized standard errors (selection score) for a given number of parameters n p as described below.
Given a set of size p of parameters q, and a number n p ≤ p, the subset selection algorithm finds a subset of parameters of size n p that minimizes a selection score as described in [6]. To implement this procedure one first needs a set of parameter estimates {q 1 , ...,q p }, with corresponding standard errors {SE 1 , ..., SE p }. One then introduces the coefficients of variation forq i   Table 3. Parameter estimates and standard errors for Patients 1, 27, and 33 computed using asymptotic theory. and take the selection score given by the Euclidean norm in R np of ν(q) = (ν(q 1 ), . . . , ν(q np )) T , i.e., as an uncertainty quantification for the estimatesq. The components of the vector ν(q) are the ratios of each standard error for a parameter to the corresponding nominal parameter value. These ratios are dimensionless numbers warranting comparison even when parameters have considerably different scales and units (e.g., in case of the HIV-model N T is on the order of 100, while k 1 is on the order of 10 −5 ). A selection score α(q) near zero indicates lower uncertainty possibilities in the estimation, while large values of α(q) suggest that one could expect to find substantial uncertainty in at least some of the components of the estimates in any parameter estimation attempt. For a given number of parameters, n p , the algorithm results in a parameter subset of size n p with a minimal selection score α min (n p ). Overall, we found that the selection score patterns for Patients 1, 27, and 33 were similar. Not surprisingly, the selection score for each patient shows the same monotone increasing trend, even in presence of the different numbers of structured treatment interruptions (STIs) between patients (Figure 9).  Figure 9. The minimum selection scores for Patients 1, 27, and 33 plotted on a log 10 scale. To further dissect how the individual parameters contributed to the increasing trend in the selection score, we analyzed the order in which the parameters were selected as the number (n p ) of parameters in the subset increased. Interestingly, we found that for each patient, the parameter subset for n p was contained in the parameter subset for n p +1 for almost all subsets, with only a few exceptions (Table  4). We observed that, in addition to the selected parameter subsets being ordered for each patient, the subsets themselves were strikingly similar between patients for any given n p (Table 5).
We next assessed whether there was a parameter subset for each patient for which we could expect reasonable standard errors using asymptotic theory. We set a threshold on the normalized standard error for each parameter in a given parameter subset in order to ensure that every 95% confidence interval was narrower than the parameter estimate itself. To accomplish this, we used a threshold of 1 1.96 = 0.5102 (see (15)). Using this criteria, the same seven parameters were found for Patient 1 and Patient 27 {N T , 2 , λ T , p T , T 2 (0), p E , T 1 (0)} (Figures 10, 11), and three parameters were found for Patient 33 {λ T , N T , T 2 (0)} ( Figure 12). Similar to the parameter nesting trend we found for the overall parameter subset selection described above, the set of three parameters found for Patient 33 is contained in the set of seven parameters found for Patients 1 and 27. These results suggest that   Table 4. The chosen parameters for each subset size n p = 1, ..., 15 for Patients 1, 27, and 33. Entries in the table represent whether the parameter in the column header was chosen for a patient at subset size n p . the number of STIs affects the ability of asymptotic theory to calculate uncertainty for parameter estimates for our model of HIV dynamics, with more STIs leading to a larger parameter subset with reasonable standard errors.
We note that the parameter subset selection algorithm does not re-estimate parameters when minimizing the selection score, which is composed of parameter standard errors. Since the standard errors themselves rely on the parameter estimates, we tested whether re-estimating only the chosen parameter subsets based on  Table 5. The number of elements in the intersection of parameter subsets selected for each patient (Patients 1, 27, or 33) for a given number of parameters n p . These data reflect the strong agreement between the chosen parameter subsets for each patient, as all entries in each row are greater than n p − 2.
our threshold criteria could change the asymptotic standard errors. We found that both the parameter estimates and the standard errors were very similar before and after the re-estimation process (Tables 6, 7 Table 6. Parameter estimates, standard errors (SE), and 95% confidence intervals for the seven parameters selected by the subset selection algorithm for Patient 1. The results are shown for parameters estimated before the subset selection algorithm was run (B.S.) and after the re-estimating those same parameters after they were selected by the subset selection algorithm (A.S).

Uncertainty analysis using bootstrapping.
We also estimated parameter uncertainty for the model (1) using bootstrapping. We used a non-parametric form of bootstrapping as given in [18,Chapter 11], [13, p.27-28], [14,21,22]. We note that asymptotic theory is computationally preferred over bootstrapping since asymptotic theory took approximately 15 minutes per patient whereas bootstrapping took approximately 350 computing hours per patient. We again investigated whether the number of structured treatment interruptions can affect uncertainty quantification by performing bootstrapping for Patients 1, 27, and 33.   Table 7. Parameter estimates, standard errors (SE), and 95% confidence intervals for the seven parameters selected by the subset selection algorithm for Patient 27. The results are shown for parameters estimated before the subset selection algorithm was run (B.S.) and after the re-estimating those same parameters after they were selected by the subset selection algorithm (A.S).
To implement bootstrapping to quantify the uncertainty in parameter estimates for equations (1), we first obtained a vector of estimated parameters,q, using the above inverse problem methodology. We then calculated the standardized residuals,   Table 8. Parameter estimates, standard errors (SE), and 95% confidence intervals for the three parameters selected by the subset selection algorithm for Patient 33. The results are shown for parameters estimated before the subset selection algorithm was run (B.S.) and after the re-estimating those same parameters after they were selected by the subset selection algorithm (A.S).
r i j for these estimates for each observable (j = 1, 2): , j = 1, 2, . . . , N 2 . Bootstrap sample points were then created by sampling residuals with replacement from the set {r i j } Nj i=1 for the observable z j , j = 1, 2 and adding them to the model solutions z 1 , z 2 evaluated atq to obtain bootstrapping sample points y i 1 = z 1 (t i 1 ;q)+ r i 1 , i = 1, 2, . . . , N 1 , y j 2 = z 2 (t j 2 ;q) + z γ 2 (t j 2 ;q)r j 2 , j = 1, 2, . . . , N 2 , (see [9, p. 94-96]). We repeated this process M = 1000 times to create 1000 data sets, each with the same number of data points for the CD4+ T cell count and viral load as in the original data set, respectively. We then conducted an inverse problem on each of these 1000 data sets and stored the parameter estimatesq m (m=1,..., 1000) in a matrix, Q BOOT . With these values, the mean, variance, and standard errors for the parameters can be calculated using the following formulas given in [8,9,13,14,18]: For each parameterq k (k = 1, ..., dim(q)) in our original parameter estimate vectorq, the 95% confidence interval was calculated directly as the range between the 25-th and 975-th entries in the ordered set of 1000 parameter estimates from bootstrapping. (We note that there are numerous ways [21] in which one might compute confidence intervals using the empirical bootstrapping distribution. We use the simple percentile method here.) The bootstrapping for Patient 1 resulted in normal distributions for most parameter estimates ( Figure 13) with narrow 95% confidence intervals ( Table 9). The distributions for the parameters a E and a T were positively skewed for Patient 1. The bootstrapping for Patient 27 resulted in normal distributions for all parameters except for a A , p T , T * 1 (0), and ε 1 which showed positive skewness, and ε 2 which showed negative skewness ( Figure 14). Skewness was most evident for the parameter ε 2 , which was estimated to be 0.98237. We note that this parameter is bounded between 0 and 1 by biological definition, so it is possible that the skewness was a result of performing the inverse problem in a bounded parameter domain. The 95% confidence intervals for Patient 27 are given in Table 10. One possible source of skewness for the other parameter bootstrapping distributions for Patient 27 is the existence of correlated parameters (recall that the bootstrapping approach tacitly treats the parameters as random variables for which the concept of correlation is meaningful). The bootstrapping for Patient 33 resulted in narrow but skewed distributions in most parameters ( Figure 15). This indicates that although the bootstrapping procedure resulted in narrow 95% confidence intervals for Patient 33 (Table 11), there may exist even more parameter correlations than for Patient 27. Taken together, these findings suggest that standard errors, as calculated from bootstrapping, do not appear to correlate solely with the number of STIs observed in a given patient.  Table 9. Standard errors and 95% confidence intervals from bootstrapping for the Patient 1 parameter estimates. These values correspond to the distributions in Figure 13.   5. Discussion. In anticipation of continuing efforts with additional data sets involving reservoir cells, we have returned to further uncertainty quantification investigations of a model [7] previously developed by members of our group. In a generalized least squares approach (in contrast to our ordinary and weighted least squares efforts in [2,7]), we first developed an acceptable statistical model (2)- (3) in which to pursue uncertainty quantification via both asymptotic and bootstrapping considerations. The resulting value of γ 2 = 1.2 for the observed viral load in place of the more usual absolute (γ = 0) or relative (γ = 1) error models is not surprising in light of other in-depth studies [10] with complex patient-dependent data sets. Our subsequent findings that the number of model parameters that have reasonable asymptotic standard errors increases as the patient data set becomes more dynamic, (i.e., as the number of treatment interruptions increases) supports our mathematical intuition. Ironically, this implies that we may be able to better predict outcomes in ART cessation in patients who go off ART either unpredictably, (i.e., individuals who don't adhere to their scheduled treatments), or through some type of planned STIs. Of course, this implication relies on the assumption of being able to continue to reliably collect viral load and CD4+ data from these patients despite their interruption to ART scheduling. Moreover, we expect this number of reliably-estimated parameters to be bounded (in our example possibly < 15)   Table 11. Standard errors and 95% confidence intervals from bootstrapping for the Patient 33 parameter estimates. These values correspond to the distributions in Figure 15.
depending on the basic identifiability of individual parameters even if one could produce data with an unlimited number of treatment interruptions. Interestingly, we also found that the minimum score parameter subset selected for the patient with one STI was contained in the subset for the patients with two and three STIs. This finding strongly suggests that there is an order in which parameters can be reliably estimated between patients, with more dynamic data leading to a larger number of parameters. Of note, these subsets appear to contain parameters describing viral reproduction and ART efficacy against viral reproduction (N T and