Mathematic modeling of viral dynamics in patients infected with HIV-1 has played an important role for understanding the pathogenesis of HIV-1 infection.^{1-6} Most HIV dynamic studies have only considered short-term viral dynamics under ideal assumptions about patient behavior, treatment potency, and drug susceptibility, however. Although these studies characterized HIV replication during a short period of antiretroviral (ARV) treatment under the ideal conditions, the effectiveness of ARV therapies and variations of HIV-1 dynamics during chronic treatment of HIV-1 infection under more practical settings have not been carefully modeled and studied quantitatively. The relation between short-term viral dynamics and longer term antiviral response has been investigated, but conclusions are contradictory.^{6-12} Some studies found no significant correlation between short-term HIV dynamics and long-term antiviral response,^{6,9,11} whereas others reported significant relations.^{7,8,10,12} Thus, more careful studies with appropriate quantitative approaches are needed to clarify these results. In the present study, we introduce an HIV dynamic model to characterize the long-term antiviral response by carefully considering drug potency, drug exposure/adherence, and drug resistance during chronic treatment of HIV-1-infected patients. We applied the proposed models and methods to data from an AIDS clinical study developed by AIDS Clinical Trials Group (ACTG).

#### A5055 CLINICAL STUDY

ACTG Protocol 5055 (A5055) was a phase 1/2, randomized, open-label, 24-week comparative study of the pharmacokinetic (PK), tolerability, safety, and ARV effects of 2 indinavir (IDV)/ritonavir (RTV) regimens plus 2 nucleoside analogue reverse transcriptase inhibitors (NRTIs) in HIV-1-infected subjects failing their first protease inhibitor (PI)-containing regimen consisting of amprenavir (APV), nelfinavir (NFV), saquinavir (SQV), or SQV/NFV combination therapy. All study sites had institutional review board approval to perform the study, and all subjects gave written informed consent before their participation. Subjects were randomized to 1 of 2 IDV/RTV regimens: 800/200 mg twice daily (arm A) or 400/400 mg twice daily (arm B) plus 2 NRTIs. To improve tolerability in arm B, RTV was administered at a rate of 300 mg every 12 hours during the initial 4 days of study treatment and then increased to 400 mg every 12 hours starting on day 5. Eligible subjects included HIV-1-infected individuals â‰¥18 years of age with a plasma HIV RNA level (Roche Amplicor HIV-1 Monitor) â‰¥500 but â‰¤100,000 copies/mL within 45 days of study entry. Subjects were naive to at least 1 NRTI and did not have any active opportunistic infections.

Study visits occurred before entry; at entry (within 14 days of visit before entry); at weeks 1, 2, and 4; and every 4 weeks thereafter through week 24. Plasma HIV RNA testing was done at each study visit. Clinical assessment and laboratory parameters, including CD4 and CD8 cell counts, were performed at all but the week 1 visits. Genotypic drug resistance testing was performed at screening, week 24 (if HIV RNA level was >1000 copies/mL), at the time of virologic failure, and at treatment discontinuation. Phenotypic determination of ARV drug resistance (ViroLogic) was performed at baseline and at the time of virologic failure. An intensive PK evaluation was performed on day 14. Plasma for intensive PK analysis was obtained before the first dose and at 0.5, 1, 2, 3, 4, 5, 6, 8, 10, and 12 hours after an observed IDV/RTV dose. PK parameters of IDV and RTV were determined using noncompartmental methods (WinNonlin, version 3.2; Pharsight Corporation, Mountain View, CA). Calculated PK parameters included area under the curve (AUC_{12h}), maximum plasma concentration (C_{max}), time to C_{max} (T_{max}), oral clearance (CL/F), terminal apparent distribution volume (Vz/F), and elimination half-life (T_{Â½}), where F, CL and Vz are the absolute bioavailability, clearance, and apparent volume of distribution, respectively. Sparse PK samples (trough drug concentrations at 10-14 hours after dosing) were also collected at each visit from weeks 4 through 24. Pill counts were performed at each study visit from weeks 2 to 24 to monitor adherence. A more detailed description of this study can be found in the article by Acosta et al.^{13}

