Identifying Longitudinal CD4:CD8 Ratio Trajectories Indicative of Chronic Renal Disease Risk among People Living with HIV: An Application of Growth Mixture Models

The incidence of chronic kidney disease (CKD) is increasing among people living with HIV (PLWH). Routine monitoring of indicators such as CD4:CD8 ratio might improve the early detection of CKD. Our objective was to identify clinically relevant CD4:CD8 ratio trajectories indicative of CKD risk. Participants were ≥ 18 years old, initiated antiretroviral therapy between 2000 and 2016, and were followed for ≥6 months until 31 March 2017 or last contact date. Outcome was incidence of CKD. Growth mixture models (GMMs) and decay models were used to compare CD4:CD8 ratio trajectories. Following GMM, 4547 (93.5%) participants were classified in Class 1 with 5.4% developing CKD, and 316 (6.5%) participants were classified in Class 2 with 20.9% developing CKD. The final model suggested that participants in Class 2 had 8.72 times the incidence rate of developing CKD than those in Class 1. Exponential decay models indicated a significant CD4:CD8 ratio decline among Class 2 participants who developed CKD. Among those who developed CKD in Class 2, starting at 5.5 years of follow-up, the slope of their ratio trajectory curve changed significantly, and the rate of decline increased dramatically. Routine monitored CD4:CD8 ratios can be an effective strategy to identify early CKD risk among PLWH.


Introduction
Since the introduction of antiretroviral therapy (ART), the life expectancy of people living with HIV (PLWH) has increased, which has led to higher incidence of age-associated conditions such as non-AIDS-related malignancies, liver, cardiovascular, cerebrovascular, and renal diseases in this population [1][2][3]. The complex process of aging, while living with HIV, presents many challenges for patients and healthcare providers, due to different psychosocial (e.g., social isolation), structural (e.g., access to healthcare resources, housing, food, and medications), and behavioral factors (e.g., substance use), and the cumulative toxicity of HIV/non-HIV medications increasing the complexity in the clinical management of individuals with multimorbidity [4]. Studies have shown that the management of individuals who present with a large number of medical issues can be quite complex, increasing the risk of polypharmacy and drug-drug interactions [5][6][7][8]. Since services and programs for PLWH with multimorbidity in different settings can be highly disconnected,

Outcome and Study Covariates
The outcome was incidence of CKD identified using a case finding algorithm published by the BC Ministry of Health based on physician claims and hospitalizations with a diagnostic code for CKD Supplementary Materials, Table S1 [60]. To estimate incidence, we applied a washout period of five years before first ART date. For participants lacking five years of information prior to their first ART date, a washout period of five years after their earliest encounter with the healthcare system was used. CD4:CD8 ratios were measured every six months.
Time-invariant covariates measured at baseline included: sex (female, male), age (continuous), ART naïve (yes, no), and a categorical indicator reflecting a cohort effect (years 2000-2004, 2005-2010, 2011-2016). Time-invariant covariates measured at the end of follow-up included: history of substance use disorder (SUD) (yes, no), duration of hypertension, diabetes mellitus, and cardiovascular disease (CVD) (no disease, < median duration, ≥ median duration), area under the HIV viral load curve [63], proportion of time on protease inhibitors (PI), proportion of time on non-nucleoside reverse transcriptase inhibitors (NNRTI), proportion of time on tenofovir disoproxil fumarate (TDF), and followup time (continuous). CD4 nadir (continuous) was measured one year before the end of follow-up. Case-finding algorithms used to identify participants with SUD, hypertension, diabetes mellitus and CVD are described in the Supplementary Materials Table S1 Please note that, beside TDF, we did not investigate other specific antiretrovirals known to increase the risk of CKD [20][21][22]. The variables included in this study were chosen a priori based on the clinical literature and treatment guidelines describing factors known to influence CKD risk.

