Influence of a priori Information, Designs, and Undetectable Data on Individual Parameters Estimation and Prediction of Hepatitis C Treatment Outcome.

Hepatitis C viral kinetic analysis based on nonlinear mixed effect models can be used to individualize treatment. For that purpose, it is necessary to obtain precise estimation of individual parameters. Here, we evaluated by simulation the influence on Bayesian individual parameter estimation and outcome prediction of a priori information on population parameters, viral load sampling designs, and methods for handling data below detection limit (BDL). We found that a precise estimation of both individual parameters and treatment outcome could be obtained using as few as six measurements in the first month of therapy. This result remained valid even when incorrect a priori information on population parameters was set as long as the parameters were identifiable and BDL data were properly handled. However, setting wrong values for a priori population parameters could lead to severe estimation/prediction errors if BDL data were ignored and not properly accounted in the likelihood function.CPT: Pharmacometrics & Systems Pharmacology (2013) 2, e56; doi:10.1038/psp.2013.31; published online 17 July 2013.

Chronic infection with hepatitis C virus (HCV) is a liver disease that affects about 150 million people worldwide and is directly responsible of about 350,000 deaths every year.[1] The goal of anti-HCV treatment is to achieve a sustained virologic response (SVR), defined as undetectable serum HCV RNA 24 weeks after treatment cessation.[2] HCV is classified into six major genotypes (GT) with HCV GT-2/3 being the second cause of chronic hepatitis C (after GT-1), accounting for ~15–20% of infection in Western countries.[3]

Since 2001, the combination of pegylated interferon (peg-IFN) and ribavirin (RBV) is the backbone of anti-HCV treatment with SVR rate of ~50 and 80% in patients infected with GT-1 and GT-2/3, respectively.[4,5,6] In 2011, the approval of two protease inhibitors marked a new era of HCV therapy with a dramatic improvement in SVR rates in patients infected with HCV GT-1.[7,8,9,10,11,12,13] However, there is no clear evidence that PIs are beneficial in GT-2/3 patients.[14,15,16,17] Even though new treatment, such as nucleotide analogs may be effective against GT-2/3,[18] their cost and the fact that peg-IFN/RBV is already efficient makes bitherapy likely to remain essential in the treatment against HCV GT-2/3.[15,17] Because peg-IFN/RBV therapy is associated with several significant side effects and high costs,[19] several efforts have been made to evaluate the possibility of treatment individualization.[20]

For that purpose, one can use viral kinetic models whose parameters have a high predictive value of treatment outcome.[21,22] However, the use of these models is limited by the fact that frequent viral load data in the first weeks following treatment initiation are required to obtain precise estimation of the parameters. One way to improve the precision of individual parameter is to consider that the population parameters are known a priori and to perform Bayesian estimation of individual parameters. Thus, this method combines a priori information of population parameters gathered from previous studies and individual viral load data prospectively obtained in a patient. This approach is similar to what is done in therapeutic drug monitoring using population pharmacokinetic models.[23,24]

However, the relevance of this approach is still contingent on the study design.[25] In practice, the difficulty to frequently assess viral load levels often gives predictions based on a limited number of viral load data within each patient. We would therefore like to evaluate the quality of individual parameter estimation and treatment outcome prediction using a realistic design based on a small number of short-term observations.

A common challenge in analyzing HCV kinetic data is the fact that a large proportion of viral load data are below the detection limit (BDL). Several studies have shown that naive approaches that omit or impute BDL data at an arbitrary value led to biased population parameter estimates, and this can be corrected by taking BDL data into account in the likelihood function.[26,27,28,29] However, these studies focused on the population parameters and did not evaluate whether and how these methods improved Bayesian individual parameter estimation.

Here, our goal is to evaluate, by simulation and in the context of HCV GT-2/3, the influence of a priori population model, of viral load sampling designs and of methods for handling BDL data when estimating individual parameters and treatment response.

Results

Description of the simulated data