Virologic failure in this study was defined as (1) failure to achieve a confirmed plasma HIV-1 RNA level of <200 copies/mL or a confirmed 1 log_{10} increase from the nadir (lowest recorded plasma HIV-1 RNA level) before week 24, (2) a confirmed plasma HIV-1 RNA level (2 consecutive measurements) of â‰¥200 copies/mL with less than a 1 log_{10} decrease in plasma HIV-1 RNA from baseline by week 8, and (3) a confirmed increase in plasma HIV-1 RNA level to â‰¥200 copies/mL after a prior confirmed level of <200 copies/mL. Otherwise, subjects were defined as having virologic success.

#### MATHEMATIC MODELS

##### HIV Dynamic Model

Mathematic models for HIV dynamics have been developed over the past 1.5 decades. Recent surveys can be found in articles by Perelson and Nelson,^{14} Nowak and May,^{15} Callaway and Perelson,^{16} and Perelson.^{17} We consider a simple HIV dynamic model with antiviral treatment as follows^{18}:

where the 3 differential equations represent 3 compartments: target uninfected cells (*T*), infected cells (*T**), and free virions (*V*). Parameter Î» (day^{âˆ’1} mm^{âˆ’3}) represents the rate at which new T cells are created from sources within the body such as the thymus, *dT* (day^{âˆ’1}) is the death rate of T cells, *k* (day^{âˆ’1} mm^{3}) is the infection rate without treatment, Î´ (day^{âˆ’1}) is the death rate of infected cells, *N* is the number of new virions produced from each infected cell during its lifetime, and *c* (day^{âˆ’1}) is the clearance rate of free virions. The time-varying parameter Î³(*t*) is the antiviral drug efficacy at treatment time *t*. In this model, the difference between the RTI and PI drug actions is not considered but is expected to have only a small effect on long-term HIV dynamics and model predictions. If we assume that the system (Equation 1) is in a steady state before initiating ARV treatment, it is easy to show that the initial conditions for the system are as follows:

As we have shown,^{18} there exists a drug efficacy threshold *ec* = 1 âˆ’ *cdT*/*kN**Î»* such that if Î³(*t*) > *ec* for all *t*, virus is eventually eradicated in theory. If Î³(*t*) < *ec* (treatment is not potent enough) or if the potency falls below *ec* before virus eradication (eg, drug resistance), however, the viral load may rebound (see article by Huang et al^{18} for a detailed discussion). It is important to estimate *ec* for each patient based on clinical data.

##### Drug Efficacy Models

Within the population of HIV virions in a human host, there is likely to be genetic diversity and corresponding diversity in sensitivity to the various ARV agents. In clinical practice, genotypic or phenotypic tests can be performed to determine the sensitivity of HIV-1 to ARV agents before a treatment regimen is selected. Here, we use the phenotypic marker, the median inhibitory concentration (IC_{50}),^{19} to quantify agent-specific drug susceptibility. To model within-host changes over time in IC_{50} attributable to the emergence of new drug-resistant mutations, we use the following function^{18}:

where *I*_{o} and *Ir* are respective values of IC_{50}(*t*) at baseline and time point *tr* at which the resistant mutations dominate. If *Ir* = *I*_{o}, no new drug-resistant mutation is developed during treatment. Although more complicated models for IC_{50} have been proposed based on the frequencies of resistant mutations and cross-resistance patterns,^{20,21} in clinical studies or clinical practice, it is common to collect IC_{50} values only at baseline and failure time as designed in the A5055 study. Thus, this function may serve as a good approximation.

Poor adherence to a treatment regimen is one of the major causes of treatment failure.^{22,23} Patients may occasionally miss doses, may misunderstand prescription instructions, or may miss multiple consecutive doses for various reasons. These deviations from prescribed dosing affect drug exposure in predictable ways. We use the following model to represent adherence:

where 0 â‰¤ *Rd* < 1 (*d* = 1,2), with *Rd* indicating the adherence rate for drug *d* (in our clinical study, we focus on the 2 PI drugs of the prescribed regimen, with *d* = 1 representing RTV and *d* = 2 representing IDV). *Tk* denotes the adherence evaluation time at the *k*th clinical visit.

Highly active antiretroviral therapy (HAART) containing 3 or more NRTIs or nonnucleoside reverse transcriptase inhibitors (NNRTIs) and PIs has proven to be effective at reducing the amount of virus in the blood and tissues of HIV-infected patients. In most viral dynamic studies,^{14,24-26} investigators assumed that the drug efficacy was constant over treatment time. Drug efficacy may actually vary, however, because the concentrations of ARV drugs and other factors (eg, emergence of drug-resistant mutations) vary during treatment.^{14,25,27,28} Also, patients' viral load may rebound because of drug resistance, nonadherence, and other factors.^{29} A simple pharmacodynamic (PD) sigmoidal E_{max} model for the dose-effect relation follows^{30}:

where E_{max} is the maximal effect that can be achieved, *C* is the drug concentration, and *EC*_{50} is the drug concentration that induced an effect equivalent to 50% of the maximal effect. Many different variations of the E_{max} model have been developed by pharmacologists to model PD effects. E_{max} models include the sigmoid E_{max} model, ordinary E_{max} model, and composite E_{max} model.^{30,31} The ordinary E_{max} model describes agonistic and antagonistic (inhibitory) effects of a drug, the sigmoid E_{max} model is more flexible for the steepness or curvature of the response-concentration curve compared with the ordinary E_{max} model, and composite E_{max} models are used for multiple drug effects. More detailed discussions on E_{max} models can be found in the book by Gabrielsson and Weiner^{30} and in the recent article by Huang et al.^{18} Here, we use the following modified E_{max} model to represent the drug efficacy for 2 ARV agents within a class (eg, the 2 PI drugs RTV and IDV):

where

denote the inhibitory quotient (IQ) for drugs 1 and 2 (RTV and IDV in our study), respectively.^{32,33} We refer to

as the adherence-adjusted inhibitory quotient (AIQ);

are the trough levels of drug concentration in plasma (measured after 12 hours from doses taken in our study) and the IC_{50} for the 2 agents, respectively. The parameter Ï† is used to quantify the conversion between in vitro and in vivo IC_{50} that can be estimated from clinical data. We call Ï† the IC_{50} conversion factor. Notice that the larger Ï† is, the lower is the drug efficacy Î³(*t*). Thus, Ï† acts like the IC_{50}. Our model (Equation 6) is a variation of the composite E_{max} models. Note that *C*_{12h}*d* in Equation 6 could be replaced by other PK parameters such as the area under the plasma concentration-time curve (AUC). Note that the drug efficacy Î³(*t*) should include both NRTI and PI drugs in the regimen. We only consider the PI drug effect, however, because the NRTI drug information is not available and the effect of NRTI drugs in this population is likely smaller compared with that of PI drugs.

Lack of adherence, quantified by Equation 4, reduces drug exposure, and thereby reduces drug efficacy, as seen in Equation 6, which, in turn, can affect viral response. Time-varying parameter Î³(*t*) ranges from 0 to 1 and is referred to as the drug efficacy index (the inhibition rate of viral replication) in the viral dynamic model (Equation 1). If Î³(*t*) = 1, the drug is 100% effective, whereas if Î³(*t*) = 0, the drug has no effect. More details regarding viral dynamic models and time-varying drug efficacy models can be found in the articles by Huang et al^{18} and Dixit and Perelson.^{27} Note that if C_{12h}*d*, *Ad*(*t*), and IC_{50}*d*(*t*) are measured from a clinical study and if Ï† can be estimated from clinical data, the time-varying drug efficacy Î³(*t*) can be estimated for the whole period of antiviral treatment.

#### STATISTICAL METHODS

##### Hierarchic Bayesian Modeling Approach

A hierarchic Bayesian model (Bayesian mixed-effects model) allows us to incorporate prior information at the population level into the estimates of dynamic parameters for individual subjects. Denote the number of subjects by *n* and the number of measurements on the *i*th subject by *mi*. For notational convenience, let Î¼ = (ln Ï†, ln *c*, ln Î´, ln Î», ln *dT*, ln *N*, ln *k*)*T*, *Î¸i* = (ln *Ï†i*, ln *ci*, ln Î´_{i}, ln Î»_{i}, ln *dTi*, ln *Ni*, ln *ki*)*T*,Î¸ = {*Î¸i*, *i* = 1, Â·Â·Â·, *n*}, Î¸_{{i}} = {*Î¸i*, *l* â‰ *i*}, and **Y** = {*yij*, *i* = 1,Â·Â·Â·, *n*; *j* = 1,Â·Â·Â·, *mi*}. Let *fij* (*Î¸i*, *tj*) be the numeric solution for the common logarithmic viral load, log_{10} (*V*(*t*)), of the differential Equation 1 for the *i*th subject at time *tj*. The repeated measurements of common logarithmic viral load for each subject, *yij*(*t*), at treatment times *tj*, *j* = 1,2, Â·Â·Â· *mi*, can be written as follows:

where *ei*(*t*) are measurement errors with a mean of 0. The hierarchic Bayesian model can be written in the following 3 stages^{34}:

Stage 1: within-subject variation in common logarithmic viral load measurements:

where

and the bracket notation [A|B] denotes the conditional distribution of A given B.

Stage 2: between-subject variation:

where Î¼ is the population mean of the parameters and *b* is the random effect for subject *i*.

Stage 3: hyperprior distribution:

where the mutually independent Gamma (*Ga*), Normal (*N*), and Wishart (*Wi*) prior distributions are chosen to facilitate computations.^{34-36} The parameters *a, b, Î·,* Î›, Î©, and *v*, which characterize the hyperprior distributions, are assumed to be known from previous studies and the literature.

In the proposed Bayesian approach, we only need to specify the priors at the population level. The population estimates of unknown parameters are easy to obtain from the previous studies or the literature, and the population estimates are usually accurate and reliable. The population prior information also helps to reduce the unidentifiability problem at the level of individuals. The priors for the hyperparameters were obtained from several published studies^{14,15,37-41} and are as follows:

We implemented our Bayesian approach using the Markov chain Monte Carlo (MCMC) procedure, consisting of a series of Gibbs sampling and Metropolis-Hastings algorithms.^{42} The convergence of the MCMC algorithms was carefully monitored.^{35,36,43-45} Sensitivity analyses were conducted for confirming our results.

The hierarchic Bayesian approach allows us not only to borrow information from previous studies but to borrow information across patients in the same study to estimate the kinetic parameters for an individual patient whose data may not be enough to identify all the parameters in the viral dynamic model (see the book by Davidian and Giltinan^{34} for more discussion on this point). Across-patient information is incorporated in stages 1 and 2 via Bayesian theories, and the prior information regarding the estimates of the viral dynamic parameters from previous studies is incorporated in stage 3 in the hierarchic model (Equations 8-10). Thus, the estimates of viral dynamic parameters for an individual patient are based on the data from this particular patient, the data from other patients in the same study, and the prior information from previous studies. The information from these 3 resources is weighted according to the uncertainty of each information component in an efficient and optimal way via Bayesian theories. In contrast, the nonlinear least squares (NLS) method fits a model to the data from an individual patient at a time. The data from an individual patient sometimes may not be enough to identify all the viral dynamic parameters reliably; hence, data from other patients are generally not used. Thus, the key advantage of the hierarchic Bayesian approach, compared with the NLS method, is its efficient utilization of all the information available at hand.

Equation 8 Image Tools |
Equation 9 Image Tools |
Equation 10 Image Tools |

##### Correlation Analysis

Under the assumption for individual Bayesian estimates, we may carry out further analysis for the estimated parameters using standard methods.^{26} We correlated the estimated viral dynamic parameters with baseline host factors, PK parameters, and drug adherence and resistance as well as summary statistics of virologic/immunologic responses. We used the Spearman rank test to evaluate the statistical significance of correlations. Simple linear relations were explored using a robust linear regression method (MM-estimator) because of the presence of outliers (although some of the relations may be nonlinear).^{46} All *P* values are 2-sided with a significance level of 0.05, and no adjustments for multiple testing were made in this exploratory analysis.

#### RESULTS

Forty-four HIV-1-infected subjects were enrolled in this study. Twenty-two subjects were randomized to the twice-daily 800/200-mg IDV/TRV arm (arm A), and 22 subjects were randomized to the twice-daily 400/400-mg IDV/RTV arm (arm B). Thirty-six (82%) of 44 subjects were male. Twenty percent were white, 43% were black, 32% were Hispanic, and 5% were of other races. The median age of subjects was 37 years. Because of toxicities, 1 subject switched treatment from arm A to arm B and 5 subjects switched treatment from arm B to arm A. Note that of the 44 subjects, 42 were included in this analysis; of the remaining 2 subjects, the first was excluded from the analysis because the PK parameters were not obtained and the second was excluded because the phenotype assay could not be completed on this subject. Although drug exposures differed between arms because of the difference in dosing of IDV/RTV, virologic/immunologic responses (HIV-1 RNA copies and CD4^{+} cell counts) and tolerability showed no difference between the 2 treatment arms. More detailed results of pharmacokinetics, ARV effects, and tolerability/safety from this study have been reported elsewhere.^{13}

##### Model Fitting and Parameter Estimation Results

The dynamic model (Equation 1) was fitted to the viral load data from patients in treatment arms A and B separately as well as from the pooled data from the 2 arms using the proposed Bayesian approach. We report the dynamic parameter estimates and their summary statistics for arms A and B in Tables 1 and 2. We observed a large between-subject variation in the estimates of all individual dynamic parameters (the coefficient of variation [CV] ranges from 37% to 380% for different dynamic parameters). For instance, the smallest virus clearance rate (*c*) was estimated as 3.6 day^{âˆ’1} with a corresponding half-life of 0.19 days or 5 hours, and the largest was 21.3 day^{âˆ’1} with a half-life of 0.03 days or 1 hour. The smallest death rate of infected cells (Î´) was estimated as 0.03 day^{âˆ’1} with a corresponding half-life of 23 days, and the largest was estimated as 1.3 day^{âˆ’1} with a half-life of 0.5 days or 12 hours. The smallest number of virions (*N*/2) produced per infected cell was estimated as 12, and the largest number of virions produced per infected cell was 1149. We did not find any significant difference in any of the viral dynamic parameters (Ï†, *Î´*, *dT*, Î», *k*, *N*, *c*, and *ec*) between the 2 treatment arms.