Statistical Analyses
We conducted bivariable analyses between risk factors and CKD using Pearson's chi-square or Fisher's exact test for categorical variables (depending on cell count), and Wilcoxon's Rank Sum test for continuous variables. This study consisted of a two-pronged analysis:

Growth Mixture Modeling
To investigate the patterns of CD4:CD8 ratio trajectories among participants with and without CKD, we applied growth mixture modeling (GMM) [52,64,65]. This type of modeling allows the identification of unobserved subpopulations or classes based on distinctive trajectories of CD4:CD8 ratios. First, we fit unconditional models with a different number of latent classes, and with linear (in case the ratio changes linearly over time) and Viruses 2023, 15, 385 4 of 14 quadratic (in case there is a non-linear pattern over time) growth factors; these terms were introduced one at a time. Second, we focused on examining the relationships between covariates and latent classes by fitting the GMM including covariates. Adding covariates to the model helps us understand the key characteristics of each group and how different covariates may affect the longitudinal trajectory within each class. To examine the fitness of the models, we used the Akaike information criteria (AIC), Bayesian information criteria (BIC), sample size-adjusted BIC (SBIC), entropy (range 0 to 1), adjusted Lo-Mendell-Rubin likelihood ratio test (LMR-LRT), and class size. A LMR-LRT with p-value of <0.05 was deemed significant while lower values of AIC, BIC and SBIC indicated better fit [66,67]. Third, we compared the incidence of CKD using two versions of the selected model, one following an unconditional approach (without covariates) and the other including covariates. To analyze the association between the outcome (CKD incidence) and the latent classes, we assigned participants to their most likely classes and then included this variable as an explanatory variable in the Poisson regression to model the incidence CKD rate per 1000 person-years between the distinct trajectories. Last, to understand the key characteristics of participants in each class and how different covariates may affect the CD4:CD8 ratio trajectories, we conducted bivariable analyses between covariates and the latent classes stratified by CKD status. Analyses were conducted with SAS version 9.4 (SAS Institute Inc., Cary, NC, USA) and Mplus v.8.1 (Los Angeles, CA, USA).

Decay Model
Based on the results from the GMM, we estimated CD4:CD8 ratio trajectories to assess whether there are early indications of participants being at a high risk of CKD and potential cut-off time points of intervention. CD4:CD8 ratio measurements were obtained longitudinally from study baseline until the end of follow-up. We modeled and compared the CD4:CD8 ratio depletion trajectories among the identified latent classes and CKD status. The natural logarithm of CD4:CD8 ratios were modeled using a non-linear mixed effects model [68,69]. The follow-up period was modeled using exponential decay function, assuming a Gaussian distribution (identity link), an AR(1) correlation structure, and a random intercept. The exponential decay model is described as follows: where t represents each of the six-month intervals, i represents each participant at the study, ε i is the random error distributed as Normal (0, D i ), where D is the covariance matrix and b 0i as N (0,

Growth Mixture Modeling
To investigate the patterns of CD4:CD8 ratio trajectories among participants with and without CKD, we applied growth mixture modeling (GMM) [52,64,65]. This type of modeling allows the identification of unobserved subpopulations or classes based on distinctive trajectories of CD4:CD8 ratios. First, we fit unconditional models with a different number of latent classes, and with linear (in case the ratio changes linearly over time) and quadratic (in case there is a non-linear pattern over time) growth factors; these terms were introduced one at a time. Second, we focused on examining the relationships between covariates and latent classes by fitting the GMM including covariates. Adding covariates to the model helps us understand the key characteristics of each group and how different covariates may affect the longitudinal trajectory within each class. To examine the fitness of the models, we used the Akaike information criteria (AIC), Bayesian information criteria (BIC), sample size-adjusted BIC (SBIC), entropy (range 0 to 1), adjusted Lo-Mendell-Rubin likelihood ratio test (LMR-LRT), and class size. A LMR-LRT with p-value of <0.05 was deemed significant while lower values of AIC, BIC and SBIC indicated better fit [66,67]. Third, we compared the incidence of CKD using two versions of the selected model, one following an unconditional approach (without covariates) and the other including covariates. To analyze the association between the outcome (CKD incidence) and the latent classes, we assigned participants to their most likely classes and then included this variable as an explanatory variable in the Poisson regression to model the incidence CKD rate per 1000 person-years between the distinct trajectories. Last, to understand the key characteristics of participants in each class and how different covariates may affect the CD4:CD8 ratio trajectories, we conducted bivariable analyses between covariates and the latent classes stratified by CKD status. Analyses were conducted with SAS version 9.4 (SAS Institute Inc., Cary, NC, USA) and Mplus v.8.1 (Los Angeles, CA, USA).

Decay Model
Based on the results from the GMM, we estimated CD4:CD8 ratio trajectories to assess whether there are early indications of participants being at a high risk of CKD and potential cut-off time points of intervention. CD4:CD8 ratio measurements were obtained longitudinally from study baseline until the end of follow-up. We modeled and compared the CD4:CD8 ratio depletion trajectories among the identified latent classes and CKD status. The natural logarithm of CD4:CD8 ratios were modeled using a non-linear mixed effects model [68,69]. The follow-up period was modeled using exponential decay function, assuming a Gaussian distribution (identity link), an AR(1) correlation structure, and a random intercept. The exponential decay model is described as follows: where t represents each of the six-month intervals, i represents each participant at the study, ε is the random error distributed as Normal (0, Di), where D is the covariance matrix and as N (0, 2 ). The exponential decay function was modeled through * , where and R are coefficients of this function. Loss rate was modeled using the following equations.
). The exponential decay function was modeled through γe −R * t , where γ and R are coefficients of this function. Loss rate was modeled using the following equations.
where r is the CD4 : CD8 ratio, dr dt is the CD4 : CD8 rate and dy dt is the ln(CD4 : CD8) rate. These analyses were conducted in R© version 4.0.5.

