Email updates

Keep up to date with the latest news and content from BMC Medical Research Methodology and BioMed Central.

Open Access Research article

Estimation of dynamical model parameters taking into account undetectable marker values

Rodolphe Thiébaut12*, Jérémie Guedj1, Hélène Jacqmin-Gadda1, Geneviève Chêne2, Pascale Trimoulet3, Didier Neau4 and Daniel Commenges1

Author Affiliations

1 INSERM E0338 Biostatistics, Bordeaux 2 University, Bordeaux, France

2 INSERM U593, Bordeaux 2 University, Bordeaux, France

3 Department of virology, Bordeaux University Hospital, Bordeaux, France

4 Department of infectious disease, Bordeaux University Hospital, Bordeaux, France

For all author emails, please log on.

BMC Medical Research Methodology 2006, 6:38  doi:10.1186/1471-2288-6-38


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2288/6/38


Received:14 February 2006
Accepted:1 August 2006
Published:1 August 2006

© 2006 Thiébaut et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Mathematical models are widely used for studying the dynamic of infectious agents such as hepatitis C virus (HCV). Most often, model parameters are estimated using standard least-square procedures for each individual. Hierarchical models have been proposed in such applications. However, another issue is the left-censoring (undetectable values) of plasma viral load due to the lack of sensitivity of assays used for quantification. A method is proposed to take into account left-censored values for estimating parameters of non linear mixed models and its impact is demonstrated through a simulation study and an actual clinical trial of anti-HCV drugs.

Methods

The method consists in a full likelihood approach distinguishing the contribution of observed and left-censored measurements assuming a lognormal distribution of the outcome. Parameters of analytical solution of system of differential equations taking into account left-censoring are estimated using standard software.

Results

A simulation study with only 14% of measurements being left-censored showed that model parameters were largely biased (from -55% to +133% according to the parameter) with the exception of the estimate of initial outcome value when left-censored viral load values are replaced by the value of the threshold. When left-censoring was taken into account, the relative bias on fixed effects was equal or less than 2%. Then, parameters were estimated using the 100 measurements of HCV RNA available (with 12% of left-censored values) during the first 4 weeks following treatment initiation in the 17 patients included in the trial. Differences between estimates according to the method used were clinically significant, particularly on the death rate of infected cells. With the crude approach the estimate was 0.13 day-1 (95% confidence interval [CI]: 0.11; 0.17) compared to 0.19 day-1 (CI: 0.14; 0.26) when taking into account left-censoring. The relative differences between estimates of individual treatment efficacy according to the method used varied from 0.001% to 37%.

Conclusion

We proposed a method that gives unbiased estimates if the assumed distribution is correct (e.g. lognormal) and that is easy to use with standard software.

Background

Dynamical models based on system of differential equations have been successfully used for a better understanding of the pathogenesis of infectious diseases [1,2]. Two landmark papers appeared in 1995 demonstrating the high turnover of the human immunodeficiency virus (HIV) and infected CD4+ T lymphocytes cells [3,4]. Using such dynamical models, Neumann et al. [5] gave some insight in the effect of interferon based therapy used to treat patients infected by hepatitis C virus (HCV). Moreover, the estimate of the percentage of virus production blocked by the therapy is now widely used in this field [6-10] to evaluate the efficacy of treatment regimens in various contexts such as patients co-infected with HIV and HCV.

Although dynamical models parameters such as virus clearance or treatment efficacy are very useful, their estimation is most often performed for individual subjects separately. The limitations of such statistical approach as well as the interest of hierarchical models have already been underlined [11,12]. The main advantage of hierarchical models (also called mixed/random effects models) is their ability to estimate all parameters at the same time, using all available data even in case of unbalanced data, i.e. the number of measurements can vary from one patient to another. Parameters can be estimated using a Bayesian approach [13-15] or other approaches [16]. Another advantage working with analytical solutions of the system of differential equations is that standard softwares for non linear mixed models can be used [17].