Table 1 Image Tools |
Table 2 Image Tools |

As examples, Figures 1 and 2 present the model-fitting results for 2 subjects from each arm, respectively: 1 subject with virologic success and 1 subject with viral rebound. For patient 11 in arm A (see Fig. 1), no drug resistance developed (IC_{50} is constant during treatment) and adherence to the prescribed regimen was close to perfect, which yielded an estimated drug efficacy above the drug efficacy threshold (*ec* = 0.81) required to suppress the viruses. Thus, this patient had a successful viral load trajectory. For patient 10 in arm A (see Fig. 1), IC_{50} had no significant increase from baseline (IC_{50} = 2.89 and 11.9 ng/mL for IDV and RTV at baseline, respectively) to viral rebound time (IC_{50} = 3.13 and 14.68 ng/mL for IDV and RTV, respectively), and this patient also had a good adherence record (no missing pills were found). The drug efficacy, ranging from 0.57 to 0.60 during the treatment period, was below the drug efficacy threshold of this patient (*ec* = 0.96) required to suppress virus, however, which resulted in a large viral load rebound at week 20. Compared with the average of the patient population, we observed high drug concentrations for these 2 patients (IDV C_{12h} = 422.2 and 2674.5 ng/mL, IDV AUC_{12h} = 69.4 and 72.1 Î¼g*h/mL, RTV C_{12h} = 761.4 and 2515.9 ng/mL, and RTV AUC_{12h} = 46.0 and 37.1 Î¼g*h/mL for the 2 patients, respectively). The parameter Ï† in the drug efficacy model (Equation 6) was estimated with a large difference between the 2 patients, however: Ï† = 29.5 for patient 11 and Ï† = 256.0 for patient 10.