Results
The STOP HIV/AIDS cohort includes 15,599 participants, of which 8286 were excluded since they were not eligible for the study. Thus, out of the remaining 7313 patients, 2450 were eliminated since they did not fulfill the data requirements for the modeling. Therefore, 4863 participants were included in our study Supplementary Materials, Figure S1. Those excluded (i.e., 2450) were more likely to be female, with a history of SUD, not naïve to ART, enrolled in the cohort between 2011 and 2016 and of an older age Supplementary Materials,  Table S2. Most of the included participants were male (82.2%), ART naïve (89.6%), without a history of SUD (57.3%), hypertension (87.2%), diabetes mellitus (92.0%), and CVD (93.1%). The median follow-up time was 6.45 (25th, 75th percentiles: 3.57, 9.89) years, and the median number of CD4:CD8 measurements was 8 (25th, 75th percentiles: 4, 15). A total of 310 (6.4%) participants developed CKD throughout the study. Bivariable analyses revealed that participants who developed CKD and those who did not differ across all measured characteristics, except age at the end of follow-up, and the proportion of time on NNRTI Table 1. Participants who developed CKD were more likely to be female, naïve to ART, enrolled in the cohort between 2000 and 2004, with a history of SUD, hypertension, diabetes and CVD, older at study entry, higher baseline HIV viral load, shorter exposure to TDF, and longer exposure to PI.