Nonetheless, a major problem arises when using viral load data. The assays used to quantify HIV or HCV RNA are limited by a detection threshold that may lead to undetectable values when the true viral load is below this threshold. From a clinical point of view, the aim of any treatment is to reduce the viral load as much as possible [18]. Therefore, the practical definition of virological response is the occurrence of sustained undetectable values. The threshold of undetectable values is changing with the improvement of the assays for quantifying the viral load. When analysing viral load as a continuous variable, the left-censored measures are most often analyzed by replacing their value by an arbitrary value (e.g. threshold or half of the threshold). Although the sensitivity of the assays is improving, this limitation still persists and has already been underlined in the context of dynamical models [12,15]. Methods to take into account left-censored repeated measures in linear mixed models [19-21] or in non linear mixed models [15] have already been proposed. In this paper, we show how such an approach can be implemented using standard software in the case of non linear mixed models. Furthermore, we evaluate the impact of not taking into account undetectable values when studying HCV dynamics in the context of a phase II randomised clinical trial for the treatment of HCV infection in HIV co-infected patients.

Methods

Study example

The motivating study was a phase II randomised clinical trial evaluating the efficacy of pegylated-interferon (PEG-IFN)-α2a and Ribavirin (RBV) for the treatment of HCV infection among 17 HIV co-infected patients who had already been treated for HCV [22]. HCV RNA quantification was performed at least three times within the first 4 weeks (W): W0 (treatment initiation), W2 and W4. In 8 patients, blood samples were collected more intensively with additional measures at 6 hours (H6), H12, day 1 (D1), D2, D4, W1 and W3. Patients were followed until W72 for final evaluation of the virological response but the study of viral dynamics was restricted to the first 4 weeks because of model assumptions (see below). The concentration of plasma HCV RNA was determined using a quantitative reverse transcription polymerase chain reaction (RT-PCR) assay (Cobas Amplicor HCV Monitor Test, version 2.0; Roche Molecular Systems). The lower detection limit of this assay was 600 IU/mL, i.e. 2.78 log10 IU/mL. Of note, one international unit (IU) equals approximately 2.2 copies/mL.

Mathematical model

The model used to estimate HCV dynamics was first described by Neumann et al. [5] with the following differential equations:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M1">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M2">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M3">View MathML</a>

where T is the number of target cells (i.e. hepatocytes), I the number of productively infected cells and V the plasma HCV viral load. Target cells are produced at rate s (per day) and die at rate μ. The number of cells which become infected per day is proportional to the number of circulating virions and available target cells with a proportionality constant β (infection rate). Infected cells die at rate δ per day. HCV virions are produced at a rate p per infected cells per day and are cleared at a rate c per day. In the present model, the HCV treatment is supposed to reduce the production of virions from infected cells by a fraction (1-ε). The possible effect of IFN as well as RBV on de novo rate of infection [5] or on infectivity by producing a fraction of non-infectious virions [23] have been discussed. For the purpose of this paper, we assumed only a combined effect of both drugs on production rate of new virions because this measure was the most widely used by other investigators [6-9].

When working on a short period of 2–4 weeks, it sounds reasonable to consider that the number of uninfected hepatocytes (T) remains constant (equal to the baseline value) because of the slow turnover of these cells [5]. Therefore, assuming a pre-treatment steady-state, the analytical solution of the equations (2) and (3) with T constant is:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M4">View MathML</a>

for t>t0, where

λ1 = {(c + δ) + <a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M5">View MathML</a>}, λ2 = {(c + δ) - <a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M5">View MathML</a>} and <a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M6">View MathML</a>. The viral decay is assumed to begin at t0 = 0.25 day (6 hours), corresponding to the drugs pharmacokinetics [5].

Hierarchical formulation