Figure 1 Image Tools |
Figure 2 Image Tools |

For patient 13 in arm B, the viral load was reduced from baseline to below detection (see Fig. 2). No drug resistance was detected. Although there was a minor adherence problem (only a couple of IDV doses were missed), the drug efficacy for this patient was estimated to be above the threshold (*ec* = 0.66) required to suppress the virus. For patient 14 with a viral load rebound in arm B, however, virus with reduced sensitivity to RTV and IDV was detected at failure time (IDV IC_{50} = 11.96 ng/mL and RTV IC_{50} = 55.04 ng/mL at failure time compared with IDV IC_{50} = 4.97 ng/mL and RTV IC_{50} = 22.66 ng/mL at baseline). Although perfect adherence was reported for this patient, the drug efficacy dropped far below the threshold (*ec* = 0.97) required to suppress virus, which resulted in a viral load rebound. Compared with the average of the patient population, the drug concentrations for these 2 patients in arm B were also high (IDV C_{12}*h* = 144.3 and 793.3 ng/mL, IDV AUC_{12}*h* = 12.3 and 20.5 Î¼g*h/mL, RTV C_{12}*h* = 258.4 and 2978.6 ng/mL, and RTV AUC_{12}*h* = 31.2 and 69.5 Î¼g*h/mL for these 2 patients, respectively).

