This study is based on the follow-up cohorts of PLWH at Nanfang Hospital and the Fifth Hospital of Zunyi, where participants have been undergoing long-term treatment and regular follow-ups at HIV clinics. A total of 1577 participants were enrolled based on defined inclusion and exclusion criteria. The inclusion criteria were: 1) a baseline CD4+ T cell counts of less than 350 cells/μL at the initiation of HAART, with continuous follow-up for 2 years, and two HIV RNA measurements of less than 50 copies/mL; 2) age 18 years or older, with complete baseline, 6-month, 12-month, and 24-month CD4+ T cell counts. The exclusion criteria included: 1) poor treatment adherence or a history of treatment interruption; 2) concurrent malignancy or long-term use of immunosuppressive medications; and 3) incomplete clinical data. As illustrated in Figure 1 , the cohort from Nanfang Hospital was divided into a training set and an internal validation set in a 7:3 ratio, while the cohort from the Fifth Hospital of Zunyi was designated as the external validation set.

The research received approval from the Institutional Ethics Committee of Nanfang Hospital (study identifier: NFEC-2021-448) and adhered to the Helsinki Declaration of 1964, along with its subsequent updates. Informed consent was obtained from all participants.

We systematically collected demographic and clinical parameters of participants including age, gender, HAART regimens, HBsAg positivity, anti-HCV positivity, HIV viral load, and laboratory measurements at baseline, 6 months, 12 months, and 24 months into treatment. These measurements encompassed CD4+ T cell counts, CD8+ T cell counts, CD4/CD8 ratios, Platelet (PLT), creatinine (CR), hemoglobin (HGB), white blood cell count (WBC), aspartate aminotransferase (AST), alanine aminotransferase (ALT), triglycerides (TG), total cholesterol (CHOL), and fasting plasma glucose (FPG). The aforementioned data were obtained from clinical records or databases.

Currently, there is no universally accepted definition for immune reconstitution failure. In this study, INR was defined as having two consecutive HIV RNA measurements <50 copies/mL after two years of HAART, still maintaining a CD4+ T lymphocyte count of <350 cells/µL (Cuzin et al., 2007; Gunda et al., 2017).

In this study, variables from the training set that demonstrated significance at a p-value <0.05 in univariate analysis were included in the model construction. We employed several machine learning algorithms to predict INR classification, including the Logistic Regression Model (LRM), Random Forest (RF), XGBoost, Support Vector Machine (SVM), Naive Bayes, Decision Trees, neural network, and k-nearest Neighbors (KNN). To prevent overfitting and enhance the generalizability of the models, a 10-fold cross-validation method was employed for model evaluation, with iterative refinements through repeated trials.

To further assess and compare the predictive performance of these models, we constructed receiver operating characteristic (ROC) curves and determined the area under the ROC curve (AUC). An AUC value closer to 1 indicates better predictive performance. Additionally, we utilized calibration curves to evaluate the consistency between the observed and predicted risks. The more the calibration curve of the model aligns with the 45 - degree line, and the closer the value of the Brier score is to 0, the more the predicted probability matches the observed event incidence. Furthermore, decision curve analysis (DCA) was used to evaluate the clinical utility of the models. By comparing the net benefits of the model with two default strategies (treating all or none), DCA provides insights into the clinical value of the models.

To improve the interpretability of machine learning models, which are often regarded as “black box” models due to their complex and opaque decision-making processes, we applied Shapley Additive Explanations (SHAP) analysis. SHAP is a cooperative game theory-based approach that quantifies each feature’s contribution by assessing its influence on model predictions. A SHAP value greater than 0 indicates a positive contribution of the feature to the prediction, while a value less than 0 indicates a negative contribution. The larger the SHAP value, the greater the feature’s influence on the prediction. In our study, we visualized these contributions using importance ranking charts, which highlight the relative weight of each feature in influencing the outcome. Additionally, we employed partial dependence plots to demonstrate how each feature affects the predicted results, illustrating the relationship between individual features and the model’s output while considering the influence of other variables.

In our analysis, datasets that conformed to a normal distribution were described using the mean ± standard deviation, and comparisons between two groups were conducted using Student’s t-test. For datasets that were non-normally distributed, comparisons were made based on the median and interquartile range, with the Mann-Whitney U test applied for statistical evaluation. Categorical variables were summarized as frequencies and percentages and analyzed using either the chi-square test or Fisher’s exact test, as appropriate. Independent risk factors for INR were identified through univariate and multivariate logistic regression analysis. To evaluate the dose-response relationship between continuous variables and INR, we employed restricted cubic splines (RCS). This method enables the visualization and quantification of potential non-linear associations, and by analyzing the shape of the dose-response curve, we can identify critical thresholds where the relationship between the predictor and the outcome changes. It is important to note that all aspects of data analysis and graphical representation were performed using R version 4.2.1. All tests conducted in this study were two-tailed, and a p-value <0.05 was considered statistically significant.

In the longitudinal cohort study of PLWH to predict the risk of INR during follow-up, we retrospectively included 903 PLWH from Nanfang Hospital and 674 PLWH from the Fifth Hospital of Zunyi University, who had been under treatment for more than two years. These cohorts served as the internal and external datasets, respectively. As shown in Table 1 , Nanfang Hospital enrolled 903 participants, with 532 achieving immune response (IR) and 371 not achieving IR, while the Fifth Hospital of Zunyi University included 674 participants, with 408 in the IR group and 266 in the INR group. In both cohorts, the INR group exhibited significantly higher ages and viral loads compared to the IR group, while CD4+ T cell counts were notably lower in the INR group. There were no significant differences between the two groups in terms of gender, HAART regimens, and the prevalence of baseline HBsAg and anti-HCV.

