Introduction

The development of ART has vastly improved morbidity and mortality associated with HIV infection. However, a long-lived proviral reservoir precludes a cure of infection (1–4). Furthermore, despite ART, PWH are still at increased risk of morbidity and mortality related to sequelae of persistent immune activation. This persistent immune activation is thought to be multifactorial: a product of residual immune pertubation from untreated HIV infection, persistent proviral expression during ART, and specific clinical risk factors for immune activation that are enriched amongst PWH. Understanding the interplay between the HIV reservoir and persistent immune activation during long-term ART is important to inform strategies for a cure.

While a number of cell types have been proposed as contributing to the HIV reservoir, the largest and most well-characterized reservoir is within CD4 T cells. Molecular analyses of residual HIV proviruses in CD4 T cells during ART have shown that the reservoir is comprised of a heterogeneous mixture of full-length replication-competent proviruses and a larger set of defective proviruses containing internal deletions or APOBEC-mediated hypermutation(s) (5–7). Infected cells can undergo clonal expansion and contraction over time, potentially in response to homeostatic cues or antigen-driven activation (8–10). Studies of SIV-infected non-human primates and HIV-infected humans show that a pool of latently-infected cells forms very early during infection (3,11). However, studies suggest the reservoir is enriched in viral variants that were circulating at the time of ART initiation, indicating that initiation of ART may promote the entry of recently infected cells into the reservoir (12, 13).

HIV infected cells are rare in peripheral blood, making measuring and characterizing the reservoir challenging. Additionally, the size of the reservoir is highly variable between individuals, ranging over multiple logs in scale. During therapy, a subset of the reservoir remains transcriptionally active, expressing viral RNA and, in some cases, viral protein, leading to persistent viremia or “blips” during therapy (14). The factors that regulate reservoir size and transcriptional activity are only beginning to become clear. Similarly, the impact of the latent reservoir on the immune system of PWH is not fully understood. Over time on therapy, intact proviruses become preferentially depleted from transcriptionally active genomic regions and enriched in heterochromatic regions, suggesting that ongoing viral expression promotes infected cell depletion, either through a cytopathic effect from viral proteins or from immune killing of infected cells with residual expression (15).

In addition to the presence of a latent HIV reservoir, PWH on therapy exhibit a persistent level of T cell immune activation and other altered phenotypes compared to uninfected individuals, involving cells of both the innate and adaptive immune systems. Specifically, HIV infection is associated with increased levels of inflammatory biomarkers and increased levels of activation for immune cells (16–19). PWH also exhibit a depleted naïve T cell compartment (20) and increased expression of activation/exhaustion and senescence markers such as PD-1 and KLRG1 (21–24). These alterations may contribute to increased levels of disease and death in PWH. The molecular mechanisms underlying persistent immune activation and dysfunction for PWH on ART are unclear, but it is possible that ongoing viral gene expression and reverse transcription generates ligands for pattern recognition receptors in the innate immune system (25). Viral protein expression during therapy likely maintains antigen-specific cells of the adaptive immune system. Given that viral protein expression may be relatively rare during therapy, another important potential cause of persistent immune activation during ART is persistent sequelae of the massive disruptions to T cell homeostasis and lymphoid structures that occur during untreated HIV infection.

It is possible that the HIV reservoir and the host immune system interact in multiple complex ways during therapy. The dynamics of how the HIV reservoir affects the immune system and vice versa are important to characterize and may lead to novel ideas to promote reservoir depletion or to inhibit immune dysfunction that may be related to ongoing HIV expression in persistently infected cells. A number of previous studies have examined the relationship between reservoir size and the immune system (26–45). Consistent major findings from these studies include enrichment of the HIV reservoir in CD4 T cells expressing activation and exhaustion markers, positive correlation of the reservoir with central memory T cells (Tcm) and transitional memory (Ttm) T cell subsets, and inverse correlation of the HIV reservoir with the CD4 nadir and the CD4/CD8 T cell ratio. One important limitation of most of these studies is that they were not able to distinguish intact proviruses from defective proviruses, making the relationship between the immune system and the intact reservoir unclear. Recent studies indicate enrichment of the intact HIV reservoir in cells expressing VLA-4, but no correlation with classic T cell activation markers (26, 37). Additional assessment of immune correlates of the intact HIV reservoir in additional cohorts is needed. Furthermore, most of these studies examined a limited number of surface markers and may have missed associations that would be revealed by a higher-resolution examination of immunophenotypes. Finally, there is a recent, growing appreciation for the profound influence of ART initiation timing and length of treatment on the size and composition of the HIV reservoir that must be accounted for in assessing persistent immune activation in PWH on therapy (46, 47).

In this study, we use a mixture of traditional statistical approaches and novel data science methods to assess the relationship between the immune system and the HIV reservoir across a cohort of 115 PWH. Collectively, this study highlights the utility of machine learning approaches to identify otherwise imperceptible global patterns in high-parameter studies of the HIV reservoir. Additionally, these approaches corroborate recent findings in the HIV persistence field regarding selective intact proviral decay and immune dynamics and highlight the complex immunologic signatures of HIV latency.

Results

Cohort description, reservoir quantification, and immunophenotyping

115 participants with HIV infection for at least one year and who had been receiving antiretroviral therapy (ART) for at least 0.9 years (median 9 years, IQR 5.2-16.6) were studied (Table 1). The study cohort was 77 percent male and median participant age was 45.

Table 1:

The Intact Proviral DNA assay (IPDA) was performed on isolated total CD4 T cells (48). The IPDA enables estimation of the frequency of total proviruses, intact proviruses, and the percent intact of total proviruses for each participant. Immunophenotyping was performed using 25-color spectral flow cytometry on peripheral blood mononuclear cells (PBMCs) from each study participant. A representative flow gating is shown in Figure S1. Briefly, live lymphocytes were first defined as T cells (CD3+/CD56-), NK cells (CD3-/CD56+), or NKT cells (CD3+/CD56+). T cells were then further subdivided as CD4+ or CD8+ T cells, and then as naïve T cells (Tn, CD45RA+/CCR7+), central memory T cells (Tcm, CD45RA-/CCR7+), effector memory cells (Tem, CD45RA-/CCR7-) and terminally differentiated effector cells (Teff, CD45RA+/CCR7-). Expression of various surface markers within each of these subsets was then examined (CD38, HLA-DR, PD-1, KLRG1, CD127, CD27, NKG2A). In addition to measuring the baseline abundance of surface proteins, we also stimulated PBMCs from each participant with a pool of HIV-derived peptides and measured intracellular cytokines TNFα, IL-2, and IFNγ as well as the surface expression of the degranulation marker CD107a following background subtraction of the signal from unstimulated cells. Frequencies of CD4 or CD8 T cells that were single- or double-positive for TNFα, IL-2, IFNγ, or CD107 were used for analysis.

These analyses resulted in the determination of a total of 144 parameters for the cohort of n = 115 PWH. These 144 parameters include 133 immunophenotypic cell frequencies, 3 HIV reservoir parameters (frequency of intact HIV DNA, frequency of total HIV DNA, and percentage intact of total HIV DNA), and 8 clinical-demographic parameters. Clinical-demographic variables (Table 1) included age, biological sex, years of ART, estimated years of HIV infection prior to ART (here, we use categorical variables: years before ART=NA, years before ART< 1, years before ART≥ 1), CD4 nadir, most recent CD4 T cell count, and race (Caucasian, African-American, other). The median and interquartile range (IQR) for all immune cells population variables are presented in Table S2. No significant differences in the reservoir based on biological sex or race were observed (Table S1).

Correlation of total, intact, and percent intact HIV DNA with immunophenotypes

For the initial analysis, we performed Spearman correlations of the clinical-demographic and immunophenotypic data with HIV reservoir metrics. 69 host variables correlated (unadjusted p < 0.05) with one or more characteristics of the HIV reservoir (Table S3). Following Benjamini-Hochberg p-value adjustment for multiple comparisons using a false discovery rate correction at 5%, 31 variables were found to be correlated with total HIV DNA, 3 with intact HIV DNA, and 4 with the percentage of intact HIV DNA (Table 2). Overall, this analysis indicated that the HIV reservoir size is significantly associated with the abundance of several immunophenotype frequencies and clinical-demographic variables. Of note, the correlations observed were generally weak to moderate in magnitude (the mean of the absolute values of Spearman correlation coefficients was 0.31 for total HIV DNA, 0.33 for intact HIV DNA, and 0.35 for percent intact HIV DNA).

Table 2:

With regard to clinical-demographic variables, CD4 nadir (Spearman r = −0.32) was negatively correlated with total HIV DNA, whereas age (r = 0.31) and years of ART (r = 0.31) were positively correlated with total reservoir frequencies (Table 2). In contrast, years of ART (r = −0.45) was negatively correlated with the percentage of intact proviruses (Table 2). Given the growing appreciation in the field for the role of ART initiation and timing on the HIV reservoir and immune recovery, we examined this association further (Figure 1). Correlation plots of total, intact, and percentage intact HIV DNA versus years of ART were visualized (Figure 1A-C). As expected, the plots showed no significant correlation for intact HIV DNA versus years of ART (Figure 1B), while total reservoir size was positively correlated with the time of ART (Figure 1A, Spearman r = 0.31). Interestingly, when we examined a plot of percent intact proviruses versus time on therapy (Figure 1C), we observed a biphasic decay pattern. Indeed, when we fitted piece-wise linear functions with allowances for up to 2 breaks, the best-fit model had decay slope changes at 1.6 and 9 years of ART (R² score of piece-wise linear model with two breaks was 44.18%, while linear model without any breaks – 16.3%). In PWH within approximately 0-6 years of ART initiation, a significant proportion of the participants exhibited high fractions of intact proviruses (50-90%). After approximately 6 years of ART, however, the majority of participants had reservoirs for which intact proviruses represented a minor fraction of the overall reservoir. Of note, several immune features also demonstrated statistically significant correlations with years of ART, including CD8 T cell expression of CD38, CD8 Tcm frequencies, and CD107a+IL-2+IFNγ-TNFα-CD8 T cells (Figures S2, S3, S4, Table S4).

Figure 1:

Examination of immunophenotypic parameter correlations with HIV DNA metrics demonstrated several notable correlations (Table 2, Figures S2, S3, S4). As previously reported, CD4 T cell frequencies were negatively correlated with total HIV DNA frequencies (r = −0.36), whereas CD8 T cells were positively correlated (r = 0.41). This relationship held true for intact HIV DNA frequencies but not for the percentage of intact HIV DNA (Table 2). Frequencies of CD4 and CD8 T naïve cells, as well as their expression of CD38, were significantly negatively correlated with total HIV DNA frequency. By contrast, CD4 and CD8 Tcm and Tem cells were weakly positively correlated with total HIV DNA. The direction of these correlations was concordant for intact HIV DNA frequencies, albeit without statistical significance. Markers of T cell activation and exhaustion were weakly but significantly positively correlated with total HIV DNA, including NKG2A+ CD4 T, PD-1+ CD4 T, HLA-DR+ CD4 T, PD-1+ CD8 T, and PD-1+/CCR7+CD8 T frequencies, among others (Table 2, Figures S2, S3, S4). Again, the direction of these correlations were conserved for intact HIV DNA frequencies, albeit without statistical significance. For both intact and total HIV DNA, an inverse correlation with the CD4/CD8 T cell ratio was observed (Figure S5).