The percentages of BDL data were equal to 57.4, 27.8, 37.7, and 38.5 with designs D, D, D, and D, respectively (see “Methods”). The SVR rate was equal to 78.4%, consistent with the rate of 80% reported in the literature for HCV GT-2/3.[5,6] Of note, 53.6% of the simulated patients had already eradicated the virus by week 12 (W12) of treatment. Non-SVR patients were classified as relapsers (undetectable HCV RNA at the end of treatment becomes detectable 24 weeks later), null-responders (HCV RNA declines from baseline <1 log10 at all time), or rebounders (HCV RNA declines from baseline >1 log10 then rebounds under treatment), accounting for 7.6, 2.3, and 11.7% of our simulated population, respectively. Furthermore, two types of rebounders were considered: patients having a breakthrough (rebounds after a period being undetectable) or a partial virologic response (rebounds after a declining period but never being undetectable), accounting for 3.3 and 8.4% of the simulated population, respectively.

Individual parameter estimation

gives the mean relative bias (MRE), relative root mean square error (RRMSE), and shrinkage of the parameters β, δ, c, and ε obtained for each design using true (M) or false population models (M, M) and methods that omit (OMIT) or handle BDL data in the likelihood function (LBA).

Even when the a priori information on population parameters was correct (i.e., M), individual parameters were estimated with some bias and with the most informative design (D), RRMSE for the infectivity rate β was still equal to 55.4%. Moreover, β consistently had an extremely high shrinkage, >80% regardless of designs ( and ). This is due to the fact that this parameter can be precisely estimated only if there is a virologic rebound, a feature that was observed in only 11.7% of patients. The three other parameters c, δ, and ε could be estimated with a shrinkage and a RRMSE <50% as long as viral load data in the first week of treatment were available. As expected, the precision of parameter estimation deteriorates with less frequent sampling designs and the RRMSE for δ, c, and log10(1-ε), increased by ~30, 60, and 15%, respectively, when using the D as opposed to the D design. Similarly, the shrinkage for δ, c, and ε increased by ~80, 90, and 300%, respectively, when using the D as opposed to the D design. Surprisingly, a correct handling of BDL data did not critically improve the estimation of individual parameters when true a priori population parameters were set: for all parameters, only a small difference was found between LBA and OMIT methods using MRE (<6%) and RRMSE (<4%) as evaluation criteria, regardless of designs.

Setting the fixed effect of δ and ε at false values (M) naturally led to larger estimation errors, especially for δ and c. Since in this model the a priori value for δ was less than half of the true value (0.14 vs. 0.32 day−1, respectively), the individual values for this parameter were consistently underestimated with an increasing relative bias, compared with that obtained with true information: for instance, using the design D, relative bias going from 3.9% with true a priori population parameters up to 20.7% with false a priori information even when correctly handling BDL data. As the rate of the first phase of viral decline is proportional to c × ε (see “Methods”), the estimates of these two parameters were correlated. This explains why setting a lower ε than the one used to generate the data (90% instead of 99.6%) led to an overestimation of c. This was particularly visible when considering designs without frequent samplings where c is poorly identifiable. For instance, the relative bias increased from 6% to 28% and RRMSE increased from 45% to 63% when comparing the D design under the true and false model. However, the quality of parameter estimates remained reasonable (RRMSE < 50%) when using designs with frequent sampling in the first week, at least as long as BDL data were correctly handled. Indeed, and unlike what was seen when using the true population parameters, omission of BDL data dramatically deteriorated parameter precision. This was particularly noticeable for δ where MRE doubled when omitting instead of correctly handling BDL data: for instance, using D design, MRE increased from −8.6% with LBA method to −15.7% when BDL data were omitted.

However, the result that setting false population parameters is not critical to achieve precise parameter estimates (at least as long as BDL data are correctly handled) becomes wrong if parameters are poorly identifiable. Indeed, extremely high errors were obtained for all parameters when setting β to a value 10-fold larger than the true one (M, ).

Prediction of treatment outcome

We then investigated the impact of individual parameter estimation on treatment outcome prediction for the simulated patients. Results are presented in . The influence of simulated individual parameters on treatment outcome was assessed in Supplementary Material 1 online.

In case of correct a priori information, the prediction of individual treatment outcome was very good: even with the most challenging design where individual parameter estimates were biased, the misclassification rate was lower than 6%. The choice of designs did not have much impact on the prediction of treatment outcome although a slightly better prediction was obtained with rich design. With less information (by shortening or reducing sampling frequency), more false-positive (FP) responses were obtained as the prediction tended to be influenced by the response of the profile using the mean parameter values, i.e., SVR with this parameter regime (Supplementary Material 2 online) Of note, the relevance of this approach can be examined by comparing misclassification rate using viral kinetic model with what would be obtained if no information was available. In that case, one would predict SVR in all patients, leading to an error rate of 21.6%, in which all are FP.

