A Prognostic Model for Estimating the Time to Virologic Failure in HIV-1 Infected Patients Undergoing a New Combination Antiretroviral Therapy Regimen

Modern combination antiretroviral therapy (cART) can suppress plasma viral load to undetectable levels in a large proportion of HIV-1 infected patients. The risk for a patient to experience virologic failure has been decreasing consistently during the last decade in high-income countries [1–5].

Part of this high rate of virologic success may be due to the increasing ability of more potent and tolerable antiretroviral drugs targeting a wider range of molecular targets, providing physicians the opportunity to tailor cART regimens according to patients' background and to manage treatment failure promptly. Additionally, better understanding of the mechanisms of drug resistance allows optimization of treatment regimens based on an individual's virus genotype. Although the prevalence of drug resistance seems to be decreasing in recent years [6], despite modern cART options [7], drug resistance remains a concern in chronically infected patients with a long treatment history and in treatment-naïve patients who have been infected with drug resistant isolates [8]. Overall, patients remain at risk of developing drug resistance.

Genotypic susceptibility scores (GSSs) have been developed for interpreting HIV-1 drug susceptibility based on the sequence of the virus genome coding for drug targets. GSSs usually consist of a set of general rules for scoring susceptibility to individual drugs. These systems are curated and updated by panels of experts and made freely available via the internet, and in some cases, as web-services. The most popular systems include the Stanford HIVdb [9], Agence Nationale de Recherche sur le Sida (ANRS) [10], and Rega [11] algorithms. Recently, machine learning methods have been also introduced to model HIV-1 drug resistance and to optimize cART design. These machine learning methods explore a larger set of variables besides the viral genotype (such as viral load, CD4+ T cell counts, demographic information, treatment history), and include techniques such as artificial neural networks [12], mutagenetic trees [13], Bayesian networks [14], and random forests [15]. Both GSSs and machine learning approaches have been proven to usefully predict virologic outcome at fixed time points (e.g. n-weeks) after cART initiation or switch [14–18]. The state-of-the art systems, available as free web-services, are able to select a set of suitable cARTs for a patient, given the patient's viral genotype and other background information, ensuring the maximal probability of viral load reduction after n-weeks of uninterrupted treatment [19, 20].

Systems that predict the actual time to virologic failure, indicated by a viral load rebound following a cART-induced viral load reduction, or no viral load reduction after a few months of uninterrupted therapy, have not been developed yet. However, GSSs have been already associated with the time for achieving an undetectable viral load [18]. In this work we explore the predictive ability of linear and non-linear survival models with respect to the time to virologic failure endpoint, along with its prognostic factors. A model that could predict an individual's duration of virologic success with a cART regimen would provide valuable information in tailoring cART choice.

We observed 1,067 virologic failures over 2,820 person years of follow-up (rate = 37.8 per 100 person years). By Kaplan-Meier estimation, median (95% CI) time to virologic failure was 659 (533-784) days for all patients, and 2,510 (1,715-N/A) days for those previously therapy-naïve. By two years, the estimated proportion of patients not experiencing virologic failure was 0.48 (0.46-0.51) when considering the whole set of patients, and 0.71 (0.67-0.75) for those previously therapy-naïve.

The evaluation of extra-sample predictive performance (via the bootstrap .632 method, on 100 runs) is shown in Figure 1. In detail, the multivariable Cox models fitted with different GSSs yielded an average (st.dev) c-index of 0.7060 (0.007) for Rega, 0.7048 (0.007) for ANRS, and 0.7068 (0.007) for HIVdb. As expected, univariable Cox models fitted with the single GSSs were greatly outperformed by their multivariable versions, yielding -respectively- an average (st.dev) c-index of 0.6277 (0.007), 0.6271 (0.007), and 0.6330 (0.007). Additionally, a likelihood-ratio-test conducted on the Cox regression, considering a null model with the sole GSS and an extended model with all covariates, reported a consistently worse fit of the null model as compared to the extended model (L_null = -7573.536, L_extended = -7361.47, χ² = 424.13 on 41 degrees of freedom, p < 0.0001).

The c-index performance of RSF, grown with the limited number of 30 trees due to the high computational burden, using the same covariate settings as in the multivariable Cox regression, were 0.7298 (0.009) for Rega, 0.7276 (0.009) for ANRS, and 0.7319 (0.008) for HIVdb.