In general, the magnitude of correlation was stronger for total HIV DNA compared to intact HIV DNA (Table 2). This pattern may reflect a larger influence of total HIV DNA on immunophenotype of the host. A notable exception to this pattern was CD127+ CD4 T cell frequencies, which were significantly negatively correlated with intact but not total HIV DNA. For percentage of intact HIV DNA, negative correlations were observed with the frequencies of 2 subsets of bifunctional HIV peptide-stimulated CD4 T cells (IL-2+TNFα+ and IFNγ+IL-2+), as well as years on ART. Changes in immune markers in this study relative to years of ART were also indicative of progressive immune recovery (Figure S2-S4 and Table S4). Specifically, years of ART negatively correlated with CD38 expression on CD8 T cells and positively correlated with the frequency of CD8 Tcm cells as well as CD107+IL-2+IFNγ-TNFα-HIV-specific CD8 T cells (Table S4). There are likely many more immune parameters that change in concert with ART as recently reported (47); however, the noise associated with cross-sectional immunophenotype data likely precludes their identification in this study.

Leave-One-Covariate-Out (LOCO) inference identifies specific HIV reservoir and clinical-demographic parameters important for the prediction of several immunophenotypes

Given the evidence that reservoir and immune recovery dynamics likely occur in concert, we next used a variable importance approach - Leave-One-Covariate-Out (LOCO) inference analysis, to account for potential confounding variables and more carefully assess the relative importance of a given variable to an immunophenotype. The LOCO inference approach is demonstrated in Figure 2A. First, a least-squares regression linear model that predicts the dependent variable is fitted on a set of independent variables, and an R² value is generated (Figure 2A Step 1). Next, one independent variable is excluded from the model, the model is refitted, and an adjusted R² value is determined (Figure 2A Step 2). Finally, the change in R² scores, ΔR², is calculated (Figure 2A Step 3). This process is then repeated for all independent variables in the model, and the ΔR² values are plotted in a heatmap for the variables of interest (Figure 2B-D).

Figure 2:

We applied the LOCO inference analysis to all 133 immunophenotype parameters (dependent variables). Three separate analyses were conducted for a set of independent variables consisting of age, sex, years of ART, CD4 nadir, recent CD4 account, years of HIV infection prior to ART (=NA, < 1, ≥ 1), and 1 of 3 HIV reservoir metrics (total HIV DNA, intact HIV DNA, or percentage intact HIV DNA).

Multicollinearity analysis confirmed that independent variables from this analysis were not correlated to a degree that would interfere with model generation (Table S5). Each analysis produced a heatmap of 133 dependent immunophenotype variables, with visualization of the variables for each immunophenotype (Figure S6, Tables S6-S8). For each of these analyses, we displayed the top 10 variables with the largest ΔR² value for LOCO inference incorporating either total (Figure 2B), intact (Figure 2C), or percentage intact (Figure 2D) HIV DNA. In general, the least-squares linear regression models were at best weakly explanatory (R² values in range [-8.38%, 31.69%]) of the global variability in immunophenotypes (Tables S6-S8, Figure 2B-D, Figure S7), however, the LOCO analysis is still valuable for understanding which clinical/demographic variables and HIV characteristics are important to the immunophenotypes.

Across all 3 analyses, age emerged as an important covariate of several immunophenotypes related to T cell subsets, particularly the naïve and central memory compartments, as well as KLRG-1, a marker of immune senescence (Figure 2 and Figure S6). This is consistent with well-described age-related changes in T cell subset populations (49, 50). Notably, biological sex was a key variable for model prediction of NK cell frequency and NK cell NKG2A expression. Years of ART was important for predicting the frequencies of HLA-DR+ CD8 T cell subsets including Tcm, Tn, and Tem, but more so for the model incorporating intact, rather than total, HIV DNA (Figure 2B and 2C). This suggests that both years on ART and total HIV DNA are important explanatory variables for the variance in CD8 T cell activation levels. Total HIV DNA was also noted to be an important predictor of PD-1+ CD4 T naïve cells, HIV-specific IL-2+TNFα-IFNγ-CD107-CD4 T cells, and CD4/CD8 T cell frequencies, among others (Figure 2B, Figure S6). Total HIV DNA and, to a lesser extent, intact HIV DNA were important predictors of CD4 T effector memory HLA-DR+CD38+ expression (Figure 2B and 2C). In general, total and intact HIV DNA were identified to be predictive of similar immune features (Figure S6, Tables S6-S8); how-ever, total HIV DNA generally had greater magnitude ΔR² than intact HIV DNA. The percentage of intact HIV DNA did not appear to contribute greatly to the predictive power of the model for the variance in the immunophenotypes examined (Figure S6C, Tables S10, S11), most likely due to its relatively higher correlation with years of ART (Figures 1C). Nonetheless, the LOCO inference findings suggest a complex biological concert of immunologic changes that occur over years of life and suppressive ART, and which may be influenced by the HIV reservoir.

Receiver operator characteristic curve analysis demonstrates distinct immune parameters associated with intact and total reservoir frequencies

The complex relationship between the HIV reservoir and the immune system led us to hypothesize that machine learning (ML) approaches that combine several parameters into models of the HIV reservoir could be a useful approach to understand this interaction. Before this could be attempted, we decided to first identify a defined set of the most valuable variables from which ML models could be built. To achieve this, we first binarized the 3 reservoir metrics (total, intact, %intact) into high (above median) or low (below median) reservoir groups and generated receiver operator (ROC) curves for the 144 clinical-demographic and immunophenotype parameters (Figure 3). The area under the curve (AUC) of the ROC curve indicates the ability of the feature (immune parameter or clinical/demographic information) to correctly identifying a participant as having qualitatively low or high total HIV DNA (Figure 3A), intact HIV DNA (Figure 3B), or percentage intact HIV DNA (Figure 3C). A model that randomly guesses high or low reservoir has an AUC of 0.5. In this analysis, 57 host variables had an AUC higher than 0.6 for one or more of total HIV DNA, intact HIV DNA, or percent intact DNA (Table S9). The most effective individual immune markers for classification of high versus low total HIV DNA included %NKG2A+ CD4 T (AUC = 0.70), %PD-1+ Tn CD4 T (AUC = 0.68), and %CD38+/HLA-DR-CD8 T (AUC = 0.68). In contrast, the most effective markers for classifying based on the frequency of intact HIV DNA included %CD127+ CD4 T (AUC = 0.71), %CD8 T (AUC = 0.66), and %Tn CD8 T (AUC = 0.65). Finally, for percentage intact HIV DNA, years of ART (AUC = 0.72), %CD107-IFNγ-IL-2+TNFα+ CD4 T (AUC = 0.65) and %KLRG1+/PD1+ CD4 T (AUC = 0.65) were the most effective (Figure 3).