In case of false a priori information about δ and ε, we found, as previously, that a correct handling of BDL data was necessary to obtain good prediction and omitting BDL data dramatically increased false-negative (FN) responses (Supplementary Material 2 online), especially when designs become less informative: for instance, in comparison with a correct handling of BDL data, FN increased from 1.6% to 23.5% and 0.5% to 71.5% when omitting BDL data in designs D and D.

If one would predict SVR using rapid virologic response (RVR), defined as undetectable HCV RNA at W4 after treatment initiation, the misclassification rate would be still very good (6.9%). This, actually, means that the use of a viral kinetic model and data in the first month improved treatment outcome prediction by only 1.4% and 1.2% if a priori information was true or not, respectively. Thus the use of viral kinetic models does not bring, with this parameter regime, a critical benefit. This may be due to the fact that we evaluated the ability of treatment outcome prediction in a very favorable case where treatment outcome was highly correlated with the virologic response at W4. This is why we also compared the predictive ability of RVR and modeling in the case of shorter treatment duration of 12 weeks. In this more challenging case, the superiority of modeling approach was clear, with misclassification rates of only 14.6% and 19.2% when using correct or false a priori information, respectively, as compared with 50% when using RVR (Supplementary Material 3 online).

Whereas using M did not much deteriorate the prediction quality, assuming false a priori information on β led to very high error rates (>20%) for all designs, even when BDL data were correctly handled, which is the consequence of a poor estimation of all individual parameters (see Supplementary Material 1 online).

We illustrated the influence of different factors on individual estimation and outcome prediction using examples of five typical patterns of virologic responses found in the simulated population (). As expected, long design allowed obtaining individual fits close to the true, i.e., simulated profiles, especially when a virologic rebound occurred. With lower δ and ε, it is not surprising that M predicted higher treatment duration to achieve SVR, compared to the true model. This explained high rate of FN obtained with this false information (Supplementary Material 2 online). Similarly, a larger number of treatment failures were predicted with higher β, which is associated with more chances of virologic rebound.

Discussion

Here we evaluated the role of a priori information, designs, and methods for handling BDL data in the prediction of individual virologic response and thus the ability of using viral kinetic models to individualize therapy in HCV GT-2/3 infected patients. Unlike what is usually done, we focused on individual parameters and assumed that the population parameters (i.e., the mean and the between-patient variability) were known to estimate individual parameters using Bayesian approach.

Using true a priori population parameters, precise estimation of individual parameters and good prediction of treatment outcome were obtained. Obviously, a major limitation of Bayesian approaches is that a priori values of the population parameters are required. This a priori information could be wrong, for instance if coming from a set of patients with different characteristics. This is why we also evaluated the robustness of this approach when the population parameters were different from the “true” parameters, i.e., those used for data simulation. Our results showed that setting wrong population parameters did not dramatically deteriorate individual parameters estimation as long as the parameters were identifiable (such as δ and ε). Importantly, satisfactory parameter estimation was obtained using a clinically relevant design with a limited number of data points (<6) and a last data point at W4. Therefore, the estimation of viral kinetic parameters with this approach could allow for a rapid prediction of treatment outcome. In the case of β, a parameter that can be correctly estimated at individual level only in patients who have a virologic rebound (11.7% of our simulated population), assuming a higher value led to very large estimation errors in all individual parameters. The combination of these largely biased parameters led to high prediction errors for treatment response (see Supplementary Material 1 online).

BDL data is another factor that could lead to substantial bias. Interestingly, ignoring BDL data did not lead to significant bias as long as a priori population parameters were true. However, a correct handling of undetectable data became critical when population parameters were wrong, and we showed that this can be achieved by likelihood-based approaches that rigorously calculate the information contained in BDL data.

Here we predicted the treatment outcome in a patient using the individual parameter estimates, which are the modes of a posteriori distribution, and therefore defining a yes/no response. In order to apply this approach to clinical practice, one would need to evaluate uncertainty of treatment outcome prediction by using the whole a posteriori distribution of individual parameters, i.e., by incorporating uncertainty on individual parameter estimates. This could be done by reconstructing the a posteriori distribution, for instance using MCMC methods.