CD4:CD8 Ratio Trajectories and CKD Status
We built different GMMs to model the CD4:CD8 ratio measurements over time. Based on the goodness-of-fit statistics found in the Supplementary Materials Table S3, the final model included two classes (linear) with time-invariant covariates Supplementary Materials Figure S2. Having two classes means that there are two groups in our population, and they differ in the CD4:CD8 ratio trajectory in a "particular way", which we try to understand based on the following analyses. We also compared the rate of renal disease of the two-class (linear) model with and without time-invariant covariates Table 2. Our results showed that in both models, participants in Class 2 had a higher rate of CKD than those in Class 1. For the final model, i.e., two-class (linear) model with time-invariant, those assigned Class 2 had 8.72 times the rate (95% CI 6.64-11.44) of developing incident CKD than those in Class 1. Based on the final model, we plotted the CD4:CD8 ratio trajectories for each class Supplementary Materials Figure S3 and performed bivariable analyses Table 3. A total of 4547 (93.5%) participants were classified in Class 1 with 244 (5.4%) of them developing incident CKD, and 316 (6.5%) were classified in Class 2 with 66 (20.9%) developing incident CKD. Our bivariable analyses suggested that among Class 2 participants, those who develop CKD were more likely to have history of SUD, be ART naïve, enrolled in the cohort between 2000 and 2004, younger at baseline, with a lower CD4 nadir closer to the end of follow-up, and a shorter time on TDF (p-values <0.05).

CD4:CD8 Ratio Decay among Classes and CKD Status
We fit an exponential decay function model to estimate the CD4:CD8 ratio trajectories by GMM class and CKD status. A total of 60,869 CD4:CD8 measurements were included with 59,191 belonging to Class 1 and 1678 to Class 2. The models and goodness-of-fit assessments can be found in the Supplementary Materials Table S3, Figure S4. Figure 1a,b present the ln(CD4 : CD8) trajectory, and Figure 1c,d present the ln(CD4 : CD8) loss rate over time, for every 6-month period. We observed that the estimated trajectory of ln(CD4 : CD8) ratio of participants who developed CKD differed significantly between Classes 1 and 2. Those in Class 1 experienced a rapid increase and after 7.5 years, the ratio reached a plateau Figure 1a. Please note that although the rate of increase between those with and without CKD was very similar, those with CKD in Class 1 had significantly lower ratios. For those in Class 2, we observed a decline in ln(CD4 : CD8) among those with and without CKD. However, among those who developed CKD, starting around 5.5-6.5 years of follow-up Figure 1b,d, the slope of their ratio trajectory curve changed significantly, and the rate of decline increased dramatically. Data for Figure 1 can be found in Tables S4 and S5, as well as the original values for the ratio in this analysis.  Figure S4. Figure 1a,b present the ( 4: 8) trajectory, and Figure 1c,d present the ( 4: 8) loss rate over time, for every 6-month period. We observed that the estimated trajectory of ( 4: 8) ratio of participants who developed CKD differed significantly between Classes 1 and 2. Those in Class 1 experienced a rapid increase and after 7.5 years, the ratio reached a plateau Figure 1a. Please note that although the rate of increase between those with and without CKD was very similar, those with CKD in Class 1 had significantly lower ratios. For those in Class 2, we observed a decline in ( 4: 8) among those with and without CKD. However, among those who developed CKD, starting around 5.5-6.5 years of follow-up Figures 1b,d, the slope of their ratio trajectory curve changed significantly, and the rate of decline increased dramatically. Data for Figure 1 can be found in Tables S4 and S5, as well as the original values for the ratio in this analysis.

Sex-Based Analysis
We conducted a sensitivity analysis by sex at birth. Please note that to run GMM or any other statistical model, we need power and since we have very few females in our study, we could only run descriptive analyses and fit the decay model for Class 2 only including sex as a covariate. We noticed that the rate of CKD is 53% higher in females (Table 4) and the rate of decline in ( 4: 8) among Class 2 with CKD was also more accelerated among females (Figure 2).

Sex-Based Analysis
We conducted a sensitivity analysis by sex at birth. Please note that to run GMM or any other statistical model, we need power and since we have very few females in our study, we could only run descriptive analyses and fit the decay model for Class 2 only including sex as a covariate. We noticed that the rate of CKD is 53% higher in females (Table 4) and the rate of decline in ln(CD4 : CD8) among Class 2 with CKD was also more accelerated among females (Figure 2).