The previous notations do not account for patient/measurement level. Most often parameters of such models are estimated patient by patient assuming Gaussian, homoskedastic measurement error. A more valid and powerful approach is based on a hierarchical formulation of the model [11] that can distinguish at least two levels of variation. Hence, for the jth measurement of a subject i performed at a time tij:

- Stage 1: intra-patient variation

yij = log10(V(tij, θi)) + e

with <a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M7">View MathML</a>

The outcome is the logarithm (base 10) of the true viral load (function of tij and θi, the p-vector of model parameters) plus a Gaussian measurement error e. Ini is a identity matrix of dimension ni × ni, ni being the number of measurements available for the subject i.

- Stage 2: inter-patient variation

θi = θ + γi

with γi ~ MVN(0, D)

θ = [V0, ε, δ, c] is the p-vector of average (fixed) effect in the whole study population and γi is a q-vector (q ≤ p) of random effects for correcting θ for each subject (random effect). Actually, θ is a log-transformation of original parameters that have several advantages including a positivity constraint for original parameters. Random effects γi were assumed to be normally distributed with a variance-covariance matrix D. θi are estimated through Empirical Bayes estimates.

Model likelihood

As presented in more details elsewhere [24], the method proposed to take into account left-censored values when estimating parameters is to maximise a full likelihood distinguishing the contribution of observed measures (<a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M8">View MathML</a> for j = 1,...<a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M9">View MathML</a>) and left-censored measures (<a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M10">View MathML</a> for j = 1,...<a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M11">View MathML</a>) of viral load. The likelihood can be written:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M12">View MathML</a>

with <a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M13">View MathML</a> being the univariate normal density of <a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M8">View MathML</a> given the random effects γi and <a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M14">View MathML</a> is the cumulative distribution function of the normal distribution of <a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M10">View MathML</a> given the random effects. The calculation of this likelihood leads to the integration over u = u1,u2,...uq, that is a multiple integral of dimension q. Therefore, with this method, rather than imputing a fixed value of undetectable viral load, one assumes that left-censored values are completing the Gaussian distribution of Yi. A crude approach assumes that left-censored values contribute like observed values, being equal to the value of the threshold or any other given value. In this case, the likelihood is simply:

<a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M15">View MathML</a>

Algorithm and implementation

The maximisation of the likelihood function can be performed with standard software such as NLMIXED in SAS® [24]. Using this procedure, the default algorithm is a quasi-newton algorithm and the calculation of the multiple integral is performed by adaptative quadrature. An example of code used for this paper is provided in appendix.

Simulation study

Simulations were performed to compare the bias on parameter estimates when taking into account left-censoring or not. Using the analytical solution (4) and allowing a random individual variation for the initial viral load and treatment efficacy, parameters to estimate were: ⌊V0, ε, δ c, <a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M16">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M17">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M18">View MathML</a>

with V0i = V0 + γ0i, εi = ε + γ1i and <a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M19">View MathML</a>

To constraint parameters to be in the correct range, estimations were performed on transformed parameters (for the study on real data, as well) using a logarithm function for δ, c and logit function for ε (bounding ε between 0 and 1). In the simulation study, we fixed t0 = 0 but results were similar with t0 = 0.25.

Values for model parameters were defined according to the results reported in the literature of HCV dynamics [6]. In our application where patients were previously treated and dually infected by HIV and HCV, the estimate of treatment efficacy is less than those usually reported in naïve patients mono-infected with HCV [5].

The steps followed for the simulations were:

1) Sample V0i = V0 0i and εi = ε + γ1i for a subject i

2) Simulate the differential equations (1)-(3) model and keep measures at the time points: keep measures at H0, H6, H12, D1, D2, D4, W1, W2, W3 and W4. Left-censor measures below 2.78 IU/mL.

3) Repeat N times (for N = 20 subjects) steps 1 and 2

4) Estimate parameters with (5) when taking into account left-censoring and with usual likelihood (6) replacing left-censored values by the value of the threshold, i.e. 2.78 IU/mL.