The ultimate goal of viral kinetic modeling is to allow for individual parameter estimation and treatment tailoring. To achieve this goal, pharmacometricians will need reliable models and appropriate sampling designs that can provide precise parameter estimation. The information contained in a design can be a priori evaluated and optimized using the Fisher information matrix. This approach became a standard tool for pharmacokinetic/pharmacodynamic studies using NLMEM and was implemented in several software such as PFIM, PopED, POPDES, and POPT.[30] However, viral kinetic models have specific features that are not always properly addressed in these software such as the large proportion of BDL data[31,32] and the need to achieve a good level of precision of both population and individual parameters.[33,34] Clearly, optimal design at population and individual levels may differ, and theoretical studies are needed to provide relevant designs at both levels.

Here we studied the ability of modeling approach in the prediction of IFN-based treatment outcome in HCV GT2/3 patients using a standard viral kinetic model. More sophisticated models are now being developed to explain complex viral kinetic patterns observed with new direct acting antiviral agents.[35] Further studies will be necessary to suggest relevant designs to precisely estimate treatment outcome with this new generation of treatments and models.

Methods

Viral kinetic model. In this study, we used a standard HCV model with liver regeneration[36,37] given by the following set of equations:

The model considers target cells T, infected cells I, and free virus V. Targets cells are produced at rate s, proliferate with rate r, die at rate d and become infected at rate β. Infected hepatocytes die at rate δ and release virions at rate p per cell. Virions are cleared from serum with rate c. Peg-IFN acts by blocking virion production, p, with effectiveness ε varying between 0 and 1. We assumed that the viral infection system before treatment initiation (t = 0) is at steady state. If only viral load data are measured, the parameters s, r, d, T, p cannot be identified individually and therefore were usually set to values found in literature in parameter estimation.[22,25,38] In this study, we investigated only four identifiable parameters β, δ, c and ε, which are also the most clinical relevant parameters. Indeed, in most patients treated with peg-IFN, viral load initially declines in a biphasic manner. This biphasic decline is captured by c, δ and ε. The first phase has a rate c × ε and leads to a viral load decline equal to log10(1−ε), which reflects drug effectiveness. Subsequently, viral load declines with a slower but persistent rate approximately equal to δ × ε. If ε is high enough, viral load will continuously decrease during treatment. Otherwise, a virologic rebound may occur,[37] and the kinetics of viral rebound will depend on the value of β: higher values of β give earlier rebounds.

Statistical model. Let i denote the ith individual (i = 1,…,N) and j the jth measurement of an individual (j = 1,…, n where n is the number of observations for individual i). The statistical model for the viral load y measured in individual i at time t is given by:

where f : nonlinear viral kinetic model, identical for all individuals,

θ: individual parameters vector,

e: residual error supposed to follow normal distribution N(0,).

Because viral kinetics is simulated on long time period, it can happen that viral load takes values very close to 0. In that case, numerical errors may give a negative value for the viral load. This is why, unlike what is usually done, we evaluated the viral load in its natural scale and not in the log10 scale.

The individual parameters θ can be decomposed into fixed effects µ = {β, δ, c, ε} representing mean effects of the population and random effects η specific for each individual. It is assumed that η ~ N(0, Ω) with Ω defined as the variance–covariance matrix in which, each diagonal elements ω representing the variance of the qth component of the random effect vector . We assumed log-distribution for β, δ, c (qth parameter):

and logit-distribution for ε:

The population parameters vector ψ is composed of {µ, Ω, σ }

Bayesian individual parameter estimation. In NLMEM, individual parameters are estimated using Bayesian approach. They are calculated as the Maximum a posteriori (MAP), i.e., the mode of the posterior distribution. Let ψ be the known population parameters vector and let be the conditional distribution of . The MAP estimate of is given by

As µ is known, once is estimated, can be easily calculated.

In presence of data below limit of detection (LOD), the distribution is calculated by:

where is the probability density function, and is the cumulative density function of normal distribution with mean m, variance v, evaluated at x.

Simulation study. This simulation study aimed at mimicking the viral response observed in HCV GT-2/3 infected patients during peg-IFN/RBV. The values for ε, δ, and c in GT-2/3 infected patients are based on literature.[39] For β, this parameter has never been estimated in this population, and therefore we use a value being in the range of possible values and close to those estimated in HCV GT-1 infected patients.[36,40] We used population parameters presented in the M column () for data simulation.

