Nonlinear Trajectory of Estimated GFR Progression in HIV-Positive Patients with Normal Renal Function on Tenofovir - based Therapy in Asia: a retrospective observational cohort study

Background fumarate (TDF) therapy has been widely estimated using linear models, but this linearity assumption is not well validated and the nonlinearity of eGFR change remains unknown in Asians with initially normal renal function. Estimated GFR trajectories were compared by one-linear and piecewise-linear mixed effects models, before and after propensity matching, respectively. Whether the incidence of renal dysfunction (reduced renal function RRF, eGFR < 90 mL/min/1.73 m² and rapid kidney function decline RKFD, eGFR > -3 mL/min/1.73 m²/year) follows nonlinearity was assessed by logistic regression. this Chinese observational cohort, examined 823 (299 of TDF users 524 of non-TDF users) treatment-naïve HIV-1-infected participants years) with first eGFR greater than 90 mL/min/1.73m². The


Abstract
Background Estimated glomerular filtration rate (eGFR) trajectory in HIV-1-infected patients on tenofovir disoproxil fumarate (TDF) therapy has been widely estimated using linear models, but this linearity assumption is not well validated and the nonlinearity of eGFR change remains unknown in Asians with initially normal renal function.

Methods
Estimated GFR trajectories were compared by one-linear and piecewise-linear mixed effects models, before and after propensity matching, respectively. Whether the incidence of renal dysfunction (reduced renal function RRF, eGFR < 90 mL/min/1.73 m² and rapid kidney function decline RKFD, eGFR > -3 mL/min/1.73 m²/year) follows nonlinearity was assessed by logistic regression.

Conclusion
For HIV-1-infected Asians initiated antiviral therapy with normal renal function, the nonlinearity of renal function do exist. Routine screen based on nonlinear progression of eGFR could be helpful for clinical management and patient care.

Background
The widespread use of combination antiretroviral therapy (cART) has essentially improved the life expectancy of human immunodeficiency virus (HIV)-positive individuals [1]. However, lifelong use of potentially nephrotoxic antiretrovirals (ARVs), especially tenofovir disoproxil fumarate (TDF), can cause or exacerbate renal impairment [2,3]. Recent publications pay more attention to improve the quality of life in HIV-1 infected patients. Accomplishing this may necessitate alliances with anticipation, early detection of risk factors through routine screening to maintain renal toxicity vigilance.
Estimated glomerular filtration rate (eGFR) is a common indicator of renal function [4,5]. Studies have consistently demonstrated that TDF, as a conventional component of cART [6,7], is associated with a decline of eGFR and renal dysfunction in a subpopulation [8][9][10][11][12][13]. Delineating exactly the eGFR progression trajectories on TDF therapy through routine screening is undoubtedly helpful in this scenario. Since a linear figure seems convenient to interpret, most of the relevant studies so far considered the decline of eGFR to be approximated linear. The real trajectory of eGFR over time is however missed in these simplified model, thus hinders the optimization of TDF therapy in terms of renal function progression. In chronic kidney disease (CKD) population, several groups have reported nonlinear trajectories of eGFR in the past few years, its implications on risk estimation have gained interest and encouraged researchers to identify time dependent factors associated with this phenomenon in CKD with different origins [14][15][16]. However, no studies from HIV-1-positive patients have yet rigorously assessed the nonlinear changes of eGFR over time, especially in patients with normal eGFR on initiation of TDF-based antiviral therapy.
The objective of this study was to comprehensively analyze the trajectory of eGFR over time, and to compare the impact of regimens with or without TDF on this trajectory, in a Chinese cohort of treatment naïve HIV-1-positive individuals who initiated cART between Jan., 2010 and Dec., 2015. We also assessed the incidence of renal dysfunction based on nonlinear changes in eGFR, by using a twopiecewise logistic regression model. The two outcome definitions of this study were reduced renal function (RRF: eGFR ≥ 90 mL/min/1.73 m 2 at baseline and eGFR < 90 mL/min/1.73 m 2 during follow-up) [20], and rapid kidney function decline (RKFD: with progression to CKD; eGFR decline > 3 mL/min/1.73 m 2 /year, estimated by least squares regression) [21].

