Share this article on:

Modeling Long-Term HIV Dynamics and Antiretroviral Response: Effects of Drug Potency, Pharmacokinetics, Adherence, and Drug Resistance

Wu, Hulin PhD*; Huang, Yangxin PhD*; Acosta, Edward P PharmD†; Rosenkranz, Susan L PhD‡; Kuritzkes, Daniel R MD§; Eron, Joseph J MD∥; Perelson, Alan S PhD¶; Gerber, John G MD#

JAIDS Journal of Acquired Immune Deficiency Syndromes: 1 July 2005 - Volume 39 - Issue 3 - pp 272-283
Basic Science

We propose a long-term HIV-1 dynamic model by considering drug potency, drug exposure, and drug susceptibility. Using a Bayesian approach, HIV-1 dynamic parameters were estimated by fitting the model to viral load data from a phase 1/2 randomized clinical study of 2 indinavir (IDV)/ritonavir (RTV)-containing highly active antiretroviral (ARV) therapy regimens in HIV-infected subjects who had previously failed protease inhibitor-containing ARV therapies. A large between-subject variation in estimated viral dynamic parameters was observed, even after accounting for variations in drug exposure and drug susceptibility, suggesting that characteristics of HIV-1 dynamics are host dependent. Significant correlations of baseline factors such as HIV-1 RNA levels and CD4+ cell counts with viral dynamic parameters were found. These correlations coincide with biologic interaction mechanisms between HIV and the host immune system and also provide an explanation for the correlations between the baseline viral load and phase 1 viral decay rate, for which inconsistent results have been reported in the literature. The relations between viral dynamic parameters and virologic response were established, and these results suggest that viral dynamic parameters may play an important role in determining treatment success or failure. In particular, we estimated a drug efficacy threshold for each patient that can be used to assess whether an ARV regimen is potent enough to suppress HIV viruses in the individual patient. Our findings indicate that it is necessary to individualize the ARV regimen to treat HIV-1-infected patients. The proposed mathematic models and statistical techniques may provide a framework to simulate and predict antiviral response for individual patients.

From the *Department of Biostatistics and Computational Biology, University of Rochester School of Medicine and Dentistry, Rochester, NY; †Division of Clinical Pharmacology, University of Alabama at Birmingham School of Medicine, Birmingham, AL; ‡Frontier Science and Technology Research Foundation, Boston, MA; §Brigham and Women's Hospital and Partners AIDS Research Center, Harvard Medical School, Cambridge, MA; University of North Carolina at Chapel Hill, Chapel Hill, NC; ¶Theoretical Biology and Biophysics Group, Los Alamos National Laboratory, Los Alamos, NM; and #Divisions of Clinical Pharmacology and Infectious Diseases, University of Colorado Health Sciences Center, Denver, CO.

Received for publication August 13, 2004; accepted March 30, 2005.

Supported by the Adult AIDS Clinical Trials Group of the National Institute of Allergy and Infectious Diseases and National Institutes of Health grants AI32775, AI052765, AI055290, AI38855, AI28433, RR06555, RR00046, and AI50410.

Clinical study A5055 was reviewed and approved by the institutional review boards at each site. Informed consent was obtained from each volunteer patient in accordance with guidelines of the US Department of Health and Human Services and those of the authors' institutions.

Reprints: Hulin Wu, Department of Biostatistics and Computational Biology, University of Rochester School of Medicine and Dentistry, 601 Elmwood Avenue, Box 630, Rochester, NY 14642 (e-mail:

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).

Back to Top | Article Outline


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 (AUC12h), maximum plasma concentration (Cmax), time to Cmax (Tmax), 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 log10 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 log10 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.

Back to Top | Article Outline


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 follows18:

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 mm3) 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 al18 for a detailed discussion). It is important to estimate ec for each patient based on clinical data.

Back to Top | Article Outline

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 (IC50),19 to quantify agent-specific drug susceptibility. To model within-host changes over time in IC50 attributable to the emergence of new drug-resistant mutations, we use the following function18:

where Io and Ir are respective values of IC50(t) at baseline and time point tr at which the resistant mutations dominate. If Ir = Io, no new drug-resistant mutation is developed during treatment. Although more complicated models for IC50 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 IC50 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 kth 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 Emax model for the dose-effect relation follows30:

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


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 IC50 for the 2 agents, respectively. The parameter φ is used to quantify the conversion between in vitro and in vivo IC50 that can be estimated from clinical data. We call φ the IC50 conversion factor. Notice that the larger φ is, the lower is the drug efficacy γ(t). Thus, φ acts like the IC50. Our model (Equation 6) is a variation of the composite Emax models. Note that C12hd 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 al18 and Dixit and Perelson.27 Note that if C12hd, Ad(t), and IC50d(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.

Back to Top | Article Outline


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 ith 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, li}, and Y = {yij, i = 1,···, n; j = 1,···, mi}. Let fij (θi, tj) be the numeric solution for the common logarithmic viral load, log10 (V(t)), of the differential Equation 1 for the ith 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 stages34:

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


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 studies14,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 Giltinan34 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.

Back to Top | Article Outline

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.

Back to Top | Article Outline


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

Back to Top | Article Outline

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.

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 (IC50 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), IC50 had no significant increase from baseline (IC50 = 2.89 and 11.9 ng/mL for IDV and RTV at baseline, respectively) to viral rebound time (IC50 = 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 C12h = 422.2 and 2674.5 ng/mL, IDV AUC12h = 69.4 and 72.1 μg*h/mL, RTV C12h = 761.4 and 2515.9 ng/mL, and RTV AUC12h = 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.

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 IC50 = 11.96 ng/mL and RTV IC50 = 55.04 ng/mL at failure time compared with IDV IC50 = 4.97 ng/mL and RTV IC50 = 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 C12h = 144.3 and 793.3 ng/mL, IDV AUC12h = 12.3 and 20.5 μg*h/mL, RTV C12h = 258.4 and 2978.6 ng/mL, and RTV AUC12h = 31.2 and 69.5 μg*h/mL for these 2 patients, respectively).

Back to Top | Article Outline

Relation Between Baseline Factors and Viral Dynamic Parameters

We have correlated the baseline factors such as baseline viral load (copies/mL), CD4 cells (cells/mm3), 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.

Back to Top | Article Outline

Correlations of Viral Dynamic Parameters With Virologic and Immunologic Responses