5) Calculate the relative bias for each parameter RB = 100*(estimate-true value)/true value

6) Repeat 1000 times steps 1 to 5 and average the relative biases

Results

Simulation study

Results of the simulations are shown in Table 1. On average, 14% of simulated measures of HIV RNA were left-censored when the treatment efficacy was ε = 80%. Crude estimates provided by standard likelihood (6) maximisation, replacing left-censored values by the value of the threshold, were dramatically biased with the exception of V0. In particular, with only 14% of left-censored measures, infected cells death rate δ and clearance of virus c were underestimated by 55 and 45 percent, respectively. Treatment efficacy (ε) was overestimated by 16%. Variances of random effects were also significantly biased: +20% and -19% for random effects on V0, and ε, respectively. The residual variance was overestimated (+133%).

Table 1. Relative bias of model parameters using non linear mixed models taking into account left-censored (undetectable) values (corrected) or not (crude) with simulated data. N = 20 patients, 1000 simulations, 14% left-censored measures in average.

When taking into account left-censoring, the relative bias on estimates was ≤ 2 % for all parameters but variance parameters. However, the biases on variance parameters significantly decreased (e.g. the bias on <a onClick="popup('http://www.biomedcentral.com/1471-2288/6/38/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/6/38/mathml/M20">View MathML</a> changed from -14% to -0.4%) when increasing the number of subjects included in the sample (e.g. N = 200).

Application

Model parameters were estimated using the HCV RNA data available during the first 4 weeks following treatment initiation in the 17 patients included in the ROCO2 trial. Among these 100 available measurements, 12 were undetectable, i.e. left-censored.

Estimates of parameters taking into account left-censoring or not are shown in Table 2. Differences between estimates according to the method used were large on the death rate infected cells, δ.

Table 2. Estimates of model parameters using non linear mixed models taking into account left-censored (undetectable) values (corrected) or not (crude) with data from ROCO2 clinical trial. N = 17 patients, 100 measures, 12% left-censored.

Using the crude approach, the estimate was 0.13 day-1 (95% confidence interval [CI]: 0.11; 0.17) compared to 0.19 day-1 (CI: 0.14; 0.26) when taking into account left-censoring. These estimates correspond to half-life (t1/2) of infected cells of 5.3 days (t1/2 = ln(2)/δ) and 3.6 days, respectively. Differences between estimates for the other fixed parameters were less important. Furthermore, the confidence intervals of the estimates were larger when taking into account left-censoring (Table 2).

The impact of the method used to estimate the parameters on individual viral load predictions is illustrated in Figure 1. For the first three patients (102, 108 and 201), the decrease of the second part of the viral load slope was more pronounced when taking into account left-censoring. Actually, left-censoring tended to occur on the last measurements depending on the treatment efficacy and the baseline level of viral load. The apparent discrepancy with observed values is obviously due to the fact that undetectable values are plotted on the detection limit (2.78 IU/mL) although the true value is below this threshold. This result is expected as the slope after the shoulder is proportional to the infected cell death rate (δ). As expected for the last patient (206) who did not have any undetectable viral load within 4 weeks, both predictions were very close.

thumbnailFigure 1. Observed and predicted HCV RNA values in four patients. Predictions came from non linear mixed effect models taking into account left-censored (undetectable) values (right side) or not (left side). Observed undetectable HCV RNA measures (<2.78 log10 IU/mL) are plotted at 2.78 log10 IU/mL.

The relative differences between estimates of individual treatment efficacy (εi = ε + γ1i) according to the method used varied from 0.001% to 37%. As expected from simulation results where treatment efficacy tended to be overestimated with the crude approach and from the estimate of the average (fixed) effect ε, the estimated effect was most often higher with the crude approach compared to the corrected one. For instance, the estimate of treatment efficacy in patient 101, was 36% and 45% when taking into account left-censoring or not, respectively. On the contrary, for the patient 201, the estimates were 97% and 93%, respectively. Of note, the model was able to predict viral load changes for the patient 201 thanks to the information provided by the other patients with more numerous measurements available. This is an illustration of the advantage of hierarchical models.

