Prediction of Tacrolimus Dose/Weight-Adjusted Trough Concentration in Pediatric Refractory Nephrotic Syndrome: A Machine Learning Approach

Introduction

Nephrotic syndrome (NS) is a common chronic disease in children, with clinical manifestations of proteinuria, hypoproteinemia, hyperlipidemia, and edema. Tacrolimus (TAC, also known as FK506), a calcineurin inhibitor, is the first-line immunosuppressant for patients with refractory NS. The pharmacokinetics (PK) characteristics of TAC are its low and variable oral bioavailability (4~89%) and extensive binding to red blood cells. TAC is a substrate of cytochrome P450 (CYP) 3A and P-glycoprotein (P-gp), and it is mainly excreted as metabolites (approximately 99%) via biliary route.¹ However, it is difficult to develop an individualized TAC dosing regimen for patients due to the narrow therapeutic window and large interindividual pharmacokinetic variability. Therapeutic drug monitoring (TDM) is required to individually adjust the TAC dose. Therefore, identifying the key factors that affect TAC PK, developing a prediction model, and achieving individualized medication is of great significance.

Interindividual differences in drug disposition and response are mainly related to genetic polymorphisms in drug-metabolizing enzymes, transporters, receptors, and other drug targets.² TAC is mainly metabolized by CYP3A, and it is a substrate of P-gp.¹ Most previous studies on TAC focused on the genetic polymorphisms of TAC PK pathway-related genes (CYP3A, MDR1, etc.) and transcriptional factor coding genes (POR, NR1I2, etc.).^3–5 And CYP3A5*3 (rs776746) is the only confirmed single nucleotide polymorphisms (SNPs) affecting TAC PK.⁶ However, CYP3A5*3 cannot fully explain the individual variations in TAC PK. On the other hand, due to the high plasma protein binding rate of TAC, physiologic function of the kidneys may also contribute to TAC interindividual variability. Alpha-actinin-4 (encoded by ACTN4) is critical for renal filtration barrier, and its mutations may cause kidney damage.^7,8 And previous studies on MYH9 (encodes non-muscle myosin II) have shown that changes in MYH9 expression may damage the glomerular filtration barrier.^9,10 Therefore, the genetic polymorphisms on ACTN4 and MYH9 may also affect the concentration of TAC.¹¹ Furthermore, most previous studies on individualized drug therapy of TAC were investigated in transplant patients, and studies on pediatric patients with NS are limited.^3–5,12 Thus, it is necessary to investigate the impact of those genetic polymorphisms on TAC concentration in refractory NS children, to clarify the reasons for the large individual variations in TAC PK.

Developing a prediction model could provide doctors with reliable and convenient drug administration decision support. However, only several studies have developed TAC population pharmacokinetic (PPK) models for pediatric refractory NS. Wang et al developed a TAC PPK model based on age, cystatin-C, and the daily dose of TAC,¹³ whereas Hao et al found that body weight and CYP3A5 significantly affected TAC PK.¹⁴ However, Wang et al did not investigate the effect of genetic polymorphisms and Hao et al have only investigated the effect of CYP3A5, which did not fully explain the individual differences in TAC PK. These limitations reduce the accuracy of the models and limit their clinical application. Thus, a simple and accurate TAC concentration prediction model based on pediatric NS is urgently needed to provide a reference for clinicians to make a drug delivery regimen. In recent years, the prediction models based on machine learning have excellent performance in the diagnosis, treatment, and prognosis of diseases.^15–17 Tang et al have developed machine-learning models to predict TAC stable dose in renal transplant recipients, and their performance is better than that of traditional statistical methods.¹⁸ However, there has been no research on the prediction of the TAC concentration in NS patients using machine learning.

Therefore, the purpose of this study is to investigate the effects of clinical features and genetic polymorphisms on TAC dose/weight-adjusted trough concentration (C₀/D) in pediatric patients with refractory NS, and then use machine learning algorithms to develop a simple and accurate TAC C₀/D prediction model to provide a powerful tool for the individualized drug therapy of TAC in pediatric patients with refractory NS.