Figure 3:

Similar to the findings in the Spearman correlations (Table 2) and LOCO inference analysis (Figure 2), there were more variables that had AUC values > 0.6 for total reservoir size (Table S9): 44 variables above this threshold for the total reservoir, but only 19 variables for the intact reservoir and 23 variables for percent intact. Overall, this approach allowed us to derive a ranked list of the most predictive immune parameters for each aspect of the HIV reservoir, and these highly ranked features were thus used for subsequent dimension reduction and ML modeling.

Dimension reduction reveals two clusters of PWH with distinct HIV reservoirs

To examine the overall structure of dataset we employed an unsupervised dimension reduction (DR) machine learning approach (PaCMAP (51)) using the 10 immune cell features with the highest AUC values for classifying participants based on total, intact, or percentage intact HIV DNA (Figure 3, Table S9). Interestingly, while no clear clustering was observed when we use intact or percentage intact HIV DNA associated features, we observed two distinct clusters of PWH (cluster 1 and cluster 2, Figure 4A) when using total HIV DNA associated features. These clusters were of roughly equal size with 50 participants in cluster 1 and 65 in cluster 2. Projecting total reservoir size (above or below the median) onto the clustering plot, we observe a strong distinction between the two clusters in terms of total HIV DNA reservoir characteristics, (Figure 4B,C). Analysis of quantitative reservoir frequencies across clusters (Figure 4C-E) demonstrated that cluster 1 is characterized by a smaller absolute intact reservoir frequency, but greater percentage of intact proviruses (Figure 4C-E). In contrast, cluster 2 is defined by a larger absolute intact reservoir frequency, but lower percentage of intact proviruses (Figure 4C-E).

Figure 4:

To gain insight into the immune cell features that distinguish these two clusters, we generated ROC curves for the 133 immune parameters based on their ability to identify membership in cluster 1 versus cluster 2 for each participant (Figure 4F). Many features that distinguished the clusters overlapped with the features selected for the clustering analysis, including expression patterns of CD38, HLA-DR, KLRG-1 and PD-1 on CD4 and CD8 T cells and CD4 T naïve cell frequencies. However, some novel features emerged, including CD4 KLRG1-CD27+ frequency, CD8 KLRG1-PD1-frequency, and CD8 Tn frequency (Figure 4F and Figure S8). Cluster 2, which had a higher total and intact absolute reservoir size but lower percentage of intact proviruses, was associated with a lower fraction of naïve CD4 and CD8 T cells. Additionally, cluster 2 had generally higher expression of immune exhaustion markers, including KLRG1 and PD-1 (Figure S8). These data indicate that cluster membership is highly associated with the overall level of immune activation and exhaustion for the individual. We have previous shown that cannabis use reduces the overall level of immune activation and inflammation in PWH, despite having minimal impact on the size of the HOV reservoir (52). Interestingly, a subset of participants in this study had clinical data regarding cannabis (CB) use. CB users were enriched in cluster 1, which had lower levels of immune exhaustion markers (Figure S8). Notably, with regard to other clinical variables, cluster 2 had an older median age, longer time on ART, lower CD4 nadir, and lower current CD4 T cell count (Figure S8).

Overall, these data demonstrate that, within the cohort of PWH on long-term ART, there were two major clusters of participants with distinct reservoir characteristics and immunophenotypes. These two clusters may reflect the results of a complex concert of age-related immune changes, immune recovery following years of ART, and decay of the intact HIV reservoir, among other possibilities.

Decision tree visualization of PWH with respect to reservoir characteristics

Since the interaction of the immune system and the HIV reservoir is multifactorial, we hypothesized that models that consider multiple parameters simultaneously could more accurately describe the overall dataset, and provide insights regarding the biology of the HIV reservoir and the host immune system. To accomplish this, we employed a decision tree approach to visualize combinations of variables that classify participants as having high (above median) or low (below median) reservoir size. As compared to the dimension reduction technique, decision tree visualization is an interpretable supervised approach and does not require post-hoc analysis.

We first selected 35 variables with the highest ROC AUC values for either total, intact, or percentage intact HIV DNA to be considered for model generation (Figure 4). Using these parameters, we fitted Generalized and Scalable Optimal Sparse Decision Trees (GOSDT) (53) to the data in order to classify high versus low total HIV DNA, intact HIV DNA, or percentage intact HIV DNA. We required the trees to achieve at least 80% accuracy for classifying PWH in the cohort, as well as have at least five PWH in each leaf. Since these trees are based on the entire dataset, these models are thus descriptive rather than predictive.

For the total HIV DNA decision tree, only four immune variables were required to accurately describe high versus low HIV DNA status (Figure 5A, Figure S9): CD8 T cell frequency, CD4 nadir, %CD38+HLA-DR-CD8 T cells, and %NKG2A+ CD4 T cells. The tree divided the cohort into five sub-groups (leaves), among which three have high total reservoir size and two have low total reservoir size. Comparing the labels provided by GOSDT model with the actual data, the tree achieved 83.5% accuracy (i.e. misclassifying 19 PWH among the overall cohort of 115 PWH). Notably, when we combined all samples from “high total reservoir” leaves and all samples from “low total reservoir” leaves, we observed a significant difference in the actual median total reservoir size for these two groups (266/M for low total and 1288.5/M for high total, Mann Whitney U test p-value is 3.56e-13, Figure 5B).