Discussion

In this paper, the impact of taking into account left-censored (undetectable) HCV RNA values was illustrated on the estimation of dynamical models based on a system of differential equations. Although, the proportion of undetectable values was quite low (12%), there were clinically significant differences, particularly in estimate of mean half-life and individual treatment efficacy. Such a result is important because all these parameters are of interest. Treatment efficacy evaluation through dynamical model is broadly used in HCV infection for instance.

We observed smaller biases from the crude approach applied to the real dataset compared to simulation results. However, some parameters values were different to those used in the simulations such as δ (0.13 vs. 0.40). Simulations using values estimated with real data led to smaller biases as observed in the present application (data not shown). The overestimation of the treatment efficacy by the crude approach may appear counter-intuitive because the imputation of the value of the threshold artificially limits the decrease of viral load. However, it is difficult to anticipate the impact of left-censoring in dynamical models because of the complex relationship between parameters, particularly between ε and δ [23]. In the present study, the imputation of the value of threshold level to undetectable viral loads led to a higher level of HCV RNA than the truth, particularly in the second part of the dynamics. The death rate of infected cells (δ) is one of the main parameters influencing viral load levels in this period [5,25]. This explains the underestimation of this parameter. On the other hand, an overestimation of treatment effect on viral production (ε) is needed to obtain a trajectory compatible with the first part of the viral dynamics (high viral load without left-censored measures), given a high infected cells death and virions clearance.

Half-life of infected cells helps in understanding how high is the turnover [3,26,27]. Previously published results [25] can be used to illustrate the size of the impact of left-censoring on HIV infected cells turnover. Differences in estimates of half-life of infected long-lived cells as large as those we reported in HCV would lead to halve the time needed to treat to achieve virus eradication (assuming no viral reservoir). Compared to results with piecewise linear mixed models commonly used with surveillance data (monthly to 6-months intervals between measurements) of HIV RNA [19,20], the estimates of the parameters are more sensitive to undetectable values in the context of dynamical models with highly repeated measurements. Moreover, confidence intervals of estimates were larger when taking into account left-censoring compared to simple imputation that tends to artificially decrease the variability, as previously reported with linear models [20].

The method presented in this paper is easy to implement in standard software. One limitation is that it is based on analytical solutions of the system of differential equations. However, looking at the applied papers on HCV infection, the authors used most often the same model with the same assumptions leading to the same analytical solution. Using hierarchical models taking into account left-censoring should improve the validity of estimation and may help in case of convergence difficulty when using individual data [9]. More complex mathematical models have been proposed to fit additional markers such as liver enzymes level [28] or pharmacokinetics data [29]. In this case, more general approaches based directly on numerical solution of the differential equations should be used [13,15]. Another limitation of the proposed methods is the assumption of log-normal distribution of viral load measures. In our experience, it is most often a reasonable assumption in the case of circulating HIV virus and this could be checked from residuals [20,30]. However, if this assumption is not tenable, extensions based on mixture distributions (log-normal and binary) can be used and are also easily implementable in software [21].

Conclusion

Imputing a single value to left-censored measures of viral load is a wrong assumption and stronger than assuming a given distribution for the whole measurements. We proposed a method that gives unbiased estimates if the assumed distribution is correct (e.g. lognormal) and that is easy to use with standard software.

Competing interests

The trial was supported by a grant from Roche Laboratories.

Authors' contributions

RT carried out the simulations and drafted the manuscript. JG participated to the work of estimation (with RT). JG, HJG and DC participated in the statistical analysis and helped to draft the manuscript. GC, DN and PT performed the clinical trial and helped to draft the manuscript. All authors read and approved the final manuscript.

Appendix

Example of code using NLMIXED to fit the model presented in the methods section taking into account left censoring.