Materials and Methods

Study Design and Patients

We retrospectively included patients who visited Guangzhou Women and Children’s Medical Center from June 2013 to December 2018 as a derivation cohort, which was used to develop the prediction model. We also prospectively collected some patients as a verification cohort to test the validity of the model.

Patients under 16 years of age, who were diagnosed with refractory NS (including steroid-dependent NS, steroid-resistant NS, and frequent relapses NS) and took TAC orally for at least 3 months were included in the study. Patients with steroid-sensitive NS and with poor compliance were excluded. Those receiving medication affecting TAC blood levels (verapamil, ketoconazole, rifampicin, Wuzhi tablet, etc.) were also excluded. All patients received TAC (Prograg™, Astellas, Killorglin, Ireland) and prednisone (Guangdong Huanan Pharmacy Ltd., Dongguan, China). The initial dose of TAC is 0.10–0.15 mg/kg twice daily, and doses were adjusted according to TDM.

This study, designed in accordance with the requirements of the Declaration of Helsinki, was approved by the Ethics Committee of Guangzhou Women and Children Medical Center (No. 201509), and it was registered at ClinicalTrials.gov (NCT02602873). Informed consent was signed by the patients’ parents or guardians. Additionally, the data used in this study were anonymous.

Clinical and Genetic Variables

We collected 82 clinical variables, including demographic data (eg age, gender), hematological characteristics (eg trough concentration, albumin level), urine features (eg urinary protein, urine red blood cells), medication data (eg Wuzhi tablet, ketoconazole), etc. (see Table S1 in Supplementary). Genetic variables, SNPs, were selected based on two criteria: (1) affecting TAC pharmacokinetics (both confirmed and unconfirmed); (2) the minor allele frequency (MAF) was ≥5% for Chinese subjects. As a result, 244 SNPs were selected (see Table S2 in Supplementary).

Concentration Assay

Blood samples for TAC concentration detection and genotyping were collected before taking TAC. The enzyme-multiplied immunoassay technique (Viva-E, Siemens, Germany) was used to determine the whole blood concentration of TAC. The calibration curve range was 2.5–30.0 ng/mL.

Genotyping

DNA for SNP analysis was isolated using the Genome TIANGEN Blood DNA Extraction Kit (DP348, Beijing, China). Polymerase chain reaction-restriction fragment length polymorphism methods were used to determine the CYP3A5 rs776746, MDR1 rs1045642, MDR1 rs1128503, MDR1 rs2032582, and HSD11B1 rs846908, and Agena Bioscience MassARRAY^® system (Agena Bioscience, San Diego, CA, USA) was used to determine other SNPs. The Hardy-Weinberg equilibrium test was performed using the chi-square test.

Machine Learning