##### Relation Between Baseline Factors and Viral Dynamic Parameters

We have correlated the baseline factors such as baseline viral load (copies/mL), CD4 cells (cells/mm^{3}), age, and weight of patients to the estimated viral dynamic parameters using the Spearman rank correlation test (multiple testing was not adjusted for this exploratory analysis). Baseline viral load and CD4 cell counts were significantly correlated with some of the viral dynamic parameters. Some interesting correlations are plotted in Figure 3, and the rest of them are summarized in Table 3. No significant correlation was observed between the age or weight of patients and viral dynamic parameters.

Table 3 Image Tools |
Figure 3 Image Tools |

##### Correlations of Viral Dynamic Parameters With Virologic and Immunologic Responses

We correlated the estimated viral dynamic parameters with the viral load (log_{10} scale) reductions from baseline to weeks 1, 2, 4, and 24. We found that the death rate of infected cells (Î´) was positively correlated with the viral load decline from baseline to week 4 (*r* = 0.33, *P* = 0.040). The number of virions produced per infected cell (*N*) was positively correlated with the viral load reductions at weeks 2 and 4 (*r* = 0.31 with *P* = 0.048 at week 2 and *r* = 0.44 with *P* = 0.007 at week 4). The death rate of target cells (*dT*) was positively correlated with viral load reduction from baseline to week 4 (*r* = 0.36 with *P* = 0.025). Some of these significant correlations were confirmed by separately fitting the data in arm A or arm B, and similar correlation trends, although not significant, also exist.

The estimated viral dynamic parameters in patients with virologic success were compared with those in patients with virologic failure using the Wilcoxon rank sum test (Fig. 4). First, we included 33 patients with confirmed virologic success or failure status. The patients with virologic success had a significantly lower death rate of infected cells (Î´) and a lower generation rate of new T cells (Î»). There was no significant difference in other viral dynamic parameters (Ï†, *c*, *Î»*, *dT*, *N*, and *k*). Intent-to-treat analysis (treating missing data cases as failure) produced similar results (data not shown).

CD4 T-cell count changes from baseline to week 24 were also correlated with the estimated viral dynamic parameters. Most of the correlations were not significant, except for some weak correlations with marginal significance.

#### DISCUSSION

In the past decade, a great deal of effort has been devoted to studying HIV dynamics.^{1-6,9,10,14-16,39} Most of these studies only modeled viral dynamic data in a short period (2-6 weeks) after initiating an ARV regimen, which was frequently assumed to be 100% effective. In this article, we established the relation of virologic response (viral load trajectory) to drug exposure (pharmacokinetics and adherence) and drug resistance (IC_{50}) via a viral dynamic model. We have proposed using all viral load data obtained during 24 weeks of treatment, incorporating drug exposure and drug susceptibility to estimate viral dynamic parameters and to study long-term HIV dynamics under more practical imperfect treatment settings. We introduced a Bayesian approach to combine prior information with current clinical data for estimating dynamic parameters in the proposed mathematic model for long-term viral dynamics. The proposed model fitted the clinical data reasonably well for most patients in our study, although the fitting for a few patients (<10%) was not completely satisfactory because of unusual viral load fluctuation patterns for these patients.

