Longitudinal Progression of Estimated GFR in HIV-1-Infected Patients with Normal Renal Function on Tenofovir-Based Therapy in China.

PURPOSE: Estimated glomerular filtration rate (eGFR) decline in HIV-1-infected patients exposure to tenofovir disoproxil fumarate (TDF) has been widely assessed using linear models, but nonlinear assumption is not well validated. We constructed a retrospective cohort study to assess whether eGFR decline follows nonlinearity during antiviral therapy. PATIENTS AND METHODS: We examined 823 (299 of TDF users and 524 of non-TDF users) treatment-naïve HIV-1-infected participants (age ≥ 17 years, initial eGFR ≥ 90 mL/min/1.73m2). Estimated GFR trajectories were compared by one-linear and piecewise-linear mixed effects models, before and after propensity score matching, respectively. Whether the incidence of renal dysfunction (reduced renal function [RRF], eGFR < 90 mL/min/1.73 m2 and rapid kidney function decline [RKFD], eGFR > -3 mL/min/1.73 m2/year) follows nonlinearity was assessed by logistic regression.
RESULTS: The median follow-up time of this study was 10 (interquartile range, 2-20) months, during which 178 (21.6%) experienced RRF, and 451 (54.8%) experienced RKFD. The slopes (mL/min/1.73 m2/year) of eGFR were -5.31 (95% CI: -6.57, -4.06) before 1.40 years, 4.83 (95% CI: 1.38, 8.28) from years 1.40 to 2.30 and -3.71 (95% CI: -5.97, -1.45) after 2.30 years among TDF users. Within years 1.40-2.30, each year of TDF exposure was associated with a 78% decreased risk of RKFD (95% CI: -91%, -49%). In comparison, eGFR increased slightly at the initiation of antiviral therapy, declined after 2.15 years (-4.96; 95% CI: -5.76, -4.17) among non-TDF users. Such a progression nonlinear trajectory was missed on the assumption of one-linearity, whether in TDF or non-TDF users.
CONCLUSION: Over the piecewise mixed-effects analyses with the advantage of revealing the true nature of the exposure outcome relationships, an interesting reverse S-shaped relationship was observed. A routine screen based on nonlinearity could be more helpful for patient management.

Introduction

The widespread use of combination antiretroviral therapy (cART) has essentially improved the life expectancy of human immunodeficiency virus (HIV)-positive individuals.1 Tenofovir disoproxil fumarate (TDF), an inhibitor of nucleotide analogue reverse transcriptase, which is widely used in most countries around the world as a conventional component of cART for HIV treatment and is considered as the most cost-effective drug against HIV.2,3 In addition, TDF has been approved as part of a pre-exposure prophylaxis (PreP) to prevent the spread of HIV in those who are at high risk for contracting this virus.4 However, TDF is similar to adefovir and cidofovir, which possess potential nephrotoxicity, lifelong use of TDF can cause or exacerbate renal impairment,5,6 and more and more concerns have been raised on renal toxicity of TDF to improve patients’ quality of life during this drug exposure. Thus, accurate predictive analyses of renal function overtime will be helpful for the management of these patients.

Estimated glomerular filtration rate (eGFR) is a common indicator of renal function.7,8 Studies have consistently demonstrated that TDF is associated with a decline of eGFR and renal dysfunction in a subpopulation.9–14 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 approximately linear. The real trajectory of eGFR over time is however missed in these simplified models, thus hinders the optimization of TDF therapy based on 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.15–17 However, no studies from HIV-1-infected 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. We also assessed the incidence of renal dysfunction based on nonlinear changes in eGFR, by using a two-piecewise logistic regression model.

Patients and Methods

Study Population

This is a retrospective, observational cohort study conducted at the infectious diseases department at Xixi Hospital of Hangzhou (Zhejiang, Southeast China). All treatment-naïve HIV-1-positive patients with records of cART initiation between January 26, 2010 and December 31, 2015 were screened for eligibility. All data were anonymized to comply with the provisions of personal data protection legislation. Due to the retrospective nature of this study and all data were collected anonymously, written informed consent was not required. This study was approved by the Institutional Review Board of Xixi Hospital.

Data Collection and Inclusion Criteria