Figure 5:

We then repeated this approach for classifying the cohort participants with respect to intact reservoir size (Figure 5C, Figure S10). Despite the overall lower correlation between the immune cell phenotypes and the intact reservoir size, this tree nevertheless achieved 82.6% accuracy. The tree had six leaves (three with low intact reservoir, three with high intact reservoir) and relied on %Tn CD8 T cells, %CD107-IFNγ+IL-2+TNFα-CD4 T cells, %CD127+ CD4 T cells, %CD38-HLA-DR+ Tcm CD8 T cells, and CD4 Nadir. When we combined the high intact leaves together and the low intact leaves together, we observed a significant difference in median intact reservoir frequency between the groups (29.5/M for low intact and 101/M for high intact, Mann Whitney U test p-value is 9.59e-08, Figure 5D).

Finally, a third iteration of GOSDT generation was performed, classifying the cohort participants with respect to percentage intact of total proviruses (Figure 5E, Figure S11). This tree achieved 82.4% accuracy and relied on years of ART treatment, %KLRG1+CD27-CD4 T cells, %CD127+ CD4 T cells, and recent clinical CD4 count. This tree has two leaves with low percent intact reservoir and three leaves with high percent intact reservoir. We observed a significant difference in the actual percent intact values between the model generated high percent intact and low percent intact leaves (median 5.7% for low percent intact leaves and median 15.3% for high percent intact leaves: Mann Whitney U test, p-value is 5.5e-10, Figure 5F).

Collectively, these decision trees highlight combinations of immune parameters that can accurately describe qualitatively high versus low total, intact, and percentage intact HIV DNA in a cohort of n = 115 PWH. This visualization serves as a basis for mechanistic hypotheses about the interactions of the immune system and HIV reservoir during long-term ART.

Machine learning algorithms identify immune feature combinations that predict high versus low reservoir metrics with ∼70% test accuracy

Although clustering and decision tree analysis permit visualization and understanding of global structures within a dataset, we were curious if combinations of immune and clinical-demographic parameters could actually accurately predict (rather than only visualize) whether a given participant had a high (above median) or low (below median) values with respect to total HIV DNA, intact HIV DNA, or percent intact HIV DNA. We considered five machine learning algorithms including Logistic Regression with L2 regularization (LR), CART, Support Vector Machines with RBF kernel (SVM), Random Forest (RF), and Gradient Boosted Trees (GBT). For reservoir size characteristics, we measured accuracy of the models over 10 random splits of the data into training and test sets (Figure 6A and Table S10). Overall, we found that for total reservoir size Logistic Regression achieved highest mean classification accuracy (69.31%) in test data. The accuracy of these models is likely limited by the sample size (n=115) and noise in the data, as evidenced by lower R² score of machine learning models for non-binarized HIV characteristics (Figure S12). To examine the contribution of individual immune features to model performance, we examined logistic regression coefficients for each immune cell variable in the model. Since Logistic Regression coefficients are associated with expected change in log odds (based on log_e), we can think about the coefficient β for variable X in the following way: increasing variable X by one unit multiplies the odds of high reservoir size (probability that the reservoir size is high divided by the probability that the reservoir size is low) by e^β.

Figure 6:

In Figure 6B we visualize the logistic regression model for one data split among ten we considered for total reservoir size. For this split we observe that higher values of %NKG2A+ CD4 T, %PD-1+ Tn CD4 T, and %Tcm CD8 T are associated with a increased probability of total reservoir size being high. On the other hand, an increase in %Tn CD4 T decreases the odds of high reservoir size. The model visualized in Figure 6B achieved 75.86% training and 75% test accuracy.

Similarly, logistic regression models performed well for predicting intact reservoir size compared to other methods (average 65.17% test accuracy, see Figure 6C). In Figure 6D we visualize the Logistic Regression model for one fixed data split. In this model we observed that higher values of %CD107-IFNγ+IL-2+TNFα-CD4 T and %CD127+ CD4 T, and % CD4 T are associated with lower probability of intact reservoir size being high, while higher values of %KLRG1+CD27+ CD8 T are associated with increased probability of intact reservoir size being high. For analysis of percentage intact HIV DNA, we display training and test accuracy values in Figure 6E. For the visualized model for one data split (Figure 6F), an increase in %CD127+ CD4 T, years of ART, %KLRG1+PD-1+ CD4 T, and %Tem CD4 T leads to a lower probability of percent intact HIV DNA being high.

Overall, these analyses demonstrate that we can use machine learning tools to construct models that can predict with approximately 70 percent accuracy whether a given PWH has qualitatively low or high total, intact or percentage of HIV DNA. The coefficients identified in each of these models, including NKG2A and CD127 expression on CD4 T cells, identify specific immunologic nodes that serve as a basis for the generation of hypotheses about the interplay between total versus intact HIV DNA and the host the immune system.

Discussion

In this cross-sectional study of peripheral blood samples from a well-characterized clinical cohort of 115 PWH on long-term ART, we used IPDA and high parameter flow cytometry to better understand the associations between the HIV reservoir and a broad range of immunophenotypes related to T cell differentiation, activation, exhaustion, senescence, HIV-specific T cell responses, as well as NK cell phenotypes. Using data science approaches, specific HIV reservoir features (total, intact, percentage intact HIV DNA) and clinical-demographic variables (age, sex, CD4 nadir, recent CD4 count, years on ART, years of HIV infection prior to initiating ART) were identified that were particularly important for predicting specific immunophenotypes in PWH on long-term ART. Distinct host immune correlates of intact versus total HIV DNA were also identified. Further, machine learning approaches enabled identification of combinations of immune parameters that accurately classified and predicted whether a given PWH is likely to have qualitatively high or low total, intact, or percentage intact HIV DNA. Collectively, this study provides novel insights into HIV reservoir biology and notably also demonstrates the utility of specific data science and machine learning approaches for HIV reservoir studies. Development and application of these analytic approaches are important to enable biological insights from high-parameter studies of the HIV reservoir and other complex biology, particularly with the advent of high-dimensional techniques such as CYTOF and scRNAseq.