From the estimated viral dynamic parameters (see Tables 1 and 2), we can see that, on average, the clearance rate of free HIV virions (*c*) was greater than those previously estimated in literature.^{1-4,6,9,10,14,47} The death rate of productively infected cells (Î´) was generally of the same magnitude as those previously estimated in literature, although Î´ tended to be smaller than the best current estimate.^{39} Moreover, the between-subject variations in all viral dynamic parameters were large (with the CV for different dynamic parameters ranging from 37% to 380%). This may indicate a significant heterogeneity in host characteristics and viral species in different hosts, suggesting the importance of individualized treatments. No significant difference was observed for the estimated viral dynamic parameters between arms A and B, a finding that is consistent with the primary analysis results of this study showing no difference in antiviral and other responses to the 2 regimens apart from anticipated differences in PK parameters (attributable to dose differences).^{13}

Attempts to correlate baseline HIV-1 RNA levels and first-phase viral decay rates have yielded conflicting results. For example, Notermans et al^{47} and Wu et al^{10} reported that plasma baseline HIV-1 RNA levels were positively correlated with first-phase viral decay rates. Notermans et al^{47} only observed the positive correlation in arm A with 3-drug HAART administered simultaneously, but no correlation was observed in arm B with RTV monotherapy in the first 21 days (zidovudine [AZT] and lamivudine [3TC] were added on day 21). In contrast, Wu et al^{6,9} found a negative correlation between baseline plasma HIV-1 RNA and first-phase viral decay rates in 2 studies. For these 2 studies, RTV monotherapy was administered to patients in the first 10 days in the first study^{6} and dual drug therapy was used in the other study.^{9} Note that in these previous studies, the first-phase viral decay rate was estimated from a biexponential viral dynamic model assuming perfect treatment. This parameter differs from the death rate of infected cells (Î´), because the decline slope is also proportional to the drug efficacy.^{14,16,48} In the current study, we estimated the death rate of infected cells (Î´) directly by accounting for the drug efficacy. We found a positive correlation between baseline viral load and the death rate of infected cells (Î´) (*r* = 0.38, *P* = 0.0162). In a recent study, Markowitz et al^{39} estimated a large first decay rate as 1.0 day^{âˆ’1} (half-life of 0.7 day) in HIV-1-infected patients treated with a potent HAART regimen, and no correlation of first-phase viral decay rates with baseline viral load or CD4^{+} T-cell counts was observed. The conflicting results of correlation between baseline viral load and the first-phase viral decay rates may need further study.

Some other correlations between baseline viral load and viral dynamic parameters are also interesting. To confirm these correlations, we also used the dynamic parameter estimates from separately fitting each of the treatment arms. Some of these significant correlations were confirmed in both arms and some in only a single arm (see Table 3). A strong negative correlation between baseline viral load and viral clearance rate (*c*) and a positive correlation between baseline viral load and the number (*N*) of virions produced per infected cell were consistent with the fact that the slower clearance rate of virions and a larger number of virions produced from each infected cell result in a higher viral load. The positive correlation between baseline viral load and the generation rate of new target cells (Î») is also interpretable. One possible interpretation is that the higher generation rate of new target cells may lead to more targets for HIV to infect, which may result in a higher baseline viral load level. These significant correlations are consistent with the steady-state condition (Equation 2) of baseline viral load. The positive correlation between baseline viral load and drug efficacy threshold (*ec*) indicates that a patient with a higher baseline viral load may require a more potent regimen to suppress viral replications and may suggest the benefit of initiating ARV therapy with a lower baseline viral load.

Baseline CD4^{+} T-cell counts had opposite relations to most of the dynamic parameters as baseline viral load had. This is presumably attributable to a negative correlation between baseline CD4^{+} T-cell count and viral load (data not shown). In addition, a negative correlation between baseline CD4^{+} T-cell count and the drug efficacy threshold (*ec*) was observed. This indicates that a patient with a lower CD4^{+} T-cell count may require a more potent regimen to suppress viral load, whereas the viral load in a patient with a higher baseline CD4^{+} T-cell count is easier to control. This result supports the concept of initiating ARV therapy at higher baseline CD4^{+} T-cell counts.

We also correlated the viral dynamic parameters to viral load reduction from baseline to weeks 1, 2, 4, and 24. We found that some of the dynamic parameters were significantly correlated with week 2 and 4 (short-term) virologic responses. The death rate of infected cells (Î´), the death rate of target cells (*dT*), and the number (*N*) of virions produced per infected cell were positively correlated with week 2 and/or week 4 viral load reduction (log_{10} scale), although a nonlinear trend of the relation could be observed. These correlations suggest that the higher death rates of infected cells and target cells and the larger number of virions produced per infected cell may result in a larger viral load decline at short term from baseline.