It is well known that CYP3A5 genotype significantly affects the concentration of TAC. In order to exclude the interference of CYP3A5 genotype, we further divided the patients into CYP3A5 expresser group (CYP3A5*1/*1+CYP3A5*1/*3) and CYP3A5 nonexpresser group (CYP3A5*3/*3), then developed three prediction models separately. The main process of machine learning can be divided into three steps: (1) data preprocessing, (2) feature selection, (3) model generation and validation. The performance of feature selection and model generation was evaluated using five-fold stratified cross-validation. The flowchart of this work is shown in Figure 1. Machine learning was executed in Python 3.6.5 using Scikit-learn 0.19.1, and the main hyperparameters of each machine learning method are as shown in Table S3. GraphPad Prism 8.0.2 and CorelDRAW were used for drawing figures. P-value <0.05 was considered statistically significant.

Figure 1 The flowchart of model generation and validation.

Data Preprocessing

Variables with missing rates >10% were removed. The missing values of continuous and categorical variables were filled with the mean and mode values, respectively. Then, the categorical variables were transformed into dummy variables, while continuous variables were normalized by Z-score normalization.

Feature Selection

Univariate analysis was used to assess the association between each variable and TAC C₀/D. The variables with p≥0.05 were excluded. Then, data transformation was carried out on the original continuous variables to form other two kinds of data: min-max normalization and L2 normalization. Four algorithms, including eXtreme Gradient Boosting (XGBoost), Extremely Randomized Trees (ET), Random Forest (RF), and Gradient Boosting Decision Tree (GBDT) were selected to develop models based on the above four forms of data (min-max normalization, z-score normalization, L2 normalization, and original, respectively). The important contribution of each variable in 16 models was ranking by the median. Then starting from the most important variable, XGBoost with five-fold cross-validation was used to iteratively generate a new model by adding one variable at a time. The R² of the model is evaluated until the last variable is added, aiming to select the optimal features.

Model Generation and Validation

The Lasso regression, ET, GBDT, RF, and XGBoost algorithms were used to analyze the final feature set and generate prediction models. During model generation, patients were randomly divided into a training set and a test set according to the ratio of 8:2. Moreover, the five-fold cross-validation method was used for internal validation.

We further recruited 30 patients from January to October 2019 for external validation of the prediction model. Patients with missing data were excluded.

Results

Patient Characteristics

A total of 171 pediatric patients with refractory NS were included in this study, including 128 males and 43 females. The mean age of the children was 5.4±3.5 years, and the mean body weight was 19.5±8.3 kg. Patients were divided into CYP3A5 expresser group (CYP3A5*1/*1+ CYP3A5*1/*3) and CYP3A5 nonexpresser group (CYP3A5*3/*3) based on CYP3A5 genotype. There were 84 patients in the CYP3A5 expresser group, and 87 patients in the CYP3A5 nonexpresser group. The characteristics of patients are shown in Table 1.

Table 1 Demographics and Clinical Characteristics of All Pediatric Patients with Primary Nephrotic Syndrome

Feature Selection

We performed feature selection in the whole group, CYP3A5 expresser group, and CYP3A5 nonexpresser group separately. After univariate analysis, all variables with p<0.05 were further ranked of the median important contribution by ET, GBDT, RF, and XGBoost algorithms. The median important ranking of each variable is shown in Figure S1a/S2a/S3a in Supplementary.

The smaller the median value, the more important the variable. The process of stepwise forward modeling by XGBoost for the best feature set is shown in Figure S1b/S2b/S3b in Supplementary. From these figures, we see that the R² of the model based on the test set reaches maximization when the most appropriate variable is incrementally added. Therefore, these variables constitute the best feature set.

In the whole group, five variables constitute the final feature set, including ALB0 (p=0.032), CYP3A5*3/*3 (p<0.001), MYH9 rs2239781_CT (p=0.027), MYH9 rs2239781_CC (p=0.001), and CTLA4 rs4553808_AA (p=0.011). The higher the ALB0, the higher the TAC C₀/D. Compared with other genotypes, the patients with CYP3A5*3/*3 had a higher TAC C₀/D. The MYH9 rs2239781_CT showed a lower TAC C₀/D than other genotypes, whereas the MYH9 rs2239781_CC showed a higher TAC C₀/D than other genotypes. The TAC C₀/D of patients carrying CTLA4 rs4553808_AA was significantly lower than that of patients with other genotypes (see Figure 2).

Figure 2 The influence of ALB0 (A), CYP3A5 (B), MYH9 rs2239781 (C and D) and CTLA4 rs4553808 (E) on TAC C0/D in whole group (* p<0.05, ** p<0.01).

In the CYP3A5 expresser group, four variables were included in the final feature set, including AGE0 (p=0.021), MYH9 rs2239781_CC (p=0.005), MAP3K11 rs7946115_CG (p=0.011), and CTLA4 rs4553808_GG (p=0.017). The older the age, the higher the TAC C₀/D. The MYH9 rs2239781_CC and CTLA4 rs4553808_GG genotypes were associated with higher TAC C₀/D, whereas the MAP3K11 rs7946115_CG genotype was associated with lower TAC C₀/D (see Figure 3).

Figure 3 The influence of AGE0 (A), MAP3K11 rs7946115 (B), MYH9 rs2239781 (C) and CTLA4 rs4553808 (D) on TAC C0/D in CYP3A5 expressers (* p<0.05, ** p<0.01).

In the CYP3A5 nonexpresser group, nine variables were included, including GENDER (p=0.033), MYH9 rs2239781_CC (p=0.013), MYH9 rs4821478_AA (p=0.008), ACTN4 rs56113315_CT (p=0.004), INF2 rs1128880_TT (p=0.005), ACTN4 rs3745859_CT (p=0.004), ACTN4 rs62121818_CC (p=0.003), INF2 rs1128880_GT (p=0.040), and IL2RA rs12722489_CT (p=0.013). The heterozygote of MYH9, INF2, and IL2RA were associated with lower TAC C₀/D, whereas the heterozygote of ACTN4 was associated with higher TAC C₀/D. Gender also affected the TAC C₀/D, which was higher in males than in females (see Figure 4).

Figure 4 The influence of GENDER (A), MYH9 rs4821478 (B), MYH9 rs2239781 (C), IL2RA rs12722489 (D), INF2 rs1128880 (E and F), ACTN4 rs56113315 (G), ACTN4 rs62121818 (H) and ACTN4 rs3745859 (I) on TAC C0/D in CYP3A5 nonexpressers (* p<0.05, ** p<0.01).

Model Performance and Comparison

We used five algorithms (including XGBoost, ET, RF, GBDT, and Lasso regression) to develop and verify the prediction models in each group based on the final features. The performance of these models of the test set is shown in Table 2. The R² is between 0.207 and 0.444, and the error values are within the acceptable range. In the whole group, GBDT algorithm has the best performance with the highest R² of 0.444, and the least MSE, MAE and MedAE of 591.032, 20.782 and 18.980, respectively. In CYP3A5 expresser group, ET algorithm has the highest R² of 0.380 and the least MSE of 477.948. However, its MAE is relatively small, ranking second (the MAE of XGBoost algorithm is the least), and the MedAE of ET algorithm is relatively large (the MedAE of Lasso algorithm is the least). In CYP3A5 nonexpresser group, GBDT algorithm has the highest R² of 0.264 and the least MSE and MAE of 1839.459 and 31.257, respectively. In addition, the MedAE is relatively small, ranking second (the MedAE of ET algorithm is the least). The model performance is further analyzed from the perspective of precision and bias,¹⁹ which is shown in Table 2. In the whole group, GBDT has the best performance. In CYP3A5 expresser group, ET algorithm has the best performance. However, in CYP3A5 nonexpresser group, the precision of GDBT algorithm is the best, but its bias is relatively poor. The bias of the RF algorithm is the best, but its precision is relatively poor. The contribution of each variable to the outcome in the five models is described in Table S3 of Supplementary. The ALB0 has the greatest contribution to the GBDT model in the whole group. The AGE0 contributes the most important to the ET model in CYP3A5 expresser group, while the most important variable is ACTN4 rs62121818_CC in CYP3A5 nonexpresser group.

Table 2 The Performances Comparison of Five Models in Test Set in Whole Group, CYP3A5 Expresser Group, and CYP3A5 Nonexpresser Group

External Verification and Clinical Application

We also additionally collected 30 patients for external validation. After excluding patients lost to follow-up and patients with missing data, 25, 9, and 15 patients were finally enrolled into the whole group, CYP3A5 expresser group, and CYP3A5 nonexpresser group separately. The results of external validation are depicted in Figure 5. The predicted C₀/D are significantly associated with the observed C₀/D whose r values are 0.632, 0.664, and 0.532 (p=0.001, p=0.051, and p=0.041) separately. In clinical practice, we can input patients’ variables into the best model and predict the precision C₀/D for them (see Table 3).

Table 3 The Clinical Application of These Prediction Models. After Inputting the Variables into the Predicted Models, These Models Finally Output the Predicted C0/D Values

Figure 5 Observed and predicted C0/D in the whole group (A), the CYP3A5 expresser group (B), and the CYP3A5 nonexpresser group (C).

Discussion

To our knowledge, this is the first machine learning-based TAC C₀/D prediction model developed for pediatric patients with refractory NS that integrates clinical and genetic variables. Recently, machine learning algorithms have been widely used in the medical field and play an important role.^11,15,17,20 Tang et al established a TAC stable dose prediction model based on machine learning algorithms in renal transplant recipients,¹⁸ whereas we first used machine learning models to predict TAC C₀/D in NS patients. Our results show that models developed by advanced machine learning (RF, ET, GBDT, and XGBoost) perform better than those developed by Lasso regression, which verifies that advanced machine learning algorithms are generally superior to traditional algorithms in model fitting²¹ (see Table 2). In addition, machine learning can handle complex and high-dimensional variables, so it is particularly suitable for analyzing clinical, genetic, and other complex variables in this study. In the field of NS, Wang et al and Hao et al have established the TAC PPK models in children with NS.^13,14 However, the small sample size, complex sampling operations, and incomplete consideration of genetic variables lead to low model prediction performance, which limits the application of these models. In general, the improvement of analysis methods can help us find out the key factors that affect the individual variations in drug concentration more fully and accurately, and establish prediction models with better performance, so as to achieve personalized administration.

In the whole group, five variables were included in the prediction model (see Figure 2). The higher the blood ALB0, the higher the TAC C₀/D. TAC is a drug with a high plasma protein binding rate, while the pathophysiology of NS is characterized by hypoproteinemia and proteinuria. A higher blood ALB level mainly means less albumin is lost through the kidney, so we suspect that less TAC is lost with albumin, which leads to higher TAC C₀/D. CYP3A5 is the most critical metabolic enzyme of TAC, and CYP3A5 genotype is recognized as an important factor affecting TAC pharmacokinetics. Our results confirmed that patients with CYP3A5*3/*3 genotype had a significantly higher TAC C₀/D than patients with other genotypes, which is consistent with the results reported by other studies.^6,13,22 We also found that MYH9 rs2239781_CC showed a higher TAC C₀/D, whereas MYH9 rs2239781_CT led to a lower TAC C₀/D, which was the same as our previous study.¹¹ The MYH9 gene encodes the myosin IIA heavy chain, which is involved in various cell functions (such as cell division, migration, and adhesion). MYH9 mutations are associated with the onset of various diseases including chronic kidney disease.²³ Based on our results, we suspect that patients with MYH9 rs2239781_CT had more severe damage to kidney podocytes. Thus, more TAC was lost with plasma protein through the kidney, which resulted in a lower TAC C₀/D. However, further verification is needed. Cytotoxic T-lymphocyte-associated protein 4, encoded by the CTLA4 gene, is a crucial negative regulator of the immune system.²⁴ In our study, the TAC C₀/D in patients carrying the CTLA4 rs4553808_AA was significantly lower than that with other genotypes in whole group and CYP3A5 expresser group, but we did not find this association in CYP3A5 nonexpresser. However, Liu et al found that the TAC post-operation concentration at 30 days in kidney transplant recipients carrying G allele of CTLA4 rs4553808 was significantly lower than that in the AA genotype in CYP3A5 expressers, while it was higher in the nonexpressers.²⁵ Our results are inconsistent with Liu’s outcomes, which may be caused by different populations, different analysis methods or different time points for concentration determination.

To exclude the interference of CYP3A5, we divided the patients into CYP3A5 expresser group and CYP3A5 nonexpresser group. In CYP3A5 expresser group, we found the older the age, the higher the TAC C₀/D (see Figure 3). Similarly, Violette Gijsen et al found that in pediatric heart recipients, younger age was associated with a lower TAC concentration.²⁶ This may be due to the changes in the expression and activity of metabolic enzymes with age, resulting in different metabolism of TAC in children, adolescents, and adults.^27–29 As with the results in whole group, MYH9 rs2239781 and CTLA4 rs4553808 were also included in the prediction model of CYP3A5 expresser group. In addition, the TAC C₀/D in patients with MAP3K11 rs7946115_CG was lower than that in patients with homozygous genotype. Mitogen-activated protein kinase kinase kinase 11 (encoded by the MAP3K11) is a member of the serine/threonine kinase family. It is a positive regulator of the JNK signaling pathway. TAC can exert immunosuppressive effects by inhibiting the activation of the JNK pathway,³⁰ so MAP3K11 may indirectly affect the pharmacokinetics of TAC.

In CYP3A5 nonexpresser group, we found that the TAC C₀/D in males was higher than that in females (see Figure 4). Considering that NS mostly occurs in boys, the gender differences in incidence may also be a partial cause of gender variations in TAC C₀/D. However, the underlying mechanism needs to be further explored. Interleukin 2 (IL2) plays an important role in homeostasis and activation of the immune system.³¹ It can only exert its biological effects by binding to the IL2 receptor, which is composed of IL2 receptor alpha chain (encoded by IL2RA gene), beta chain, and gamma chain.³² Zea et al found that murine podocytes express the IL2R, and activation of IL2R will cause podocyte injury, resulting in increased podocyte apoptosis and diminished filtration barrier function.³³ Our data showed that the TAC C₀/D in patients with IL2RA rs12722489_CT was lower than that in patients with homozygous genotype. The mutation of IL2RA rs12722489 may affect the binding of IL2 and IL2R, thereby affecting the function of the filtration barrier, and ultimately the TAC C₀/D. INF2 encodes a member of the formin family, which plays an important role in remodeling the actin and microtubule cytoskeleton.³⁴ Patients with INF2 rs1128880_TT had a higher TAC C₀/D than G allele carriers, whereas INF2 rs1128880_GT showed a lower TAC C₀/D than homozygous genotype. INF2 gene is associated with severe pathological subtypes of NS (focal segmental glomerulosclerosis, FSGS).³⁵ We suspect that the mutation of INF2 rs1128880 may affect the integrity of the filtration membrane, thereafter affecting the excretion of TAC. ACTN4 encoded α-actinin 4, which is an important podocyte cytoskeletal protein and plays a critical role in maintaining podocyte structure and function.³⁶ The mutation of ACTN4 gene is associated with FSGS.³⁷ We found that the TAC C₀/D in patients with CT genotype of ACTN4 rs56113315 and rs3745859 was higher than that in patients with homozygous genotype. Whereas the ACTN4 rs62121818_CC genotype showed a lower TAC C₀/D than T allele carriers, which is consistent with our previous study.¹¹

In addition, we also prospectively recruited 25 children with refractory NS for external validation. The key variables were detected and input into the final models to predict TAC C₀/D (see Table 3 and Figure 4). These predicted C₀/D correlate linearly with the observed C₀/D. This means that these models have a good ability of extrapolation, which can effectively assist doctors and pharmacists in making treatment decisions before administration. This study was limited by the sample size of NS patients and single center. Therefore, we hope to enroll more NS patients to perform multi-center research to test and verify our models. In addition, proteomics and metabolomics can be combined to further explore biomarkers that affect TAC concentration to improve the prediction model’s performance.

Conclusion

In summary, the influence of clinical variables (ALB0, AGE0, and gender) and genetic variables (ACTN4 rs3745859, ACTN4 rs56113315, ACTN4 rs62121818, CTLA4 rs4553808, CYP3A5 rs776746, IL2RA rs12722489, INF2 rs1128880, MAP3K11 rs7946115, MYH9 rs2239781, and MYH9 rs4821478) on TAC C₀/D was described, and the association between genetic polymorphisms of INF2, IL2RA, and MAP3K11 and TAC C₀/D was first reported. In addition, three machine learning-based TAC C₀/D prediction models that integrated clinical and genetic variables were developed, which may support precision dosing of TAC in pediatric patients with refractory NS.