Previous studies have sought to identify immune and other correlates of the total HIV DNA reservoir (26–45). These studies have identified several activation/exhaustion markers that correlate with HIV reservoir frequency, including T cell expression of HLA-DR, CD38, PD-1, CTLA-4, and LAG-3, the CD4/CD8 ratio, and NK cell KLRG1 expression. However, not all studies have found associations of total HIV DNA and immune subsets, particularly for CD4 T cell HLA-DR and CD38 expression(40, 54–56) during long-term ART. The correlative results from this present study corroborate many of these studies, and provide additional insights. A number of correlations between total HIV DNA and immune activation and exhaustion were identified. There were weak to moderate positive correlations observed between total HIV DNA and CD4 T cell expression of NKG2A, PD-1, KLRG1, CD27, HLA-DR as well as CD8 T cell expression of PD-1, CCR7, and CD38. Total HIV DNA was inversely correlated with CD8 T cell expression of CD38, CD4 T cell expression of CD27, the frequency of naïve CD4 and CD8 T cells, and the CD4/CD8 ratio.

Interestingly, when we examined immune correlates of intact HIV DNA measured by the IPDA, many of the same relationships as described above were noted, but with decreased magnitude and lacking statistical significance, with the exception of the inverse correlation observed for CD4/CD8 T cell ratio and intact HIV DNA. This is consistent with a recent study using the IPDA that did not find associations of intact HIV DNA with CD4 T cell HLA-DR and CD38 co-expression (57), as well as a study using nearfull-length sequencing of single cells that found only slight enrichment of HLA-DR or PD-1 in cells with genetically intact, inducible proviruses (37).

Uniquely, however, CD127 expression on CD4 T cells was significantly inversely associated with intact reservoir frequency. CD127 is the alpha-chain IL-7 receptor and is expressed on long-lived memory T cells that are depleted during untreated HIV infection. Interestingly CD127+ cells have recently been found to be highly susceptible to latent HIV infection in tissues (58), and IL-7 is associated with slower natural reservoir decay (27). Furthermore, IL-7 signaling is known to drive expression of VLA-4, which was recently reported to highly enrich for genetically intact, inducible proviruses (37). This observational finding that intact but not total HIV DNA inversely correlates with CD127 expression on CD4 T cells requires further investigation. Finally, this initial immune correlate analysis demonstrated that the percentage of intact HIV DNA was negatively correlated with years of ART as well as and HIV-specific CD4 T cell responses. These findings are consistent with known selective decay of intact proviruses over years of ART therapy (57, 59, 60).

Dimension reduction machine learning approaches identified two robust clusters of PWH when using total HIV DNA reservoir-associated immune cell frequencies (Figure 4A), but not for intact or percentage intact HIV DNA (Figure 4B and 4C). This pattern is likely indicative of a more profound association of total HIV DNA with host immunophenotype relative to intact HIV DNA. This is certainly plausible giving the growing appreciation for the immunogenicity and high frequency of defective proviruses within the HIV reservoir (61–64). These clusters exhibited distinct reservoir characteristics, with one cluster being enriched with participants with a larger total reservoir size, a lower percentage of intact viruses, a lower frequency of naïve T cells and elevated expression of activation and exhaustion markers, while the other contained PWH with smaller total reservoirs, higher percentage of intact viruses, higher frequency of naïve T cells and lower expression of activation and exhaustion markers. The basis for the two distinct clusters is unclear, but the cluster membership was notably associated with age and time on therapy. The existence of two distinct clusters of PWH with different immune features and reservoir characteristics could have important implications for HIV cure strategies – these two groups may respond differently to a given approach, and cluster membership may need to be considered to optimize a given strategy. It is tempting to speculate that these two clusters represent differing severity of the initial immune insult of untreated HIV infection, and possibly dichotomous degrees of immune recovery following initiation of ART. These two clusters are also interesting in the context of a previous study examining reservoir characteristics which proposed two major types of viral reservoir within PWH. One type with smaller reservoirs that were enriched in Tcm cells and associated with lower levels of Ki67+ cells and immune activation, and another type with larger reservoirs that were enriched in Ttm cells and associated with higher levels of Ki67+ cells and immune activation (27). It is possible that these two proposed types of reservoirs represent the two clusters we observe using immune markers and dimension reduction, although we were unable to examine the frequency of HIV proviruses specifically within Tcm and Ttm cells in this study.

In addition to unsupervised dimension reduction approaches, we also found that we were able to construct simple, interpretable decision trees that could describe our cohort of PWH with more than 80% accuracy with respect to the size of the reservoir. Furthermore, when we constructed decision trees to classify PWH separately based on total reservoir, intact reservoir and percentage intact, distinct sets of features were important (though some were the same) for each of these models. The decision tree analyses highlight complex immune parameter combinations that serve as a basis for unique hypotheses about HIV reservoir biology. For instance, lower levels of HIV-specific CD4 T cell responses and lower levels of CD8 Tcm HLA-DR expression were useful for identifying individuals with higher frequencies of intact HIV DNA. This suggests that individuals with higher levels of intact proviruses may have lower levels of HIV-specific immune responses and memory CD8 T cell activation that could facilitate clearance of intact proviruses. Notably, the same relationships were not observed for the total HIV DNA decision tree, highlighting the likely more potent immunogenicity of intact proviruses.Another striking example of machine learning approaches identifying potentially insightful reservoir biology is the percentage intact HIV DNA decision tree. As expected, years of ART was critical to predict the percentage of intact proviruses, but CD4 T cell CD127 expression, CD4 T cell KLRG expression and recent CD4 count also played a role. Notably, following approximately 7 years of ART, lower CD4 T cell CD127 expression or low recent CD4 T cell count were useful to classify individuals with higher percentages of intact proviruses. We hypothesize that this finding suggests incomplete immune recovery (low CD4 T cell CD127 expression and low CD4 T cell count) is associated with a larger percentage of intact proviruses, perhaps because of a lack of dilution with de novo naïve CD4 T cells over many years of ART.