Based on the protocol, the patient's virologic response status was classified as success, failure, or missing for a long-term (24 weeks) response. The patients with virologic success had a significantly lower death rate of infected cells (Î´) and a lower generation rate of new target cells (Î»). These results were not always consistent with the analysis results of short-term response correlations, which may suggest a discrepancy between the short-term viral response and the longer term response. For example, the higher death rate of infected cells (Î´) may result in a larger viral load reduction in the short term (2 to 4 weeks) but may not increase the likelihood of virologic success in the longer term (24 weeks). Although this is perplexing, if immune responses are important in virologic success, perhaps such responses are continued if infected cells live longer.

To examine the effect of dynamic parameter estimates against the prior distributions and initial values, we conducted the sensitivity analyses using the different mean vector Î· of prior distributions and different initial values (data not shown, see the article by Huang and Wu^{42} for details). The sensitivity analysis results can be summarized as follows: (1) the estimated dynamic parameters were not sensitive to both priors and/or the initial values, and final results are reasonable and robust; and (2) when different priors and/or different initials were used, the results follow the same patterns as those presented in this study. The conclusions of our analysis remain unchanged.

In summary, the mechanism-based dynamic model is powerful and efficient to characterize relations of antiviral response to drug exposure and drug susceptibility, although some biologic assumptions are required. Long-term HIV dynamics can be reasonably modeled with careful consideration of the effects of pharmacokinetics, adherence, and drug resistance. Dynamic parameters for individual subjects can be estimated by borrowing information from prior population estimates and across subjects in the same patient population using the novel Bayesian approach. The established models may also be used to simulate antiviral responses of new antiviral agents and new treatment strategies.

Some limitations exist for the proposed modeling method. One problem is that the mathematic model (Equation 1) is a simplified model and there are many possible variations.^{14-16} We did not separately consider the compartments of short-lived and productively infected cells or long-lived and latently infected cells.^{4} Instead, we examined a pooled productively infected cell population. The virus compartment was not further decomposed into infectious virions and noninfectious virions as in the study by Perelson et al.^{3} Thus, different mechanisms of NRTI and PI drug effects were not modeled. In fact, we only considered PI drug effects in the drug efficacy model (Equation 6), because the information on NRTI drugs was not collected in our study and the effect of NRTI drugs was considered less important compared with that of PI drugs. There is a large variation in plasma protein binding for different ARV agents. For example, the IDV is not highly protein bound (60%), but RTV is as high as 98% to 99% protein bound.^{49-51} The pharmacologic activity of ARV drugs depends on unbound drug entering cells that harbor the HIV. Thus, the protein binding difference may need to be considered to refine the drug efficacy model further. We modeled drug resistance using phenotype IC_{50} values instead of modeling genotype viral species separately.^{15} Although a more elaborate model with consideration of more infected cell and virus compartments, more detailed drug effects, and various drug-resistant viral species may provide a more accurate description for underlying long-term HIV dynamics, it may limit the identification of the model parameters because of the complexity of the model; thus, it may limit the usefulness of the more sophisticated models. The trade-off between the complexity and utility of HIV dynamic models should be carefully considered. Further studies on these issues are definitely needed. Another problem is that as measurements of adherence, pill counts may not reflect actual adherence profiles for individual patients. The data quality would affect our estimation results for viral dynamic parameters. More accurate measurements for adherence such as electronic compliance monitoring caps (Medication Event Monitoring System [MEMS]) may improve data quality and model identification. Nevertheless, these limitations would not offset the major findings from our modeling approach, although further improvement may be warranted.

#### ACKNOWLEDGMENTS

*The authors are extremely grateful to the editor and a reviewer for their insightful comments and suggestions, which led to marked improvement of this article. They appreciate the assistance of Ann Walawander, Song Yu, and Jeong-Gun Park for data preparation and computation, and they thank other A5055 and DACS210 study team members and participating investigators: Carl J. Fichtenbaum, Carla Pettinelli, Denise Neath, Elaine Ferguson, Alfred J. Saah, Charles Flexner, Scott Hammer, Michael Saag, and Robert Schooley*.