Discussion
Using routinely monitored tests such as the CD4:CD8 ratio, this study identified clinically relevant subpopulations of PLWH who are at increased risk of incident CKD. Our novel use of GMM allowed us identify higher CKD rates among participants classified into the Class 2 trajectories. Furthermore, our subsequent decay functions and CD4:CD8 ratio loss rates also suggested that, among Class 2 participants, there was a significant CD4:CD8 ratio decline among those who developed CKD starting at 5.5-6.5 years of follow-up (a CD4:CD8 ratio <0.10); or a decline of 45% from baseline values. This finding is interesting and puzzling. The usual pattern for patients is to improve their CD4:CD8 ratio once they start ART [70]. For those with advanced HIV prior to the start of ART (e.g., low CD4 nadir, prior opportunistic infection), the improvement usually seen in the CD4:CD8 ratio is minor [2,70]. In this study, CKD incidence was higher among participants with the poorest ratio response in both classes, especially those in Class 2. We also observed that in our study, females had a higher CKD incidence, and the rate of decline in CD4:CD8 was more accelerated than in males. Some studies have shown that females, especially done of child-bearing age may experience higher rates of CKD due to biological factors as well as inequities in accessing care for early diagnosis of kidney disease. [71,72] One question that comes to mind is whether HIV-related factors drove the poor response seen in the CD4:CD8 ratio. In our study, the AUCVL in all groups did not indicate that poor virologic control influenced our results. Additionally, the effect may have been driven by a loss of memory CD4 from living with HIV for a longer time, as indicated in our results. Third, we controlled for the effect of CD4 nadir and antiretrovirals known to be risk factors for CKD, and the our CD4:CD8 GMM model was still able to distinguish those more likely to develop CKD. Last, suppose that non-HIV-related factors are involved in the observed trends in the CD4:CD8 ratio, especially in Class 2. In this case, the question that remains is, as CKD progresses, is there a biological mechanism that predicts immunologic non-response among PLWH?
The use of eGFR and ACR has been identified as gold standard to identify and monitor CKD [73]. However, these measurements can yield low sensitivity in detecting early CKD, especially among PLWH, participants with certain comorbidities (e.g., chronic liver disease, obesity) [30,74,75], and external factors such as muscle mass, nutritional status, age, and race. Among PLWH, routinely monitored CD4:CD8 ratio measurements have been identified as emerging biomarkers due to its low values being associated with increased risk of morbidity and mortality [39,50,[76][77][78]. Previous studies have focused on using the CD4:CD8 ratio as a marker of greater risk of cardiovascular events, non-AIDS defining cancers, diabetes, and hypertriglyceridemia [37,40,42,46]. These studies have also

Discussion
Using routinely monitored tests such as the CD4:CD8 ratio, this study identified clinically relevant subpopulations of PLWH who are at increased risk of incident CKD. Our novel use of GMM allowed us identify higher CKD rates among participants classified into the Class 2 trajectories. Furthermore, our subsequent decay functions and CD4:CD8 ratio loss rates also suggested that, among Class 2 participants, there was a significant CD4:CD8 ratio decline among those who developed CKD starting at 5.5-6.5 years of follow-up (a CD4:CD8 ratio < 0.10); or a decline of 45% from baseline values. This finding is interesting and puzzling. The usual pattern for patients is to improve their CD4:CD8 ratio once they start ART [70]. For those with advanced HIV prior to the start of ART (e.g., low CD4 nadir, prior opportunistic infection), the improvement usually seen in the CD4:CD8 ratio is minor [2,70]. In this study, CKD incidence was higher among participants with the poorest ratio response in both classes, especially those in Class 2. We also observed that in our study, females had a higher CKD incidence, and the rate of decline in CD4:CD8 was more accelerated than in males. Some studies have shown that females, especially done of child-bearing age may experience higher rates of CKD due to biological factors as well as inequities in accessing care for early diagnosis of kidney disease [71,72]. One question that comes to mind is whether HIV-related factors drove the poor response seen in the CD4:CD8 ratio. In our study, the AUCVL in all groups did not indicate that poor virologic control influenced our results. Additionally, the effect may have been driven by a loss of memory CD4 from living with HIV for a longer time, as indicated in our results. Third, we controlled for the effect of CD4 nadir and antiretrovirals known to be risk factors for CKD, and the our CD4:CD8 GMM model was still able to distinguish those more likely to develop CKD. Last, suppose that non-HIV-related factors are involved in the observed trends in the CD4:CD8 ratio, especially in Class 2. In this case, the question that remains is, as CKD progresses, is there a biological mechanism that predicts immunologic non-response among PLWH?
The use of eGFR and ACR has been identified as gold standard to identify and monitor CKD [73]. However, these measurements can yield low sensitivity in detecting early CKD, especially among PLWH, participants with certain comorbidities (e.g., chronic liver disease, obesity) [30,74,75], and external factors such as muscle mass, nutritional status, age, and race. Among PLWH, routinely monitored CD4:CD8 ratio measurements have been identified as emerging biomarkers due to its low values being associated with increased risk of morbidity and mortality [39,50,[76][77][78]. Previous studies have focused on using the CD4:CD8 ratio as a marker of greater risk of cardiovascular events, non-AIDS defining cancers, diabetes, and hypertriglyceridemia [37,40,42,46]. These studies have also used predetermined CD4:CD8 ratio cut-offs to assess the risk of morbidity and mortality [79]. However, given our PLWH population and their variability in demographic, clinical, and behavioral characteristics, it is key to use an approach that accounts for the heterogeneity of participants, such as GMM [52]. Our GMM allowed us to examine the dynamic changes in CD4:CD8 ratio during treatment and identify the diversity in trajectories while accounting for traditional risk factors such as sex, history of SUD, ART regimen and presence of comorbidities. This model aided in the identification of clinically relevant subpopulations of PLWH with higher CKD risk [46], through the different trajectories in CD4:CD8 ratio. Please note that one of the factors that we explored in this study was time exposed to different classes of ART and TDF. As expected, longer time exposed to PIs was associated with a higher risk of CKD. Please note that, during the study, a large proportion of PLWH received ritonavir-boosted atazanavir followed by ritonavir-boosted lopinavir. We also found that those who developed CKD received TDF for a shorter time. This finding is related to confounding by indication, i.e., for some participants, TDF is rapidly discontinued likely due to abnormal renal test results.
Our study has some important strengths. First, the use of linked datasets comprised of different administrative health databases allowed us to obtain clinical information necessary to account for population characteristics, using GMM, and estimate CD4:CD8 ratio trajectories. Second, we used a relatively long washout period (i.e., five years) to estimate CKD incidence and minimize the risk of capturing prevalent cases. Third, we were able to evaluate cumulative exposure to ART known to increase the risk of CKD. However, there are also some limitations to consider. First, our estimates relied on at least two measurements of CD4:CD8 ratio per participant. Thus, participants with only one measurement were omitted from the study. Second, we were unable to adjust for covariates in our CD4:CD8 ratio decay model due to sample size limitations. However, we conducted goodness-of-fit assessments to choose our models, and the model showed good performance with the standardized residuals and the fitted values being correlated with the observed values. Third, we used AUCVL to summarize the overall viral load trajectories of each participant in our study. The measure as calculated may be subjected to selection bias since participants have different measures of viral load over time [80]. However, we showed in the table comparing the characteristics of those belonging to each of the classes, that the AUCVL was not statistically significantly different between the comparison groups, and therefore, we believe that this bias did not play a significant role in our results. Last, although our models yield good estimates, we did not have access to other factors that may explain the risk of CKD (e.g., race/ethnicity, hepatitis C co-infection, smoking) in our population [30].

Conclusions
This study demonstrated the potential use of the CD4:CD8 ratio as a complementary measure to eGFR and ACR to aid the early identification of PLWH who may be at elevated risk of developing CKD, and could therefore benefit from closer longitudinal monitoring of their renal function.