proc nlmixed data = roco2 OPTCHECK;/* option for checking convergence at the optimum */

/* declare the model parameters to estimate */

   parms beta0 = 10 beta1 = -1.0 beta2 = 1 beta3 = 0.8 s2b0 = 1 s2b3 = 0.1 s2 = 0.1;

/* declare constraints for variance parameters */

   bounds s2,s2b0,s2b3 > 0;

pi = 2*arsin(1);

/* model definition */

      V0 = exp(beta0+b0);

      d = exp(beta1);

      c = exp(beta2);

      e = beta3+b3;

      t0 = 0.25;/* 6 hours */

      th = sqrt((c-d)*(c-d)+4*(1-e)*c*d);

      l1 = 0.5*(c+d+th);

      l2 = 0.5*(c+d-th);

      if tps le t0 then pred = V0;

      if tps gt t0 then

      pred = 0.5*V0*((1-(c+d-2*e*c)/th) * exp(-l1*((tps-t0)))+

         (1+(c+d-2*e*c)/th) * exp(-l2*((tps-t0))));

      logpred = log10(pred);

/* likelihood contribution according to the observed/censored status */

* observed ;

if detec = 1 then ll = (1/(sqrt(2*pi*s2)))

         *exp(-(logCV-logpred)**2/(2*s2));

* censored ;

if detec = 0 then ll = probnorm((logCV-logpred)/sqrt(s2));

L = log(ll);

      model logCV ~ general(L);

/* definition of the random effects */

      random b0 b3 ~ normal([0,0], [s2b0,0,s2b3]) subject = id;

Example of code used for simulating data from dynamical model

%do sim = 1 %to &S;

%do id = 1 %to &N;

Data_null_;

logCV0 = 6.16+0.70*rannor(-1);

CV0 = 10**(logCV0);call symput('CV0',CV0);

kmax = 1.39+1.64*rannor(-1);

e = exp(kmax)/(1+exp(kmax));call symput('e',e);

run;

Data sim; do time = 0 to 672 by 1;output;end; run;

Proc model data = sim;

dependent T I CV ;

parm b 0.00000003 d 0.0167 e &e p 4.16 c 0.0833;

if time = 0 then do;

CV = &CV0;

T = (c*d)/(p*b);

I = (c*CV)/p;

      end;

if time ne 0 then do;

dert.T = 0;

dert.I = b*CV*T-d*I;

dert.CV = (1-e)*p*I-c*CV;

end;

solve T I CV/dynamic out = simul;

run ; quit;

Data pat;set simul;tps = time/24;id = &id;CV0 = &CV0;e = &e;

if round(time) in (0,6,12,24,48,96,168,336,504,672);run;

%if &id = 1 %then %do;

Data file; set pat;error = 0.2*rannor(-1);if CV gt 0 then logCV = log10(CV)+error;run;

%end;

%else %do;

Data file; set file pat;error = 0.2*rannor(-1);if CV gt 0 then logCV = log10(CV)+error;run;

%end;

/* truncation */

Data file; set file;

      if logCV lt 2.778 then do;

         logCV = 2.778;detec = 0;end;

      else detec = 1;

run;

ods exclude none;

%end;/* end of patients */

Acknowledgements

The authors thank the investigators and the patients of the ROCO2 clinical trials.