By executing all pairwise Student's t-tests (adjusted for multiple comparisons with the Benjamini-Hochberg method) comparing the c-index distributions, we find statistically significant differences between univariable vs. multivariable Cox models, univariable Cox models vs. RSF, and multivariable Cox models vs. RSF (all p < 0.0001). When comparing the performance of Rega, ANRS and HIVdb within the same model settings, we found evidence of a significantly better performance of HIVdb as compared to ANRS in the univariable Cox (p < 0.0001), multivariable Cox and (p = 0.053), and RSF (p = 0.0005). HIVdb outperformed Rega only under the univariable Cox modeling. Conversely, Rega and ANRS did not show any appreciable difference in any of the within-model c-index distributions. However, it has to be noted that all the GSS had the same average c-index values up to the third decimal.

When executing a sensitivity analysis on the subset of treatment-naïve patients (n = 733), multivariable Cox regression confirmed the relative hazards of the GSS of the regimen (RH = 0.50, 95%CI 0.36-0.70, p < 0.0001, per one point increase of HIVdb), and of HIV-1 RNA (RH = 1.51, 95%CI 1.15-1.97, p = 0.0026, per one Log10 copies/ml higher). Other factors independently associated with the endpoint were being non-Italian born (RH = 1.87, 95%CI 1.09-3.22, p = 0.023, as compared to Italian born), subtype C (RH = 2.16, 95%CI 1.05-4.44, p = 0.036, as compared to subtype B), and calendar year 2004 (RH = 2.08, 95%CI 1.31-3.31, p = 0.0019, as compared to 2007 and after). Differences among cART regimens were not relevant.

In a second sensitivity analysis, only patients who reached an HIV-1 RNA < = 50 copies/ml from the cART start date onwards were selected (n = 1,578, of which 622 treatment-naïve). By fitting a multivariable Cox regression, variables independently associated with the endpoint were consistent with those obtained for the subset of treatment-naïve patients (data not shown).

In this study we investigated linear and non-linear survival models for predicting the time to virologic failure in HIV-1-infected patients undergoing a new cART regimen, with the aim to assess both prognostic factors of virologic failure and performance of predictions, in light of the development of a reliable expert system.

In contrast to 2NRTI+1PI/r, older age, higher HIV-1 RNA, and lower CD4+ counts, an increased hazard of virologic failure was associated with low GSSs of cART, a less recent calendar year, administration of 2NRTI+1NNRTI, and the ART-naïve status. HIV-1 RNA and GSS remained associated with the endpoint when considering only treatment-naïve patients in a sensitivity analysis.

When looking at the goodness-of-fit of the models, the inclusion of additional covariates besides the GSS in a multivariable Cox regression yielded a significant improvement in the likelihood. Furthermore, RSF proved to be a promising approach, improving performance over that obtained by using the Cox method.

This study has some limitations. First, there was a study selection bias, since only cART regimens that were undertaken for at least 90 days were considered. Early-switches and simplifications were excluded from the analysis. In addition, the main study endpoint was a pure virologic criterion and did not include stops due to other reasons. Although we selected only clinics with the ability to perform ultra-sensitive assays, the viral load threshold at > 400 copies/ml was arbitrary and might not capture all the actual virologic failures. Other thresholds for defining virologic failure, such as a viral load > 50 copies/ml or > 1,000 copies/ml might overestimate or underestimate true failure, respectively. We observed a lower hazard of virologic failure for regimens containing ritonavir-boosted PI as compared to those containing NNRTI, and this might be a reflection of both the selection bias and the endpoint definition.

In regards to the statistical methods, we did not investigate the potential benefit in terms of likelihood fit given by the inclusion of interaction terms in the Cox regression. It is possible that a Cox model with higher-order interactions is able to reach the performance obtained by using the RSF. For the practical perspective of a time to event prediction model, Cox regression presents some problems with the baseline hazard function estimation, while the RSF gives output in terms of mortality ensembles. We also tested accelerated failure time models, which are able to give reliable predictions in terms of actual time scales, and their performance was comparable to the results of the Cox regression (data not shown).

This study demonstrates the feasibility of an expert system that predicts the actual time course of a new cART as the estimate of time to virologic failure, and the reliability of the system's predictions can be improved both by including additional covariates besides the viral genotype and by using non-linear regression techniques. The implementation of such a system would create a more clinically-oriented treatment decision tool and more effectively tailor patient cART regimens.