Logistic regression models that we trained to classify high versus low total reservoir, intact reservoir, and percentage intact were able to predict out-of-sample with ∼70% accuracy on average. Each of these models relied on five features only and were simple enough for us to visualize and assess how important each variable is for every characteristic of reservoir size. These models corroborated the potential biological relationships identified in the decision tree analysis, and confirmed that these variables could actually be used to predict out-of-sample whether total, intact or percentage intact HIV DNA was qualitatively high or low.

Our findings should be considered in the light of some inherent limitations and caveats. This study is cross-sectional in nature and is primarily observational, so caution should be used interpreting findings associated with time on therapy. As with any observational study, confounding bias may influence our results, particularly for correlative analyses. Another limitation to this study is that our analyses are limited to peripheral blood, and may not represent HIV reservoirs and immune cells located in tissues known to harbor latent HIV such as gut, spleen and brain.

Nevertheless, our study does have some strengths compared to previous work in this area, specifically a highly detailed flow panel, allowing us to quantify numerous immune subsets with a high degree of resolution. Also, the use of the intact proviral DNA assay allows us to separately examine the association of the immune system with intact and total viral DNA. It is however, important to note that IPDA does not measure all proviral sequences in PWH and is likely failing to identify some sequences due to proviral polymorphisms or dual deletions in the two IPDA amplicons. Also, a minor fraction of IPDA+ proviruses are actually defective due to the limited specificity of using two short amplicons to select for full-length intact proviruses (48, 65).

Overall, these findings suggest a complex concert of immune recovery, HIV reservoir dynamics, and intrinsic host factors (age, biological sex) that shapes host immunophenotype, even after years of ART. Mechanistic work will be needed to fully dissect the dynamic relationship between the immune system and the reservoir. Nevertheless, these findings support a model in which ongoing interaction between the HIV reservoir and the host immune cells continue to drive an association, albeit a relatively minor one, with persistent CD4 T cell activation during long-term ART. Notably, both intact and defective proviruses appear to contribute to this immune signature of HIV latency during long-term ART.

This study identifies several areas for future investigation. The question of whether the HIV reservoir directly drives persistent immune activation could be addressed directly by specifically suppressing viral transcription during therapy and examining the impact on the immune system. Future clinical trials with effective latency reversal and clearance regimens could include measurement of HIV-specific T cell responses, as this analysis suggests that lower levels may be associated with reduced intact HIV DNA clearance.

There is also growing appreciation for additional factors that correlate with the HIV reservoir. Specifically, the size of the reservoir has been linked to the abundance of specific bacteria species in the intestinal microbiome (66). Also, genetic studies have indicated the existence of single nucleotide polymorphisms (SNPs) that are associated with reservoir size, including SNPs in MX1 and IRF7 (67–69). These, along with future studies exploring the interplay of the host immune system and HIV reservoir during long-term might leverage new high-dimensional technologies such as scRNAseq and mass cytometry. The machine learning approaches described and applied in this study may be particularly useful for gaining biological insights into these high parameter datasets, particularly using analysis of global patterns in the data that may be otherwise imperceptible to humans.

Methods

Cohort and sample collection

115 people with HIV (PWH) were recruited from two clinical sites. 66 PWH were recruited at Duke University medical center and 49 PWH at UNC Chapel Hill. Peripheral blood was collected and peripheral blood mononuclear cells (PBMCs) were obtained by Ficoll separation, then frozen in 90% fetal bovine serum 10% dimethylsulfoxide (DMSO) (52, 59, 70, 71).

Flow cytometry

Cell preparation and staining followed previously described methods (72). All antibodies were titrated to optimize signal-to-noise ratio on PBMCs prior to use, assuming a 50μL staining volume. Cell viability was examined using Zombie-NIR Fixable Viability Dye (0.4μL per 50μL staining volume; Biolegend). Antibodies used for surface staining were as follows: KLRG1-BV421 (SA231A2 clone, Biolegend, CD45RA-PacBlue (H100 clone, BL), CD8-BV570 (RPA-T8 clone, Biolegend), CD127-BV605 (A019D5 clone, Biolegend), CD56-BV650 (5.1H11 clone, Biolegend), CCR7-BV711 (G043H7 clone, Biolegend), CD27-BV750 (O323 clone, Biolegend), PD1-VioBright515 (REA165 clone, Miltenyi), NKG2A-PE-Vio615 (REA110 clone, Miltenyi), CD16-PerCP-Cy5.5 (33G8 clone, Biolegend), CD38-PCPeF710HB7 clone, TF), CD14-SparkNIR685 (63D3 clone, Biolegend), CD19-SparkNIR685 (HIB19 clone, Biolegend), HLA-DR-APC-F750 (L243 clone, Biolegend). Antibodies used for intracellular staining were as follows: CD3-BV480 (UCHT1 clone, BD), CD4-PerCP (L200 clone, BD), IFN-g-PE-Cy7 (4S.B3 clone, Biolegend), IL-2-APC (MQ1-17H12 clone, Biolegend), TNFα-AF700 (Mab11 clone, Biolegend). Samples were analyzed using a Cytek Aurora spectral flow cytometer. Flow cytometry analysis was performed in FlowJo v10.8 software.