Data extracted from the medical records included demographic parameters, date of cART initiation, details of the cART regimens, route of HIV-1 transmission, comorbidities, laboratory variables (HIV-1 RNA viral load, CD4+ lymphocyte cell count, and serum creatinine [SCr]) at baseline, and SCr at 2 weeks, 1 month, 2 months, 3 months, and every 3 months thereafter until January 2017. Isotope dilution mass spectrometry traceable calibration method was used to standardize the measurement of SCr. Baseline was defined as the date of starting cART. Each enrolled patient was 17 years old or more, had a normal baseline eGFR, and had at least one additional eGFR measurement since January 2010. The flowchart is detailed in Figure 1.

Figure 1

Study flow diagram.

Quantitative Variables

The three-variable Modification of Diet in Renal Disease (MDRD) formula adjusted for Chinese populations was used to calculate the eGFR values, as the Chinese eGFR investigation collaboration recommend the use of MDRD equation for Chinese, rather than CKD-EPI.18–20

Combination ART was defined as the combined use of three or more ARVs from any drug class. Patients who took TDF alone or any TDF-containing regimen (TDF + lamivudine [3TC], or emtricitabine [FTC] + nevirapine [NVP], or efavirenz [EFV], or zidovudine [AZT]) were classified as TDF users. Patients exposed to any ARVs except TDF (AZT, or stavudine [d4T] + 3TC + NVP, or EFV) were classified as non-TDF users.

The two outcome definitions of this study were reduced renal function (RRF: eGFR ≥ 90 mL/min/1.73 m2 at baseline and eGFR < 90 mL/min/1.73 m2 during follow-up),21 and rapid kidney function decline (RKFD: with progression to CKD; eGFR decline > 3 mL/min/1.73 m2/year, estimated by least-squares regression).22

Statistical Analyses

Baseline characteristics were compared between TDF users and non-TDF users. Three models were used to analyze eGFR progression over time since ART initiation in each group (Table 1). 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 at baseline. Model 3 used propensity score matching (PSM) to reduce preexisting imbalances in the covariates and potential confounding,23,24 and a covariate was considered well balanced when the P value was more than 0.05 (Table 2), more technical details were as in additional .

Table 1

The Difference of Slopes Before and After Cutoff Times and Comparison of One-Linear and Piecewise-Linear Models

Without TDF		With TDF
Model 1
Comparison of slopes	Exp(β)a (95% CI)	Comparison of slopes	Exp(β)a (95% CI)
<2.55 y	−4.79 (−5.84, −3.74) <0.001	<1.40 y	−8.47 (−11.56, −5.37) <0.001
≥2.55 y		≥1.40 y, < 3.20 y
Comparison of models	Log likelihood ratio testb	Comparison of slopes
One-linear model	<0.001	≥1.40 y, <3.20 y	−9.22 (−12.52, −5.92) <0.001
Non-linear model		≥3.20 y
		Comparison of models	Log likelihood ratio testb
		One-linear model	<0.001
		Non-linear model
Model 2
Comparison of slopes	Exp(β)a (95% CI)	Comparison of slopes	Exp(β)a (95% CI)
<2.15 y	−5.43 (−6.47, −4.40) <0.001	<1.40 y	−10.14 (−14.44, −5.85) <0.001
≥2.15 y		≥1.40 y, <2.30 y
Comparison of models	Log likelihood ratio testb	Comparison of slopes
One-linear model	<0.001	≥1.40 y, <2.30 y	−8.54 (−12.67, −4.41) <0.0001
Non-linear model		≥2.30 y
		Comparison of models	Log likelihood ratio testb
		One-linear model	<0.001
		Non-linear model
Model 3
Comparison of slopes	Exp(β)a (95% CI)	Comparison of slopes	Exp(β)a (95% CI)
<2.15 y	−4.28 (−6.24, −2.33) <0.001	<1.30 y	−7.09 (−13.99, −0.20) 0.044
≥2.15 y		≥1.30 y, <2.10 y
Comparison of models	Log likelihood ratio testb	Comparison of slopes
One-linear model	<0.001	≥1.30 y, <2.10 y	−8.82 (−14.89, −2.76) 0.004
Non-linear model		≥2.10 y
		Comparison of models	Log likelihood ratio testb
		One-linear model	<0.001
		Non-linear model

Notes: aExp(β) represents the difference of segmented slopes (mL/min/1.73 m2/year), along with a p value from Wald test. bLog likelihood ratio test was used to compare one-linear regression model with two-piecewise regression model, below 0.05 indicates two-piecewise regression model was a better fit to the data than the one-linear model that assumed a single slope across the entire period of observation. Model 1: unadjusted for any variables at baseline. Model 2: 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 at baseline. Model 3: propensity score-matched sample.