Statistical Analyses
Baseline characteristics were recorded and compared between TDF users and non-TDF users, and electronically matched as necessary. Three models were used to analyze eGFR progression over time since ART initiation in each group (see Table 3). Model 1, the crude one, was not adjusted for any covariates. Model 2 was adjusted for age, sex, weight, height, body mass index (BMI), CD4 count, eGFR, dyslipidemia, HIV/AIDS risk factors (sexual orientation and intravenous drug use), WHO stage III/IV HIV/AIDS, hepatitis B positivity, hepatitis C positivity, anemia, diabetes, and HIV-1 RNA viral load.
Model 3 used propensity score matching (PSM) to reduce preexisting imbalances in the covariates and potential confounding [22,23], and a covariate was considered well balanced when theP value was more than 0.05 (see Table 1), more technical details were as in additional Table S1.
The nonlinear trajectories of eGFR were determined by smooth curve fitting using a generalized additive model (GAM). Two methods were used to identify significant time points (inflection points on the smooth curves): one determined whether the difference of segmented slopes was equal to zero by the Wald test; the other applied a log likelihood ratio test to compare a nonlinear regression model with a one-linear regression model (see Table 2). Eventually, the time points were determined by constructing a maximum likelihood model using a recursion method. A two-piecewise linear mixed effects model, with random intercepts, was applied to quantify the average change per year of eGFR during different periods on cART (see Table 3). In addition, a two-piecewise logistic regression model based on Generalized Estimating Equation (GEE) was used to estimate the relationship of cART duration with RRF and with RKFD (see Table 4). All multivariate regression models were adjusted for the covariates used in Model 2.
Data on HIV-1 RNA viral load were not available in up to 50% of patients, so a missing value category was used in the main analyses [24,25]. In addition, to reduce bias caused by exclusion of individuals with any missing data at baseline, 5 imputed datasets (established by multiple imputation with chained equations) were developed and run separately, and the results were combined using Rubin's method [26,27] (see Additional file 1: Table S2, S3). Another sensitivity analysis was conducted to exclude patients receiving protease inhibitors (PIs), because of the possible association of these drugs with nephrotoxicity and impaired renal function [28][29][30] (see Additional file 1: Table S4, S5).
All analyses were performed using the R software, version 3.3.1 (http://www.R-project.org). A result was considered statistically significant when the two-tailed P value was below 0.05.

Patient Selection and Propensity Score Matching
As shown in the flowchart (Fig.1), a total of 1065 patients were enrolled and 823 patients were eligible for participation, 299 of whom (36.3%) started a TDF-containing cART. Table 1

Comparison of one linear and piecewise linear mixed effects models.
We compared eGFR trajectories by one linear and piecewise linear models (see Table 2), with the piecewise model allowing a change of the eGFR slope at a given time point. Log likelihood ratio test between the two models indicated that nonlinear trajectory of eGFR was a better fit than the traditional one assuming a single linear process across the entire period of observation (P < 0.001 for all).

Time points on nonlinear trajectories of eGFR.
For non-TDF users, the time points were 2.55 years (see Table 2, model 1), 2.15 years (see Table 2, model 2), and 2.15 years (see Table 2 years. Similar results were obtained in other models (see Table 2).

The relationship between eGFR and duration of cART.
The eGFR changed over time in both groups (Fig.2, Additional file 1: Fig.S1, S2). There was a reverse S-shaped relationship between eGFR and duration of cART for TDF users, but a different temporal trajectory for non-TDF users, in all three models. The S-shaped trajectory was observed markedly in model 1 (Additional file 1: Fig.S1B) and model 2 (Fig.2B). Table 3 showed average eGFR changes per year for the two groups according to cART duration. For TDF users, we obtained different results when the duration of cART was categorized using different time points in all three models. For cART less than 1.40 years, the exp (β) was -5.31 (95% CI: -6.57, -4.06) and for 2.30 years or more was -3.71 (95% CI: -5.97, -1.45). However, the exp (β) was 4.83 (95% CI: 1.38, 8.28) for 1.40 to 2.30 years. For models 1 and 3, these time points were nearly the same, and similar trends were indicated in eGFR with increasing duration of cART (see Table 3).

Average changes in eGFR over time on different cART duration among TDF or non-TDF users.
For non-TDF users, before the time points, a longer duration of cART was associated with an increased eGFR in all three models; after the time points, there was an inverse association between eGFR and Nonlinear progression of renal function over time.
To assess whether renal dysfunction progression consists with the nonlinear trajectory of eGFR, we introduced two different outcome definitions, namely, RRF and RKFD (see Fig 4) There was no increased risk of RKFD among non-TDF users who received cART for 2.15 years or more, nor among TDF users who received cART for less than 1.40 years. But, each additional one year of TDF exposure was associated with a 78% (95% CI: -91%, -49%) decreased risk of RKFD from 1.40 to 2.30 years, and a nearly 3-fold (95% CI: 1.08, 7.27) increased risk of RKFD for those on TDF for more than 2.30 years. Similar trends were observed in PSM data (see Table 4).

Sensitivity Analyses
Two sensitivity analyses, one conducted with imputed datasets and the other with patients not using PIs, indicated these results were robust (see Additional file 1: S2 to S5).

Discussion
This was the first study, to our knowledge, to investigate whether eGFR decline follows a nonlinear trajectory as renal dysfunction occurs, in HIV-1-positive patients initiating cART with normal eGFR. We present evidence from two analyses (the piecewise linear or logistic regression model) that the traditional assumption of a steady, linear decline does not apply to HIV-1 infected patients, especially those on TDF therapies. Our results showed that these patients experienced periods of acceleration or deceleration of kidney function decline. Analyses over nonlinear patterns speak to the true nature of the exposure-outcome relationships.
The comparison of one-linear and piecewise linear models suggested that nonlinear trajectory of eGFR was more accurate than a single linear process (log likelihood ratio test: P< 0.001 for all). When a single slope was fitted to the data, eGFR decline was either over-or under-estimated during partial period of cART. Intriguingly, nonlinear trajectories accurately depicted the periods of acceleration or deceleration of renal function decline, especially in TDF users who had an obvious heterogeneity in eGFR over time. This acceleration or deceleration, which was quantified by the piecewise linear mixed effects model, could be clearly identified from the data and smooth curves (see Table 3 and Fig.2). As illustrated for TDF users in model 2 (Fig.2), there was an increase of eGFR for intermediate cART durations (1.40-2.30 years), comparing markedly with the significant decline of eGFR either for short (<1.40 years) or long cART durations (>2.30 years). Certainly, these findings were similar in model 1 and model 3.
As expected, effects of nonlinearity of eGFR on renal dysfunction progression were well supported by the results of RRF and RKFD. In particular, the trends over time of RRF were completely consistent with nonlinear changes of eGFR (see Table 4). This finding was also robust, based on a range of sensitivity analyses.
Among TDF users, during the increasing period (1.40-2.30 years) of eGFR, the incidences of both outcomes, especially RKFD definitely declined (suggesting a recovery of renal function), even though TDF continued. This is consistent with previous studies suggested an overall limited effect of TDF on renal function decline. For example, a meta-analysis that compared ART regimens with or without TDF demonstrated a mean difference in eGFR of only 3.92 mL/min/1.73 m 2 on a short-term follow-up [9].
Interestingly, a cohort study reported the cumulative decline of eGFR attributable to TDF was 3.05, 4.05 and 2.42 (mL/min/1.73 m 2 ) at year 1, 2, 3, respectively; this indicates that the eGFR decline attributable to TDF was lower 3 years after than that of before, suggesting a partial eGFR recovery from years 2 to 3 [20]. However, specific time points for renal function recovery are difficult to obtain by their one-linear analysis of eGFR.
We also found that continuous TDF exposure inevitably led to renal impairment in a substantial population. The incidence of RRF -but not the severe RKFD -increased during the initial use of TDF, incidences of both outcomes increased significantly later, suggesting that persistent TDF exposure can lead to cumulative and irreversible renal impairment, even in those with a normal baseline renal function. This was in agreement with that of the prospective international cohort study published recently, the increased incidence of CKD per year of exposure to TDF was initially small (14%; 95% CI: 10%, 19%), yet doubled for a treatment period of 5 years [2]. Regrettably, the authors used also the conventional linear analysis to address this issue, thereby the nonlinear trajectories of eGFR progression, if exist, remains unknown. As suggested by studies from CKD cohorts, linear regression methods do not exactly estimate kidney function trajectories [16], considering the big heterogeneity with respect to kidney function, dropout and number of kidney function estimates [31]. Nonlinear statistical methods, such as piecewise-linear mixed effects model [15], are able to better characterize the different profiles of renal function progression, as well as to investigate specific risk factors associated with each profile [14,16]. Therefore, our study provides a new avenue for this difficult task, at least in HIV patients with normal renal function, although having less-strength data. Future extravalidation with prospective international cohort like D:A:D Study would benefit a lot to characterize the real trajectories of eGFR progression, as well as the potential time window to salvage renal function and to investigate the underlying mechanisms of TDF related nephrotoxicity.
This present study has several implications for our understanding of renal dysfunction progression in HIV-1 infected patients during cART with initial normal renal function. First, periods of slight increasing eGFR followed by periods of eGFR decline and increasing risk of adverse events in non-TDF users suggesting that irrespective of the cART regimen (with or without TDF), loss of renal function to some extent seems inevitable following prolonged use of these drugs, especially after 2 years exposure or more. Screening frequencies on renal function should be planned according to this finding. Second, for TDF users, periods of rapid eGFR decline followed by periods of eGFR improvement, indicating that eGFR decline may sometimes be ameliorated over a given extended periods. One should aware of this early loss of renal function may not reflect future loss of renal function. The S-shaped nonlinear trajectory of eGFR may also open new avenues of diagnostic and treatment options so as to delay the progression of renal impairment among these long-term users of TDF.
This study has several strengths. First, the research has longitudinal data for up to 7 years of followup and regular eGFR assessments every 3 months for characterizing nonlinear trajectories of eGFR during cART. Second, by using of PSM, we were able to reduce confounding bias and balance the baseline characteristics of TDF exposure and non-exposure group. The results of this emulation of a randomized controlled trial were similar with model 1 and model 2, suggesting that our findings were robust. Third, the time points suggested by our study were determined by a range of powerful statistical analyses (Wald test, piecewise linear mixed effects model along with maximum likelihood model and recursion method), together with two robust sensitivity analyses, thus is more accurate and powerful than the traditional paradigm based on clinical experience [2,13,20].
Our study has several limitations. First, the inherent shortcomings belong to retrospective observational single-center study, small sample size and short term follow-up make it difficult to address the causality between TDF and CKD and reach a firm conclusion, the powerful statistical analysis thus is a trade-off to minimize these biases and confounding. Second, the patients in this study come exclusively from China, the findings may not simply apply to other populations and thus further validations from different races are warranted. Third, nonlinear trajectory of eGFR progression in patients complicated with CKD at baseline initiating TDF needs further investigation, after all, an interesting curve has already been identified by our population characterized by normal renal function. Fourth, this study did not investigate the predictive factors that may contribute to nonlinearity patterns of renal function, as well as TDF induced nephrotoxicity other than glomerular filtration function. All above limitations require further study to be overcome, nonetheless, our primary results provided moderate yet important illumination for this topic.
In conclusion, the present study suggests that renal function progression exists heterogeneity inHIVpositive patients with a normal eGFR initiating ART in Asians. There are significant differences in renal function trajectories between TDF and non-TDF therapy. Continuous TDF exposure inevitably led to renal impairment in a substantial population, but the changes of eGFR was inconsistent over time. An interesting reverse S-shaped nonlinear trajectory, the transient yet definitely recovery of renal impairment about 1.4 years after TDF initiation, do exist and could be helpful for the management of HIV-1-infected patients on TDF, which may provide a new avenue for patients care to improve the quality of their life.

Additional Files
Additional file 1: Table S1. Propensity score parameter list. Table S2. Predicted eGFR change rates in the piecewise linear mixed effects model in multiple imputation sample. Table S3. Association of antiretroviral exposure (in different time ranges) with risk of renal impairment outcomes in multiple imputation sample. Table S4. Predicted eGFR change rates among patients without receiving protease inhibitors in the piecewise linear mixed effects model. Table S5. Association of antiretroviral exposure (in different time ranges) with risk of renal impairment outcomes among patients without receiving protease inhibitors. Figure S1. Nonlinear trajectory of eGFR among HIV-1-infected patients with or without TDF in unadjusted model. Figure

Acknowledgement
We thank all study participants and staff of all participating sites.

Funding
The study did not need and was not funded.

Availability of data and materials
The data set used for this manuscript will be available from the corresponding author upon reasonable request.

Authors' Contributors
LF and XYH conceived, designed, and organised the study, interpreted the results, and drafted the manuscript. LF and XYH contributed equally to this manuscript. ZHQ analysed the data. YZX, CC, BR helped supervise the study, and revised the manuscript. The other authors contributed to collecting the data on site.

Ethics approval and consent to participate
This study was approved by the Institutional Review Board of Xixi Hospital. All data were anonymized to comply with the provisions of personal data protection legislation. Due to the retrospective nature of this study and due the fact that only historical medical data were collected, written informed consent was not required.    Abbreviations: Exp(β), the rate of change in eGFR per year, obtained with the interaction term between TDF using status and time since cART initiation.   Nonlinear trajectory of eGFR among HIV-1-infected patients with or without TDF.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.