References

  1. Perelson AS: Modelling viral and immune system dynamics.

    Nature Reviews Immunology 2002, 2:28-36. PubMed Abstract | Publisher Full Text OpenURL

  2. Perelson AS, Herrmann E, Micol F, Zeuzem S: New kinetic models for the hepatitis C virus.

    Hepatology 2005, 42:749-754. PubMed Abstract | Publisher Full Text OpenURL

  3. Ho DD, Neumann AU, Perelson AS, Chen W, Leonard JM, Markowitz M: Rapid turnover of plasma virions and CD4 lymphocytes in HIV-1 infection.

    Nature 1995, 373:123-126. PubMed Abstract | Publisher Full Text OpenURL

  4. Wei X, Ghosh SK, Taylor ME, Johnson VA, Emini EA, Deutsch P, Lifson JD, Bonhoeffer S, Nowak MA, Hahn BH, et al.: Viral dynamics in human immunodeficiency virus type 1 infection.

    Nature 1995, 373:117-122. PubMed Abstract | Publisher Full Text OpenURL

  5. Neumann AU, Lam NP, Dahari H, Gretch DR, Wiley TE, Layden TJ, Perelson AS: Hepatitis C viral dynamics in vivo and the antiviral efficacy of interferon-alpha therapy.

    Science 1998, 282:103-107. PubMed Abstract | Publisher Full Text OpenURL

  6. Torriani FJ, Ribeiro RM, Gilbert TL, Schrenk UM, Clauson M, Pacheco DM, Perelson AS: Hepatitis C virus (HCV) and human immunodeficiency virus (HIV) dynamics during HCV treatment in HCV/HIV coinfection.

    The Journal of Infectious Diseases 2003, 188:1498-1507. PubMed Abstract | Publisher Full Text OpenURL

  7. Layden-Almer JE, Ribeiro RM, Wiley T, Perelson AS, Layden TJ: Viral dynamics and response differences in HCV-infected African American and white patients treated with IFN and ribavirin.

    Hepatology 2003, 37:1343-1350. PubMed Abstract | Publisher Full Text OpenURL

  8. Talal AH, Shata MT, Markatou M, Dorante G, Chadburn A, Koch R, Neumann AU, Ribeiro RM, Perelson AS: Virus dynamics and immune responses during treatment in patients coinfected with hepatitis C and HIV.

    Journal of Acquired Immune Deficiency Syndromes 2004, 35:103-113. PubMed Abstract | Publisher Full Text OpenURL

  9. Sherman KE, Shire NJ, Rouster SD, Peters MG, Koziel MJ, Chung RT, Horn PS: Viral kinetics in hepatitis C or hepatitis C/human immunodeficiency virus-infected patients.

    Gastroenterology 2005, 128:313-327. PubMed Abstract | Publisher Full Text OpenURL

  10. Herrmann E, Zeuzem S, Sarrazin C, Hinrichsen H, Benhamou Y, Manns MP, Reiser M, Reesink H, Calleja JL, Forns X, Steinmann GG, Nehmiz G: Viral kinetics in patients with chronic hepatitis C treated with the serine protease inhibitor BILN 2061.

    Antiviral Therapy 2006, 11:371-376. PubMed Abstract OpenURL

  11. Wu HL, 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. PubMed Abstract | Publisher Full Text OpenURL

  12. Wu HL: Statistical methods for HIV dynamic studies in AIDS clinical trials.

    Statistical Methods in Medical Research 2005, 14 (2):171-192. Publisher Full Text OpenURL

  13. Wu HL, Ding AA, DeGruttola V: Estimation of HIV dynamic parameters.

    Statistics in Medicine 1998, 17:2463-2485. PubMed Abstract | Publisher Full Text OpenURL

  14. Putter H, Heisterkamp SH, Lange JM, de Wolf F: A Bayesian approach to parameter estimation in HIV dynamical models.

    Statistics in Medicine 2002, 21:2199-2214. PubMed Abstract | Publisher Full Text OpenURL

  15. Banks HT, Grove S, Hu S, Ma Y: A hierarchical Bayesian approach for parameter estimation in HIV models.

    Inverse Problems 2005, 21:1803-1822. Publisher Full Text OpenURL

  16. Kuhn E, Lavielle M: Maximum likelihood estimation in nonlinear mixed effects models.

    Computational Statistics & Data Analysis 2005, 49:1020-1038. Publisher Full Text OpenURL

  17. SAS Institute Inc : The NLMIXED Procedure. In SAS/STAT User's Guide, Version 8. Carry, NC, SAS Institute Inc.; 2000:2419-2504. OpenURL

  18. Alberti A, Clumeck N, Collins S, Gerlich W, Lundgren J, Palu G, Reiss P, Thiébaut R, Weiland O, Yazdanpanah Y, Zeuzem S: Short statement of the first European consensus conference on the treatment of chronic hepatitis B and C in HIV co-infected patients.

    Journal of Hepatology 2005, 42:615-624. PubMed Abstract | Publisher Full Text OpenURL

  19. Hughes JP: Mixed effects models with censored data with application to HIV RNA levels.

    Biometrics 1999, 55:625-629. PubMed Abstract | Publisher Full Text OpenURL

  20. Jacqmin-Gadda H, Thiébaut R, Chêne G, Commenges D: Analysis of left-censored longitudinal data with application to viral load in HIV infection.

    Biostatistics 2000, 1:355-368. PubMed Abstract | Publisher Full Text OpenURL

  21. Berk KN, Lachenbruch PA: Repeated measures with zeros.

    Statistical Methods in Medical Research 2002, 11:303-316. PubMed Abstract | Publisher Full Text OpenURL

  22. Breilh D, Neau D, Djabarouti S, Trimoulet P, Pellegrin JL, Duprat C, Ragnaud JM, Dupon M, Saux MC: Plasma Target Concentration of Ribavirin in HCV/HIV Co-infected Patients: ; Boston, February 22-25, 2005. ; 2005.

  23. Dixit NM, Layden-Almer JE, Layden TJ, Perelson AS: Modelling how ribavirin improves interferon response rates in hepatitis C virus infection.

    Nature 2004, 432:922-924. PubMed Abstract | Publisher Full Text OpenURL

  24. Thiébaut R, Jacqmin-Gadda H: Mixed models for longitudinal left-censored repeated measures.

    Comput Methods Programs Biomed 2004, 74:255-260. PubMed Abstract | Publisher Full Text OpenURL

  25. Perelson AS, Essunger P, Cao Y, Vesanen M, Hurley A, Saksela K, Markowitz M, Ho DD: Decay characteristics of HIV-1-infected compartments during combination therapy.

    Nature 1997, 387:188-191. PubMed Abstract | Publisher Full Text OpenURL

  26. Emery VC, Cope AV, Bowen EF, Gor D, Griffiths PD: The dynamics of human cytomegalovirus replication in vivo.

    Journal of Experimental Medicine 1999, 190:177-182. PubMed Abstract | Publisher Full Text OpenURL

  27. Whalley SA, Murray JM, Brown D, Webster GJ, Emery VC, Dusheiko GM, Perelson AS: Kinetics of acute hepatitis B virus infection in humans.

    Journal of Experimental Medicine 2001, 193:847-854. PubMed Abstract | Publisher Full Text OpenURL

  28. Ribeiro RM, Layden-Almer J, Powers KA, Layden TJ, Perelson AS: Dynamics of alanine aminotransferase during hepatitis C virus treatment.

    Hepatology 2003, 38:509-517. PubMed Abstract | Publisher Full Text OpenURL

  29. Powers KA, Dixit NM, Ribeiro RM, Golia P, Talal AH, Perelson AS: Modeling viral and drug kinetics: hepatitis C virus treatment with pegylated interferon alfa-2b.

    Seminars in Liver Disease 2003, 23 Suppl 1:13-18. PubMed Abstract | Publisher Full Text OpenURL

  30. Lyles RH, Lyles CM, Taylor DJ: Random regression models for human immunodeficiency virus ribonucleic acid data subject to left censoring and informative drop-outs.

    Journal of the Royal Statistical Society: Series C 2000, 49:485-497. Publisher Full Text OpenURL

Pre-publication history

The pre-publication history for this paper can be accessed here:

http://www.biomedcentral.com/1471-2288/6/38/prepub