We visualized the clinical characteristics of PLWH at each follow-up point using line graphs ( Figure 2 ) and compared the levels between the IR group and the INR group. We observed that at each follow-up point, the IR group exhibited higher levels of CD4+ T cells, CD4/CD8 ratio, WBC counts, HGB levels, and PLT levels compared to the INR group. However, differences in CD8+ T cells, liver function markers such as ALT and AST, lipid levels including TG and CHOL, renal function as indicated by CR, and FPG were only present at certain follow-up points. A similar analysis was conducted in the external dataset ( Supplementary Figure 1 ), and the results were consistent. The only exception was that the CD8+ T cell levels were also higher in the IR group compared to the INR group.

To investigate the factors influencing INR, we conducted a univariate logistic analysis that identified 20 significant variables ( Figure 3 ). Given the potential for multicollinearity among these variables, we conducted a collinearity test on variables with a p-value < 0.05 from the logistic univariate analysis by calculating the variance inflation factor (VIF) ( Supplementary Figure 2 ). Since all parameters had a VIF value less than 10, all were included in the multivariate analysis and identified independent factors for INR as Baseline-CD4 (OR = 0.995, P = 0.030), 6M-CD4 (OR = 0.992, P < 0.001), 12M-CD4 (OR = 0.993, P < 0.001), Baseline-HGB (OR = 1.023, P = 0.002), and 6M-HGB (OR = 0.968, P = 0.014).

To further analyze the relationship between baseline parameters and INR, we conducted the same analysis and found that in the multivariate analysis ( Supplementary Figure 3 ), age (OR = 1.021, P = 0.010), HIV load (OR = 0.725, P = 0.009), baseline CD4 (OR = 0.983, P < 0.001), baseline WBC (OR = 0.842, P = 0.008) and baseline HGB (OR = 1.012, P = 0.014) were independently associated with INR.

Through RCS analysis, we further investigated the relationship between independent factors and INR incidence ( Figure 4 ). We observed that 6M-CD4 and 6M-HGB showed a linear relationship with INR (overall p<0.05, nonlinearity p>0.05), with threshold concentrations of 273 cells/μL and 127.47 g/L, respectively. Conversely, a nonlinear relationship was evident between Baseline-CD4, 12M-CD4, Baseline-HGB, and INR (overall p<0.05, nonlinearity p<0.05). The risk of INR rapidly increased when Baseline-CD4 was below 165 cells/μl, 12M-CD4 was below 293 cells/μl, and Baseline-HGB was less than 125.23 g/L.

We divided the internal dataset into a training set for model construction and an internal validation set following a 7:3 split, while the external dataset served as the models’ external validation set. We compared the baseline clinical characteristics across the three datasets ( Supplementary Table 1 ). The median age of PLWH in all three datasets was 32 years old, and the proportion of INR was similar across the datasets. Notably, the external validation set had a higher proportion of female PLWH and a lower proportion using INSTI-based treatment regimens.

Subsequently, we incorporated significant variables in the univariate analysis in Figure 3 into model construction, including Baseline and 6-month/12-month CD4+ T cells, CD4/CD8 ratio, WBC, HGB, PLT and etc. Using these variables, we developed eight predictive models employing machine learning methods. We then validated the stability and generalizability of these eight models across the training, internal, and external validation sets. Ultimately, the RF model exhibited the best clinical predictive performance across all datasets, with AUROC values of 0.866, 0.943, and 0.897, respectively ( Figures 5A–C ). In terms of calibration, the RF model outperformed other models in all three datasets, with Brier scores of 0.136, 0.102, and 0.126 ( Figures 5D–F ). In clinical utility assessment, the DCA curves of the RF model were consistently higher than the “treat all” and most other model lines across the majority of threshold probabilities, indicating significant clinical application value ( Figures 5G–I ).

Given the RF model’s outstanding predictive capability across both internal and external validation datasets, we ultimately designated it as the best-performing model. To clarify the clinical relevance of specific features, this research quantified their importance using SHAP values. The variables were prioritized by their impact on predicting INR risk ( Figure 6A ), identifying the top five predictors in PLWH after two years of HAART as 6-month CD4+ T cells, 12-month CD4+ T cells, baseline CD4+ T cells, 6-month CD4/CD8 ratio, and 12-month CD4/CD8 ratio. Consequently, CD4+ T cell counts measured between 6 and 12 months post-treatment are critical for assessing immune reconstitution.

Through the summary plot ( Figure 6B ), we detailed the positive and negative relationships between features and outcomes, finding that higher CD4+T cell counts were associated with a lower probability of INR, and older age correlated with a higher probability of INR. Subsequently, we illustrated the impact of model variables in predictions for an example of PLWH with IR and INR respectively ( Figures 6C, D ). Finally, we generated a partial dependence plot ( Figure 6E ). Specifically, the critical threshold for CD4+ T cell counts was observed around 350 cells/µL at 12 months, 250 cells/µL at 6 months, and 150 cells/µL at baseline. For the 6-month CD4/CD8 ratio, maintaining a value near 0.5 was associated with minimizing INR risk. When the parameter values fall below these critical thresholds, the risk of INR increases. Nevertheless, it is noteworthy that the partial dependence analysis did not detect significant correlations between variables and age.