Intact Proviral DNA Assay (IPDA)

Cryopreserved samples of PBMCs from each study participant were viably thawed. A portion was used for immunophenotyping as described above, and the remainder were subjected to total CD4 T cell negative selection with the StemCell Technologies EasySep^TM Human CD4+ T Cell Enrichment Kit (Cat# 19052). CD4 T cell DNA was extracted using the QIAamp DNA Mini Kit and quantified on a NanoDrop 1000 (Thermo Fisher Scientific). IPDA was performed as originally described (48), with a validated PCR annealing temperature modification to increase signal to noise ratio (59). Gating for positive droplets was set using negative (DNA elution buffer and HIV-seronegative CD4 T-cell DNA), and positive (Integrated DNA Technologies gblock amplicon) positive control wells processed in parallel. DNA shearing index values were similar to those reported previously (median, 0.33; interquartile range [IQR], 0.31–0.34) (59, 70, 71). A median of 1.04 x 106 (Q1 8.71 x 105, Q3 1.16 x x 106) cell equivalents were assessed for each donor.

Statistics

Variable data analysis

To create binarized labels that represent reservoir characteristics, we split the data at the median, which is 553/M for total reservoir size, 53/M for intact reservoir size, and 8.64 percent for percent intact. For every variable (immune cell frequency, demographics, and clinical information), we computed Spearman correlation and an Area Under the Curve (AUC) value. To compute the AUC value, we first created a Receiver Operating Characteristic (ROC) curve by plotting the true positive rate versus the true negative rate for every cell subset frequency. AUC is then computed as the area under the ROC curve using the trapezoidal rule.

To assess the importance of clinical, demographic information and HIV reservoir characteristics in describing the immune markers, we performed a leave-one-covariate-out (LOCO) analysis. We trained a linear regression model (M₁) with intercept based on all variables. We removed each variable at a time and trained a new model without this variable. Then we computed the difference in R² score before and after dropping the variable. Variance Inflation Factors (VIFs) for age, sex, years of ART, CD4 nadir, recent CD4 account, and years of HIV infection prior to HIV infection were less than 2, indicating an acceptable level of correlation between these independent variables (73). However, higher VIF values were observed for intact and total HIV DNA, therefore separate LOCO inference analyses were conducted in order to avoid artificial fluctuations in model fit due to multicollinearity.

Data visualization

We used PaCMAP (51) to reduce the dimensionality of the dataset to a 2-dimensional space. Due to the relatively small sample size of our dataset (n=115), we set the number of neighbors to 5. Additionally, we used GOSDT (53) which computes sparse optimal trees to identify patterns in the dataset. For all the visualization trees we set the depth budget parameter to 4 and the regularization parameter to 0.02.

Training and generalization

Our training procedure for classification is described in Algorithm 1 in Supplemental Materials, which is a standard machine learning pipeline with an added search for the number of variables. For the full hyperparameters list see Table S8. We normalized continuous variables to have values between 0 and 1. For logistics regression models, we set the regularization parameter of .e − 2 regularization based on cross-validation. The values are 10 for total reservoir size and percent intact and 100 for intact reservoir size.

Study approval

Written informed consent was obtained for all study participants. The study design was reviewed and approved by IRB for both Duke University and UNC Chapel Hill.

Data availability

Deidentified datasets will be made publically available upon manuscript publication.

Author contributions

Designed the study – EPB, DMM, CDR. Conducted experiments – SDF, ADVC. Analyzed the data – LS, YW, SDF. Wrote the manuscript: All authors.

Conflict of interest statement

The authors have declared that no conflict of interest exists.

Acknowledgements

This work was supported by the following grants from the National Institutes of Health: NIAID R01 AI143381 (EPB), NIAID UM1 AI164567 (DM Margolis), NIDA R61 DA047023 (EPB), NIDA R01 DA054994 (CDR), NIAID F30 AI145588 (SDF).

Supplementary Methods

CD4/CD8 ratio regression analysis

To examine the connection between CD4/CD8 and (CD127+ CD4)/CD8 ratio and HIV characteristics (Figure S5), we first computed the natural logarithm (log_e) of the intact, total and percent intact as the outcome and normalized outcome of each immune feature (meaning that our values of ratio and outcome are between 0 and 1). We removed outliers using the DBSCAN algorithm (1), which is an unsupervised clustering algorithm that groups data points based on density into a single cluster. For DBSCAN, we set the parameter that determines radius of a circle around each point (that is used to compute density) as 0.15 and minimum number of samples in each cluster as 10. Then we fitted the remaining datapoints using Linear Regression.

Regression machine learning models to predict reservoir size from immune cell frequencies

To predict reservoir size based on immune cell frequencies, clinical and demographics information (Figure S12), similar to classification, we sought to find the set of model parameters and types of models that would be able to fit the dataset without overfitting. We used several approaches: Linear Regression (LR), Ridge Regression (RR), Kernel Regression with RBF kernel (KR), Decision Tree Regressor (DT), Random Forest (RF), and Gradient Boosted Trees (GBT). We followed the procedure described in Algorithm 1, optimized the R² score, returned the mean and standard deviation of the R² score of training and test sets, and chose features based on the absolute value of the correlation coefficient. We normalized variables to have values between 0 and 1. We then fitted the models to the normalized natural logarithm (log_e) of total, intact reservoir size and percent intact and then further normalization makes outcomes to be between 0 and 1. For total and intact reservoir size we visualized the Ridge Regression model, which performed the best based on mean test R² score (Figure S12 D, E, G, H). We set the regularization parameter to 1, as this value was best on average according to cross-validation for total and intact reservoir size. For percent intact, we visualized Linear Regression based on 4 variables (Figure S12 F, I) as it achieved the second highest R² score after Kernel Ridge based on 25 variables.