Table 2

Characteristics at Cohort Entry Stratified by Tenofovir Disoproxil Fumarate Before and After Propensity Score Matching

Characteristics	Before Matching			After Matching
	Without TDF	With TDF	P value	Without TDF	With TDF	P value
Overall	(n=524,63.7%)	(n=299,36.3%)		(n=260,66.7%)	(n=130,33.3%)
Age (years)	27 (24–32)	30 (25–36)	<0.001	27 (25–32)	27 (25–33)	0.638
Female	20 (3.8%)	18 (6.0%)	0.147	5 (1.9%)	5 (3.8%)	0.428
Weight (kg)	63 (57–70)	63 (56–67)	0.185	62 (57–70)	63 (58–68)	0.810
Height (cm)	172 (169–175)	172 (169–175)	0.546	172 (170–175)	172 (170–175)	0.790
BMI (kg/m2)	21.1 (19.5–23.1)	21.0 (19.4–22.7)	0.240	21.0 (19.4–23.1)	21.2 (19.5–22.9)	0.637
CD4 (cells/μL)	323 (246–423)	247 (117–359)	<0.001	326 (262–420)	335 (246–414)	0.988
Triglycerides (mmol/L)	1.1 (0.8–1.7)	1.3 (0.9–1.7)	0.715	1.2 (0.8–1.7)	1.2 (0.8–1.7)	0.739
Total cholesterol (mmol/L)	4.0 (3.5–4.5)	3.8 (3.3–4.4)	0.006	4.0 (3.5–4.5)	3.9 (3.5–4.4)	0.676
eGFR (mL/min per1.73m2)	111 (102–121)	112 (103–126)	0.426	112 (101–122)	111 (103–120)	0.767
Dyslipidemia	72 (13.7%)	42 (14.0%)	0.879	34 (13.1%)	17 (13.1%)	1.000
Risk Factors			0.005			0.689
Homosexual	413 (78.8%)	203 (67.9%)		205 (78.8%)	100 (76.9%)
Heterosexual	69 (13.2%)	60 (20.1%)		31 (11.9%)	20 (15.4%)
Injection drug user	1 (0.2%)	0 (0.0%)		1 (0.4%)	0 (0.0%)
Other	41 (7.8%)	36 (12.0%)		23 (8.8%)	10 (7.7%)
WHO stage III/IV	84 (16.0%)	103 (34.4%)	<0.001	31 (11.9%)	13 (10.0%)	0.692
Hepatitis B Status			<0.001			NA
Positive	9 (1.7%)	43 (14.4%)		0 (0.0%)	0 (0.0%)
Negative	476 (90.8%)	236 (78.9%)		260 (100.0%)	130 (100.0%)
Unknown	39 (7.4%)	20 (6.7%)		0 (0.0%)	0 (0.0%)
Hepatitis C Status			0.003			NA
Positive	3 (0.6%)	11 (3.7%)		0 (0.0%)	0 (0.0%)
Negative	472 (90.1%)	255 (85.3%)		260 (100.0%)	130 (100.0%)
Unknown	49 (9.4%)	33 (11.0%)		0 (0.0%)	0 (0.0%)
Anaemia	11 (2.1%)	40 (13.4%)	<0.001	3 (1.2%)	0 (0.0%)	0.539
Diabetes	16 (3.1%)	10 (3.3%)	0.818	8 (3.1%)	2 (1.5%)	0.571
Viral Load (Copies per mL)			0.012			0.712
<400	39 (7.4%)	17 (5.7%)		26 (10.0%)	9 (6.9%)
≥400, <10,000	127 (24.2%)	45 (15.1%)		65 (25.0%)	29 (22.3%)
≥10,000, <100,000	92 (17.6%)	56 (18.7%)		49 (18.8%)	23 (17.7%)
≥100,000	29 (5.5%)	15 (5.0%)		15 (5.8%)	8 (6.2%)
Unknown	237 (45.2%)	166 (55.5%)		105 (40.4%)	61 (46.9%)
Protease inhibitors	8 (1.5%)	34 (11.4%)	<0.001	2 (0.1%)	1 (0.1%)	1.000