Four nested designs D ⊂ D ⊂ D ⊂ D were considered. The richest design, D, contains 12 sampling times until the end of the treatment (days 0, 1, 4, 7, W2, 4, 6, 8, 12, 16, 20, 24). D (days 0, 1, 4, 7, W2, 4) does not have any data after W4, D (days 0, 7, W2, 4) does not have any data after W4 and during the first week and D has only one measurement at baseline and W4 (day 0, W4).

We simulated N = 1,000 vectors of random effects to calculate individual parameters θ from the population parameters of the true model. We then predicted the viral kinetics in these 1,000 in silico patients using the mathematical model (Eq.1). Residual errors were then simulated and added to viral load data for all sampling times of D to obtain dataset for D design. As the four designs are nested, to generate datasets for other designs, we selected the simulated viral loads corresponding to the sampling times of the design of interest (Supplementary Material 4 online). As the same data of the same 1,000 simulated individuals were used to generate all the datasets, all the differences obtained in results were not due to Monte Carlo errors.

We considered a limit of detection LOD = 45 IU/ml. Since successful therapy leads to SVR, i.e., the eradication of HCV, we hypothesized the existence of a threshold, called “cure boundary,” under which infection is considered eradicated. Following what has been done previously, we assumed that this threshold corresponds to having <1 infected cell in the whole extracellular body fluid, assumed to be 15 L.[22] Although this cure boundary has not yet proven biologically correct, prediction of treatment outcome using this threshold has been shown to match SVR rate observed during peg-IFN based therapy.[22] If the course of infected cells crossed the cure boundary, we considered that SVR was achieved and no viral relapse could occur afterwards. All viral loads obtained after achieving the cure boundary are considered undetectable.

To evaluate the impact of a priori information on individual parameter estimation and outcome prediction, we modified fixed effects of some parameters to generate two false models (): for M, we set δ and ε to the values obtained in GT-1 infected patients; for M we set β to the upper limit of possible values found in literature.[40]

To estimate individual parameters for 1,000 simulated patients, individual random effects η were first estimated and then used to calculate individual parameters. We supposed that population parameters were already known and fixed them at different sets of a priori values corresponding to M, M, M to perform MAP. BDL data were handled using two approaches: taking into account BDL data in the likelihood function (likelihood-based approach or LBA) or omitting all BDL observations in the dataset (OMIT).

To evaluate the ability of estimation of the component q of an individual parameter, we used the following criteria: Mean Relative Error (MRE) and Relative Root Mean Square Error (RRMSE) defined by:t

For ε, we evaluate the MRE and RRMSE of log10(1−ε).

We also used shrinkage (Shr) to evaluate the quality of individual parameter estimates (Eq. 9). A low shrinkage (<50%) indicates that individual parameter is correctly estimated.

Once estimated, individual parameters were then used to simulate the individual predicted profiles of target cells, infected cells, and viral loads assuming a 24- or 12-week treatment. If infected cells were predicted to achieve the cure boundary during treatment, the predicted treatment outcome is having SVR and vice versa.

By comparing the predicted treatment outcome with the true (simulated) response, we calculated the percentage of false-negative FN (patients having SVR predicted to remain infected), false-positive FP (uncured patients predicted to have SVR) and of misclassification MC (sum of FN and FP). Using these criteria, we compared the quality of treatment outcome prediction between different scenarios and with the case where no individual observation is available. To evaluate the gain of using modeling approach compared to empirical rules used in clinics, we also evaluated the predictive ability for SVR of “rapid virologic response,” RVR—defined as undetectable viral load at W4. We also considered a shorter treatment that was stopped at W12 in order to evaluate the ability of outcome prediction in a challenging case.

All simulation was done using R 2.14.0. Parameter estimation (Supplementary Material 5 online) was performed using MONOLIX 4.1.2.

Conflict of Interest

T.H.T.N. has a research grant from Novartis. J.Y. and M.L. are working for Novartis, in Modeling and Simulation department. As an associate editor for CPT: PSP, F.M. was not involved in the review or decision process for this paper.

Study Highlights

Table 1

Mean relative bias, RRMSE, and shrinkage (in %) of β, δ, c, ε for different scenariosa

Table 2

Parameter values of the different viral kinetic models in response to treatment (M, M, M)