We correlated the estimated viral dynamic parameters with the viral load (log10 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.

Back to Top | Article Outline


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 (IC50) 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 al47 and Wu et al10 reported that plasma baseline HIV-1 RNA levels were positively correlated with first-phase viral decay rates. Notermans et al47 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 al6,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 study6 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 al39 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 (log10 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 Wu42 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 IC50 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.

Back to Top | Article Outline


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.

Back to Top | Article Outline


1. Ho DD, Neumann AU, Perelson AS, et al. Rapid turnover of plasma virions and CD4 lymphocytes in HIV-1 infection. Nature. 1995;373:123-126.
2. Wei X, Ghosh SK, Taylor ME, et al. Viral dynamics in human immunodeficiency virus type 1 infection. Nature. 1995;373:117-122.
3. Perelson AS, Neumann AU, Markowitz M, et al. HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time. Science. 1996;271:1582-1586.
4. Perelson AS, Essunger P, Cao Y, et al. Decay characteristics of HIV-1-infected compartments during combination therapy. Nature. 1997;387:188-191.
5. Nowak MA, Bonhoeffer S, Clive L, et al. HIV results in the frame. Nature. 1995;375:193.
6. Wu H, Kuritzkes DR, McClernon DR, et al. Characterization of viral dynamics in human immunodeficiency virus type 1-infected patients treated with combination antiretroviral therapy: relationships to host factors, cellular restoration and virological endpoints. J Infect Dis. 1999;179:799-807.
7. Mueller BU, Zeichner SL, Kuznetsov VA, et al. Individual prognoses of long-term responses to antiretroviral treatment based on virological, immunological and pharmacological parameters measured during the first week under therapy. AIDS. 1998;12(Suppl):F191-F196.
8. Mittler J, Essunger P, Yuen GJ, et al. Short-term measures of relative efficacy predict longer-term reductions in human immunodeficiency virus type 1 RNA levels following nelfinavir monotherapy. Antimicrob Agents Chemother. 2001;45:1438-1443.
9. Wu H, Mellors J, Ruan P, et al. Viral dynamics and their relations to baseline factors and longer-term virologic responses in treatment-naive HIV-1 infected patients receiving abacavir in combination with HIV-1 protease inhibitors. J Acquir Immune Defic Syndr. 2003;32:557-564.
10. Wu H, Lathey J, Ruan P, et al. Relationship of plasma HIV-1 RNA dynamics to baseline factors and virological responses to highly active antiretroviral therapy (HAART) in adolescents (aged 12-22 years) infected through risk behavior. J Infect Dis. 2004;189:593-601.
11. Huang W, De Gruttola V, Fischl M, et al. Pattern of plasma human immunodeficiency virus type 1 RNA response to antiretroviral therapy. J Infect Dis. 2001;183:1455-1465.
12. Polis MA, Sidorov IA, Yoder C, et al. Correlation between reduction in plasma HIV-1 RNA concentration 1 week after start of antiretroviral treatment and longer-term efficacy. Lancet. 2001;358:1760-1765.
13. Acosta EP, Wu H, Walawander A, et al, for the Adult ACTG 5055 Protocol Team. Comparison of two indinavir/ritonavir regimens in treatment-experienced HIV-infected individuals. J Acquir Immune Defic Syndr. 2004;37:1358-1366.
14. Perelson AS, Nelson PW. Mathematical analysis of HIV-1 dynamics in vivo. SIAM Review. 1999;41:3-44.
15. Nowak MA, May RM. Virus Dynamics: Mathematical Principles of Immunology and Virology. Oxford: Oxford University Press; 2000.
16. Callaway DS, Perelson AS. HIV-1 infection and low steady state viral loads. Bull Math Biol. 2002;64:29-64.
17. Perelson AS. Modelling viral and immune system dynamics. Nat Rev Immunol. 2002;2:28-36.
18. Huang Y, Rosenkranz SL, Wu H. Modeling HIV dynamics and antiviral responses with consideration of time-varying drug exposures, sensitives and adherence. Math Biosci. 2003;184:165-186.
19. Molla A, Korneyeva M, Gao Q, et al. Ordered accumulation of mutations in HIV protease confers resistance to ritonavir. Nat Med. 1996;2:760-766.
20. Wainberg MA, Hsu M, Gu Z, et al. Effectiveness of 3TC in HIV clinical trials may be due in part to the M184V substitution in 3TC-resistant HIV-1 reverse transcriptase. AIDS. 1996;10(Suppl):S3-S10.
21. Bonhoeffer S, Lipsitch M, Levin BR. Evaluating treatment protocols to prevent antibiotic resistance. Proc Natl Acad Sci USA. 1997;94:12106-12111.
22. Besch CL. Compliance in clinical trials. AIDS. 1995;9:1-10.
23. Ickovics JR, Meisler AW. Adherence in AIDS clinical trial: a framework for clinical research and clinical care. J Clin Epidemiol. 1997;50:385-391.
24. Wu H, Ding AA. Population HIV-1 dynamics in vivo: applicable models and inferential tools for virological data from AIDS clinical trials. Biometrics. 1999;55:410-418.
25. Ding AA, Wu H. Relationships between antiviral treatment effects and biphasic viral decay rates in modeling HIV dynamics. Math Biosci. 1999;160:63-82.
26. Ding AA, Wu H. Assessing antiviral potency of anti-HIV therapies in vivo by comparing viral decay rates in viral dynamic models. Biostatistics. 2001;2:13-29.
27. Dixit NM, Perelson AS. Complex patterns of viral load decay during antiretroviral therapy: influence of pharmacokinetics and intracellular delay. J Theor Biol. 2004;226:95-109.
28. Dixit NM, Markowitz M, Ho DD, et al. Estimates of intracellular delay and average drug efficacy from viral load data of HIV-infected individuals under antiretroviral therapy. Antivir Ther. 2004;9:237-246.
29. Fitzgerald AP, DeGruttola VG, Vaida F. Modelling HIV viral rebound using non-linear mixed effects models. Stat Med. 2002;21:2093-2108.
30. Gabrielsson J, Weiner D. Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications. Stockholm: Apotekarsocieteten; 2000.
31. Sheiner LB. Modeling pharmacodynamics: parametric and nonparametric approaches. In: Rowland M, Sheiner LB, Steimer JL, eds. Variability in Drug Therapy: Description, Estimation, and Control. New York: Raven Press; 1985:139-152.
32. Hsu A, Isaacson J, Kempf DJ, et al. Trough concentrations-EC50 relationship as a predictor of viral response for ABT-378/ritonavir in treatment-experienced patients [poster session 171]. Presented at: 40th Interscience Conference on Antimicrobial Agents and Chemotherapy; 2000; San Francisco.
33. Kempf DJ, Hsu A, Jiang P, et al. Response to ritonavir intensification in indinavir recipients is highly correlated with virtual inhibitory quotient [abstract 523]. Presented at: Eighth Conference on Retroviruses and Opportunistic Infections; 2001; Chicago.
34. Davidian M, Giltinan DM. Nonlinear Models for Repeated Measurement Data. London: Chapman & Hall; 1995.
35. Gelfand AE, Hills SE, Racine-Poon A, et al. Illustration of Bayesian inference in normal data models using Gibbs sampling. J Am Stat Assoc. 1990;85:972-985.
36. Carlin BP, Louis TA. Bayes and Empirical Bayes Methods for Data Analysis. London: Chapman & Hall; 1996.
37. Perelson AS, Kirschener DE, De Boer R. Dynamics of HIV infection of CD4+ T cells. Math Biosci. 1993;114:81-125.
38. Han C, Chaloner K, Perelson AS. Bayesian analysis of a population HIV dynamic model. In: Gatsonis C, Carriquiry A, Gelman A, et al, eds. Case Studies in Bayesian Statistics, vol. 6. New York: Springer-Verlag; 2002: 223-237.
39. Markowitz M, Louie M, Hurley A, et al. A novel antiviral intervention results in more accurate assessment of human immunodeficiency virus type 1 replication dynamics and T-cell decay in vivo. J Virol. 2003;77:5037-5038.
40. Zhang L, Dailey PD, He T, et al. Rapid clearance of simian immunodeficiency virus particles from plasma of rhesus macaques. J Virol. 1999;73:855-860.
41. Ramratnam B, Bonhoeffer S, Binley J, et al. Rapid production and clearance of HIV-1 and hepatitis C virus assessed by large volume plasma apheresis. Lancet. 1999;354:1782-1785.
42. Huang Y, Wu H. Mechanistic PK/PD modeling of antiretroviral therapies in AIDS clinical trials. In: D'Argenio DZ, ed. Advanced Methods of Pharmacokinetic and Pharmacodynamic Systems Analysis, vol. III. New York: Kluwer Academic Publishers; 2004:221-237.
43. Wakefield JC. The Bayesian approach to population pharmacokinetic models. J Am Stat Assoc. 1996;91:61-76.
44. Wakefield JC, Smith AFM, Racine-Poon A, et al. Bayesian analysis of linear and non-linear population models using the Gibbs sampler. Appl Stat. 1994;43:201-221.
45. Gamerman D. Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference. London: Chapman & Hall; 1997.
46. Venables WN, Ripley BD. Modern Applied Statistics with S-Plus. 3rd ed. New York: Springer; 1999.
47. Notermans DW, Goudsmit J, Danner SA, et al. Rate of HIV-1 decline following antiretroviral therapy is related to viral load at baseline and drug regimen. AIDS. 1998;12:1483-1490.
48. Ding AA, Wu H. Relationships between antiviral treatment effects and biphasic viral decay rates in modeling HIV dynamics. Math Biosci. 1999;160:63-82.
49. Zhang XQ, Schooley RT, Gerber JG. The effect of increasing α1-acid glycoprotein concentration on the antiviral efficacy of human immunodeficiency virus protease inhibitors. J Infect Dis. 1999;180:1833-1837.
50. Shulman N, Zolopa A, Havlir D, et al. Virtual inhibitory quotient predicts response to ritonavir boosting of indinavir-based therapy in human immunodeficiency virus-infected patients with ongoing viremia. Antimicrob Agents Chemother. 2002;46:3907-3916.
51. Boffito M, Back DJ, Blaschke TF, et al. Protein binding in antiretroviral therapies. AIDS Res Hum Retroviruses. 2003;19:825-835.

Cited By:

This article has been cited 1 time(s).

JAIDS Journal of Acquired Immune Deficiency Syndromes
Density-Dependent Decay in HIV-1 Dynamics
Holte, SE; Melvin, AJ; Mullins, JI; Tobin, NH; Frenkel, LM
JAIDS Journal of Acquired Immune Deficiency Syndromes, 41(3): 266-276.
PDF (202) | CrossRef
Back to Top | Article Outline

correlation analysis; drug exposures; drug resistance; dynamic parameter estimation; HIV dynamic modeling; treatment efficacy

© 2005 Lippincott Williams & Wilkins, Inc.