Notes: Data are n (%) or median (IQR) unless otherwise indicated. Baseline was defined as the date of starting antiretroviral therapy on or after January 2010. After matching, P value > 0.05 indicates a relatively small baseline imbalance between TDF and non-TDF users. Diabetes and dyslipidemia defined by the diagnosis or related medication. Anemia was defined as hemoglobin <12.0 g/dL in women and <13.0 g/dL in men. Coinfection with hepatitis B defined by positive hepatitis B surface antigen, coinfection with hepatitis C defined by positive HCV viral load.

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 (Table 1). 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 (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 (Table 4). All multivariate regression models were adjusted for the covariates used in Model 2.

Table 3

Predicted eGFR Change Rates in the Piecewise-Linear Mixed Effects Model

Without TDF			With TDF
	Exp(β) (95% CI)	P value		Exp(β) (95% CI)	P value
Model 1 (n=823 Patients, 11,422 Measurements)
Time as linear trend	−1.29 (−1.58, −1.00)	<0.001	Time as linear trend	−1.46 (−1.94, −0.98)	<0.001
Fitted Groups			Fitted Groups
<2.55 y (n=6098 measurements)	0.74 (0.21, 1.28)	0.006	<1.40 y (n=3172 measurements)	−4.73 (−6.09, −3.37)	<0.001
≥2.55 y (n=900 measurements)	−4.04 (−4.72, −3.37)	<0.001	≥1.40 y, <3.20 y (n=996 measurements)	3.74 (1.64, 5.84)	0.004
–	–	–	≥3.20 y (n=256 measurements)	−5.48 (−8.03, −2.93)	<0.001
Model 2 (n=707 Patients, 8507 Measurements)
Time as linear trend	−1.20 (−1.54, −0.85)	<0.001	Time as linear trend	−2.56 (−3.19, −1.94)	<0.001
Fitted Groups			Fitted Groups
<2.15 y (n=4857 measurements)	0.47 (0.00, 0.94)	0.049	<1.40 y (n=2395 measurements)	−5.31 (−6.57, −4.06)	<0.001
≥2.15 y (n=492 measurements)	−4.96 (−5.76, −4.17)	<0.001	≥1.40 y, <2.30 y (n=551 measurements)	4.83 (1.38, 8.28)	0.006
–	–	–	≥2.30 y (n=212 measurements)	−3.71 (−5.97, −1.45)	0.001
Model 3 (n=390 Patients, 4663 Measurements)
Time as linear trend	−0.47 (−1.09, 0.15)	0.139	Time as linear trend	−1.77 (−2.60, −0.94)	<0.001
Fitted Groups			Fitted Groups
<2.15 y (n=2794 measurements)	0.77 (−0.07, 1.60)	0.072	<1.30 y (n=1124 measurements)	−2.78 (−4.73, −0.83)	0.005
≥2.15 y (n=271 measurements)	−3.51 (−5.04, −1.99)	<0.001	≥1.30 y, <2.10 y (n=306 measurements)	4.31 (−1.28, 9.90)	0.131
–	–	–	≥2.10 y (n=168 measurements)	−4.51 (−6.86, −2.17)	<0.001

Notes: Exp(β), the rate of change in eGFR (mL/min/1.73 m2) per year, obtained with the interaction term between TDF using status and time since cART initiation. Model 1: unadjusted for any variables at baseline. Model 2: 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 at baseline. Model 3: propensity score matched sample.

Table 4

Association of Antiretroviral Exposure (in Different Time Ranges) with Risk of Renal Impairment Outcomes

Unmatched Samplea
Without TDF			With TDF
	OR (95% CI)	P value		OR (95% CI)	P value
Reduced Kidney Functionb
Time as linear trend	1.67 (1.42, 1.98)	<0.001	Time as linear trend	1.80 (1.54, 2.09)	<0.001
Fitted Groups			Fitted Groups
<2.15 y	1.33 (0.97, 1.81)	0.074	<1.40 y	3.33 (2.34, 4.75)	<0.001
≥2.15 y	2.05 (1.54, 2.71)	<0.001	≥1.40 y, <2.30 y	0.59 (0.25, 1.39)	0.229
–	–	–	≥2.30 y	1.58 (1.03, 2.43)	0.035
Rapid Kidney Function Declinec
Time as linear trend	0.91 (0.84, 0.98)	0.020	Time as linear trend	1.05 (0.93, 1.18)	0.418
Fitted Groups			Fitted Groups
<2.15 y	0.89 (0.80, 1.00)	0.048	<1.40 y	1.07 (0.87, 1.32)	0.512
≥2.15 y	0.94 (0.77, 1.14)	0.524	≥1.40 y, <2.30 y	0.22 (0.09, 0.51)	<0.001
–	–	–	≥2.30 y	2.80 (1.08, 7.27)	0.034
Matched Sampled
Without TDF			With TDF
	OR (95% CI)	P value		OR (95% CI)	P value
Reduced Kidney Functionb
Time as linear trend	1.38 (1.12, 1.70)	0.003	Time as linear trend	1.49 (1.25, 1.78)	<0.001
Fitted Groups			Fitted Groups
<2.15 y	1.23 (0.84, 1.79)	0.287	<1.30 y	2.62 (1.50, 4.59)	<0.001
≥2.15 y	1.54 (1.08, 2.20)	0.017	≥1.30 y, <2.10 y	0.56 (0.14, 2.33)	0.429
–	–	–	≥2.10 y	1.34 (0.90, 1.99)	0.152
Rapid Kidney Function Declinec
Time as linear trend	1.01 (0.92, 1.11)	0.834	Time as linear trend	1.15 (0.99, 1.34)	0.064
Fitted Groups			Fitted Groups
<2.15 y	0.94 (0.82, 1.08)	0.396	<1.30 y	1.19 (0.87, 1.62)	0.275
≥2.15 y	1.17 (0.94, 1.45)	0.171	≥1.30 y, <2.10 y	0.19 (0.07, 0.56)	0.002
–	–	–	≥2.10 y	12.43 (0.78, 197.43)	0.074

Notes: aRepresents the model 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 at baseline. bReduced kidney function was defined as the development of an eGFR< 90mL/min/1.73m2 during follow-up among patients who had an eGFR greater than or equal to 90 mL/min/1.73m2 at baseline. cRapid kidney function decline was defined as an annual decline of 3 mL/min/1.73m2 or more. dRepresents the propensity score-matched model.

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.25,26 In addition, to reduce bias caused by exclusion of individuals with any missing data at baseline, five imputed datasets (established by multiple imputation with chained equations) were developed and run separately, and the results were combined using Rubin’s method ( and ).27,28 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 ( and ).29–31

All analyses were performed using the R software, version 3.3.1 (). A result was considered statistically significant when the two-tailed P value was below 0.05.

Results

Patient Selection and Propensity Score Matching

As shown in the flowchart (Figure 1), a total of 1065 patients were screened and 823 patients were eligible for participation, 299 of whom (36.3%) started a TDF-containing cART. Table 2 shows the baseline characteristics of TDF users and non-TDF users before and after PSM. After matching, there were 130 (33.3%) patients in the TDF group, and all baseline variables were well balanced (P > 0.05 for all).

The median age was 30 years among TDF users, and 27 years among non-TDF users. Most enrolled patients were male and were infected via male-male sex. Of 823 patients, 178 (21.6%) experienced RRF, and 451 (54.8%) experienced RKFD over a median follow-up of 10 (interquartile range [IQR], 2–20; maximum 90) months. In TDF users, 97 experienced (32.4%) RRF, and 183 (61.2%) experienced RKFD. There were 4424 eGFR measurements for TDF users. For each group, there was a median of seven eGFR measurements per person (IQR, 3–11) and the median interval between eGFR measurements was 90 (IQR, 30–90) days.

Main Analyses

Comparison of One-Linear and Piecewise-Linear Mixed Effects Models

We compared eGFR trajectories by one-linear and piecewise-linear models (Table 1), 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 the 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 (Table 1, model 1), 2.15 years (Table 1, model 2), and 2.15 years (Table 1, model 3). The difference of eGFR slopes was −4.79 (−5.84, −3.74), −5.43 (−6.47, −4.40) and −4.28 (−6.24, −2.33), respectively.

For TDF users, the time points on the nonlinear trajectory of eGFR were significantly different from that of non-TDF users. For example, the time points were 1.40 years and 2.30 years in model 2. The difference of eGFR slopes was −10.14 (−14.44, −5.85) at 1.40 years and −8.54 (−12.67, −4.41) at 2.30 years. Similar results were obtained in models 1 and 3 (Table 1).

The Relationship Between eGFR and Duration of cART

The eGFR changed over time in both groups (Figure 2, and ). 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 () and model 2 (Figure 2B).

Figure 2

Nonlinear trajectory of eGFR among HIV-1-infected patients with or without TDF.

Notes: Nonlinear eGFR changes over time can be approximated with a piecewise-linear mixed effects model. (A)  and (B) show the adjusted smooth fit of eGFR data. (C) and (D) show the fit from the adjusted one linear and adjusted piecewise-linear mixed effects models. Models adjusted for age, sex, weight, height, 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 at baseline.

Average Changes in eGFR Over Time on Different cART Duration Among TDF or Non-TDF Users

Table 3 shows 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. The exp(β) was −5.31 (95% CI: −6.57, −4.06) for cART less than 1.40 years and −3.71 (95% CI: −5.97, −1.45) for 2.30 years or more. However, the exp(β) was reverse, 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 (Table 3).

For non-TDF users, before the time points, a longer duration of cART was associated with a slight increased eGFR in all three models; after the time points, there was an inverse association between eGFR and duration of cART (Table 3).

Nonlinear Progression of Renal Function Over Time

Two outcome definitions, RRF and RKFD, were used to assess whether renal dysfunction progression consists of the nonlinear trajectory of eGFR (Table 4). For patients without TDF exposure, use of cART for 2.15 years or more, the risk of RRF increased steadily to 2.05 per year (95% CI: 1.54, 2.71). For patients using TDF, there was an increased risk of RRF for those using cART less than 1.40 years (adjusted odds ratio [aOR]: 3.33 per year; 95% CI: 2.34, 4.75) and for those using cART for 2.30 years or more (aOR: 1.58 per year; 95% CI: 1.03, 2.43). However, those using TDF for 1.40 to 2.30 years had a decreased risk of RRF (41% decrease per year; 95% CI: −75%, 39%).

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 1 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 threefold (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 (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 ().

Discussion

This was the first study, to our knowledge, to investigate whether eGFR progression follows a nonlinear trajectory in HIV-1-infected patients initiating cART with normal eGFR. We present evidence from two analyses (the piecewise-linear and logistic regression model) that the traditional assumption of a steady, linear decline does not apply to HIV-1 infected patients on treatment, especially those on TDF-based therapies. Our results showed that these patients experienced periods of acceleration or deceleration of kidney function decline. Analyses over nonlinear patterns seemly speak to the true nature of the exposure–outcome relationships.

The comparison of one-linear and piecewise-linear models suggested that the 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 the 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 (Table 3 and Figure 2). As illustrated for TDF users in model 2 (Table 3), 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, the 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 (Table 4). This finding was also robust enough based on a range of sensitivity analyses. This phenomenon can not be explained explicitly thus far.32 A speculation of far from mature is that TDF, as a well-known nephrotoxic antiretroviral, causes a rapid stress in renal tubular at the beginning exposure followed by a transient recovery possibly from the self-repairing mechanisms of kidney; then, an inevitable damage occurs if beyond the ability of self repairment over time.33

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.10,21 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 m2 on a short-term follow-up.10 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 m2) 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.21 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. TDF-induced nephrotoxicity was reported in 0.5–45% of HIV-positive patients.6 The wide range of prevalence is attributed to different populations and definitions of TDF-induced nephrotoxicity and duration of follow-up. Renal function assessment and monitoring at baseline and during TDF treatment is the main approach of prevention of TDF-induced nephrotoxicity. But how to monitor appropriately is a challenging issue in daily practice. 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.5 Regrettably, the authors used also the conventional linear analysis to address this issue, thereby the nonlinear trajectories of eGFR progression, if exist, remain unknown. As suggested by studies from CKD cohorts, linear regression methods do not exactly estimate kidney function trajectories,17 considering the big heterogeneity with respect to kidney function, dropout and number of kidney function estimates.34 Nonlinear statistical methods, such as piecewise-linear mixed effects model,16 are able to better characterize the different profiles of renal function progression, as well as to investigate specific risk factors associated with each profile.15,17 Therefore, our study provides a new avenue for this difficult task, at least in HIV patients with normal renal function. Future external validation 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 period. One should be aware of early loss of renal function may not reflect permanent 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 follow-up and regular eGFR assessments every 3 months for characterizing nonlinear trajectories of eGFR during cART. Second, by using 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.5,14,21

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 came exclusively from China and mainly with no history of drug abusing which is a risk factor for HIV, 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 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.

Conclusion

The present study suggests that renal function progression exists heterogeneity in HIV-infected patients with a normal eGFR initiating ART in Chinese. 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 in eGFR were inconsistent over time. Analyses assuming nonlinear patterns over piecewise mixed effects models speak to the true nature of the exposure–outcome relationships in this scenario. 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.