Predicting Blood Concentration of Tacrolimus in Patients With Autoimmune Diseases Using Machine Learning Techniques Based on Real-World Evidence

Tacrolimus is a widely used immunosuppressive drug in patients with autoimmune diseases. It has a narrow therapeutic window, thus requiring therapeutic drug monitoring (TDM) to guide the clinical regimen. This study included 193 cases of tacrolimus TDM data in patients with autoimmune diseases at Southern Medical University Nanfang Hospital from June 7, 2018, to December 31, 2020. The study identified nine important variables for tacrolimus concentration using sequential forward selection, including height, tacrolimus daily dose, other immunosuppressants, low-density lipoprotein cholesterol, mean corpuscular volume, mean corpuscular hemoglobin, white blood cell count, direct bilirubin, and hematocrit. The prediction abilities of 14 models based on regression analysis or machine learning algorithms were compared. Ultimately, a prediction model of tacrolimus concentration was established through eXtreme Gradient Boosting (XGBoost) algorithm with the best predictive ability (R 2 = 0.54, mean absolute error = 0.25, and root mean square error = 0.33). Then, SHapley Additive exPlanations was used to visually interpret the variable’s impacts on tacrolimus concentration. In conclusion, the XGBoost model for predicting blood concentration of tacrolimus on the basis of real-world evidence has good predictive performance, providing guidance for the adjustment of regimen in clinical practice.


INTRODUCTION
Tacrolimus is a calcineurin inhibitor and widely used immunosuppressive drug in solid organ transplant recipients and patients with autoimmune diseases (e.g., lupus nephritis, systemic lupus erythematosus [SLE], and nephrotic syndrome). (Johnston, 2013;Hannah et al., 2016;Mok, 2017a;Hao et al., 2018) The therapeutic window of tacrolimus is narrow, and interindividual variability in dose requirement is high. (Johnston, 2013;Andrews et al., 2017;Mok, 2017b) Therapeutic drug monitoring (TDM) is an important approach for personalized medicine of immunosuppressive drugs, which helps to control drug dosage, minimize drug toxicity, guide therapies, and improve patient care. (Cremers et al., 2016;Zhang and Zhang, 2018) Reasonable plasma tacrolimus concentration can effectively reduce the incidence of adverse drug reactions. According to previous studies, tacrolimus can lead to renal dysfunction, new onset hypertension, and hyperglycemia in patients with lupus nephritis. (Miyasaka et al., 2009;Hannah et al., 2016) In patients with SLE, tacrolimus can result in infections, nephrotoxicity, liver function disorders, nausea, hypertension, anemia, leukopenia, tremors, and itching. (Huang et al., 2019) Furthermore, blood concentrations of tacrolimus appear to be related to acute nephrotoxicity, neurotoxicity, diabetogenicity, and infections in the treatment of lupus nephritis. (Mok, 2017b;Chen et al., 2020) Thus, TDM of tacrolimus is important to reduce adverse effects and ensure adequate drug exposure.
With the rapid development of information technology, the real-world evidence (RWE) derived from medical records has become an important data source for clinical research. RWE is rooted in real clinical practice and comes from a wide range of sources, including electronic medical record examination and imaging data follow-up records during diagnosis and treatment. The research of RWE is a process of data mining, model building, and clinical feature data extraction. Machine learning is suitable for processing a large volume of real-world data, dealing with missing value and high-dimensional data, and capturing complicated relationships between variables, especially for retrospective studies. Recently, some algorithms with more sophisticated principles have been developed, such as eXtreme Gradient Boosting (XGBoost), random forest, K-nearest neighbor (KNN), light gradient boosting machine (LightGBM), and Categorical Boosting (CatBoost), which have been highly recognized in algorithm competitions. (Aha et al., 1991;Breiman, 2001;Chen and Guestrin, 2016;Ke et al., 2017;Prokhorenkova et al., 2017) XGBoost is one of the intelligent classifier construction algorithms, seen in various classification or regression studies with promising prediction results, such as the risk prediction model for type 2 diabetes. (Xu and Wang, 2019) With the increasing number of input subject data, the machine learning model can continually optimize parameters to achieve better accuracy and practicality.
We aimed to identify important influencing variables for tacrolimus concentration on the basis of real-world data and establish a prediction model of tacrolimus concentration in patients with autoimmune disease via machine learning techniques to assist clinical regimen decisions.

Study Population
We enrolled inpatients who were diagnosed with autoimmune diseases and treated with tacrolimus at Southern Medical University Nanfang Hospital from June 7, 2018, to December 31, 2020. The inclusion criteria included the following: 1) patients who were administered tacrolimus 15 days before TDM; 2) patients who were diagnosed with autoimmune diseases, including lupus nephritis, SLE, and nephrotic syndrome. The exclusion criteria included the following: 1) patients with more than 50% absence of major research data; 2) patients who had undertaken transplantation. The workflow of patient inclusion is illustrated in Figure 1. The final dataset included 123 patients with 193 cases of tacrolimus TDM information, and the final dataset was divided into training and testing groups according to the ratio of 8:2.
Study data have been fully deidentified, and confidential information of patients has been deleted, in accordance with the CIOMS/WHO International Ethical Guidelines for Healthrelated Research Involving Humans (2016). Consequently, the study was deemed exempt from informed consent by study participants.

Data Processing
Tacrolimus patient data were extracted from the hospital information system. The interval between two TDM tests of the same patient was mainly distributed within 7-15 days and over 1 month. Hence, samples could be considered independent of each other. We firstly extracted the frequency and dose of tacrolimus ordered by the doctor closest to the TDM test. Drug names (dosage) of tacrolimus included Procofol (0.5 mg*50 tablets), Hongshin (0.5 mg*50 tablets), Serfol (0.5 mg*50 tablets), Neoprocofol (1 mg*50 tablets; sustained-release capsules), and Formexin (1 mg*50 capsules), all in capsule form. Additionally, the last laboratory data before tacrolimus TDM were extracted as essay index, and variables with missing value over 50%, variables completely unrelated to tacrolimus TDM clinically, and classified variables with severe imbalance were deleted. In terms of combination medication, medical order information within 15 days before tacrolimus TDM was extracted, and combination medication was extracted according to drug classification. Given that some drugs did not show in the dataset, the combination medication actually extracted included glucocorticoid (dexamethasone, methylprednisolone, prednisone, cortisone), proton pump inhibitors (PPIs; omeprazole, pantoprazole, ilaprazole, rabeprazole, lansoprazole, ranitidine), calcium channel blockers (CCBs; nifedipine, amlodipine, nitrendipine, felodipine, diltiazem), other immune inhibitors (ciclosporin, mycophenolate mofetil, cyclophosphamide, azathioprine, methotrexate), clarithromycin, and azithromycin. Considering the large amount of missing genetic information, it was not considered in this study.

Variable Selection
On the basis of tacrolimus patient record data, the relevant important variables were screened from multiple influencing factors. To be specific, we included patient's demographic information (e.g., gender, age, height, and weight), medication information (e.g., dosage, frequency, and daily dose of tacrolimus), assay index (e.g., liver function index, kidney function index, and routine blood test), and combination medication (e.g., glucocorticoids, PPIs, CCBs, other immune inhibitors, and clarithromycin/azithromycin). The TDM value of tacrolimus was set as the target variable.
Variable selection was performed using sequential forward selection (SFS) algorithm based on XGBoost to select variable subsets with the minimum size and optimum performance. First, data transformation was carried out on continuous variables using logarithm. Then, all 52 variables were selected using XGBoost model with comparison of model performance, which was measured by R-square (R 2 ), mean absolute error (MAE), and root mean square error (RMSE). R 2 is the squared correlation between predicted and actual TDM values (square root of tacrolimus TDM value), with higher value indicating better predictive ability. MAE is the average of the absolute difference between actual and predicted TDM values, and RMSE is the square root of the mean square error between predicted and actual TDM values; hereon, lower values indicate better predictive ability. We examined the XGBoost model performance of the most prominent variable and identified the point at which there was no considerable gain in R 2 and no considerable loss in MAE and RMSE when adding the next variable to the model. Random forest was used for the interpolation of missing values.

Model Establishment
Using the selected important variables as covariates, we established and analyzed five linear models (linear regression, LASSO regression, ridge regression, elastic net regression, Bayesian ridge regression) and nine machine learning models (KNN, support vector regression [SVR], random forest, XGBoost, LightGBM, CatBoost, NGBoost, AdaBoost, and GradientBoosting). The prediction performance of all models was evaluated through R 2 , MAE, and RMSE. Ultimately, the model with the highest R 2 and lowest MAE and RMSE was selected as the final model to predict tacrolimus TDM value. In the testing cohort, the relationship between the predicted value and the true value of tacrolimus TDM value in the final model was obtained.

Clinical Interpretation
The importance of variables refers to the degree to which each variable in the model contributes to improving the predictive power of the whole model. Herein, we calculated and ranked the importance scores of variables using the algorithm with the best predictive performance. Variables with higher importance scores were more closely related to the accurate prediction of tacrolimus TDM value. Afterward, we used the SHapley Additive exPlanations (SHAP) to visually interpret the impacts of important variables on the model output. (Lundberg and Lee, 2017) Specifically, SHAP could help to explain which variables have positive or negative impacts on predicting tacrolimus TDM value.

Baseline Information
A total of 193 cases of tacrolimus TDM data from 123 patients were included in the study. The study population's baseline information is shown in Table 1. The continuous variables were described by "median  (interquartile range, IQR)," while the classified variables were described by "frequency (percentage, %)." The median tacrolimus TDM value was 4 (IQR 3-6) ng/ml. The median dosage of tacrolimus was 1 mg, the median frequency was 2 times/day, and the median daily dose was 2 mg. The median age of cases in this study was 18 (IQR 8-46) years, the median height was 157 (IQR 122-166) cm, the median weight was 53 (IQR 33-65) kg, the median body mass index was 21 (IQR 17-24) kg/m 2 , and the proportion of male patients was 57.51%. The percentage of using combination medication was 85.49% for glucocorticoids, 48.7% for PPIs, 22.28% for CCBs, 17.1% for other immune inhibitors, and 9.33% for clarithromycin/azithromycin.

Variable Analysis
In Figure 2, the three measurements reached the optimum (R 2 0.42, MAE 0.27, and RMSE 0.37) when nine of the most prominent variables were selected. Although the curves of R 2 , MAE, and RMSE had slight fluctuations after the subset of nine variables, what we pursued was a concise and accurate model with minimal variables but high predictive performance. Hence, the first nine important variables were selected to establish the prediction model, including height, tacrolimus daily dose, other immunosuppressants, low-density lipoprotein cholesterol (LDL-C), mean corpuscular volume (MCV), mean corpuscular hemoglobin (MCH), white blood cell (WBC) count, direct bilirubin (DBIL), and hematocrit (HCT).

Model Performance
The prediction performance of the 14 models is shown in Table 2.
As can be seen, XGBoost had the highest R 2 (0.54) and lowest MAE (0.25) and RMSE (0.33) among all models. We also calculated the accuracy of the predicted concentration within ±30% of the actual concentration of the 14 models to validate the outcome, and XGBoost model gets the highest accuracy, which is 74.4%. Thus, the XGBoost model had optimal ability to predict tacrolimus TDM value. The relationship between the predicted and true values of the tacrolimus TDM value in the XGBoost model is depicted in Figure 3. The red line means predicted value equals to true value. When the blue spots get closer to the red line, the predicted results are more accurate. As illustrated, the majority of blue spots are distributed close to the red line, which can visually demonstrate the predictive accuracy for tacrolimus TDM value. Table 3 shows the ranking of variable importance scores using the XGBoost model, including height, other immunosuppressants, HCT, WBC, tacrolimus daily dose, MCV, MCH, LDL-C, and DBIL, in descending order. According to the score distribution, the impacts of the nine variables on predicting tacrolimus TDM value were balanced, and the importance of height and other immunosuppressants was slightly higher than those of other variables. Furthermore, in Figure 4, SHAP values represent the impacts on model output, which is the prediction of tacrolimus TDM value. Feature value means the contribution of each variable to the predictive power of the model. For variables including tacrolimus daily dose, WBC, MCH, and DBIL, the dot color is redder when the SHAP value becomes larger, while it is bluer when the SHAP value becomes smaller, thus showing the positive impacts of these variables on tacrolimus TDM value. On the contrary, in terms of height, LDL-C, and other immunosuppressants, the dot color is bluer when the SHAP value becomes larger, while it is redder when the SHAP value becomes smaller, thus showing the negative impacts of these variables on tacrolimus TDM value. HCT and MCV did not show significant and regular positive or negative impacts on tacrolimus TDM value.

DISCUSSION
In this study, a machine learning model was established to predict the blood drug concentration of tacrolimus in patients with autoimmune diseases. A feature SFS algorithm was applied to screen model variables, and five linear models and nine machine learning models were used for comparison. Ultimately, the XGBoost algorithm was selected to build the prediction model with an R 2 of 0.54, which achieved a high prediction accuracy in the absence of genetic information. This finding indicated that the XGBoost model has good predictive ability and can be used to predict the blood concentration of patients after taking tacrolimus.
Among the selected important variables in the prediction model, some were identified by previous studies, and some were newfound influencing variables. For instance, height was rarely found to have a relationship with tacrolimus concentration in other studies, whereas it had high importance in our prediction model (importance score 0.1218). Only one population pharmacokinetic model in Asian liver transplant patients mentioned that height influenced the apparent volume of tacrolimus distribution on the basis of whole blood concentrations after oral administration. (Sam et al., 2006) In our study, the sample age distribution was relatively discrete. Among the 193 cases, 100 individuals were aged below 18 years, and their height data were greatly influenced by individual variability. Hence, height was important in the present model. However, future prediction models could be established according to age stratification to avoid bias.  Other immunosuppressants, including ciclosporin, mycophenolate mofetil, cyclophosphamide, azathioprine, and methotrexate, were identified as important variables with negative impacts on tacrolimus concentration. Previous studies mostly discussed the comparison of the therapeutic effects of these drugs. Few investigated the impacts of other immunosuppressants on tacrolimus concentration. As calcineurin inhibitors, the main effect of tacrolimus and ciclosporin was the inhibition of T-cell function. (Schreiber and Crabtree, 1992) Combination drug therapy with ciclosporin may reduce the absorption of tacrolimus in the blood, leading to low concentration. Similarly, other immunosuppressants may influence tacrolimus concentration due to combined therapy. However, data are lacking to explain this relationship, and more research is needed in the future.
The most commonly discussed variable influencing tacrolimus concentration was HCT. Several studies on transplant recipients have identified HCT to affect tacrolimus clearance. (Undre and Schäfer, 1998;Lee et al., 2006;Brunet et al., 2019) To be more specific, Staatz C, et al. and Sam W, et al. have quantified the impact of HCT on tacrolimus apparent clearance; they found that an increase in tacrolimus clearance was associated with decreased HCT levels. (Sam et al., 2000;Staatz et al., 2002) Delays in drug clearance indicate the accumulation of tacrolimus in the body, leading to higher concentration. In addition, low level of HCT may reduce the proportion of tacrolimus bound to red blood cells and increase the proportion bound to plasma, which can be more easily metabolized by the liver. (Li et al., 2011;Brunet et al., 2019) Therefore, low HCT level could result in a decrease in total tacrolimus concentration in the whole blood. This could also interpret the relationship between MCV and tacrolimus blood concentration. Nevertheless, the specific quantified impacts of HCT and MCV on tacrolimus blood concentration in our study were unclear. Larger samples may be needed for investigation in the future. Regarding MCH, low HCT level may lead to low MCH value. Moreover, it was positively associated with tacrolimus concentration in the present study.
A study on Chinese renal transplant recipients taking tacrolimus as immunosuppressive drug found that total bilirubin was a significant influencing factor on maintenance tacrolimus dose. (Li et al., 2011) Our results indicate that DBIL was significantly related to tacrolimus concentration and showed a reverse association. According to the findings of Bellarosa, et al., bilirubin and tacrolimus competitively bind with P-glycoprotein; thus, high bilirubin level may delay the transport and metabolism of tacrolimus, resulting in high blood concentration. (Bellarosa et al., 2009) Additionally, high bilirubin level indicates impaired liver function and poor liver metabolism, which may lead to increased tacrolimus concentration in the blood.
In terms of other influencing variables, daily dose shows obvious positive relationship with tacrolimus concentration, illustrating that using high dose of tacrolimus daily could contribute to high tacrolimus blood concentration. WBC was positively related to tacrolimus concentration, and LDL-C was negatively related to tacrolimus concentration, both of which were rarely investigated in previous studies. Increased WBC counts imply that the body is attacked by an infection, requiring intensive drug treatment. Thus, the use of tacrolimus is probably increased, leading to higher blood concentrations. High level of LDL-C could represent liver dysfunction and was theoretically associated with high blood concentration of tacrolimus. However, our study shows the opposite findings. These results should be validated in the future with larger sample size and higher quality data. Previous studies found that some factors, including CYP3A5 genotype, bodyweight, liver and renal function, serum albumin, and coadministration of diltiazem and fluconazole, are related to the concentration, dose, and clearance of tacrolimus; however, they were not significant in our model. (Zahir et al., 2005;Sam et al., 2006;Li et al., 2011;Brunet et al., 2019) In this study, genetic data were too scarce to be used for exploration, and further research should include more factors in the future.
Most previous research established pharmacokinetic models. In real-world studies, variables are not always independent to each other, and the majority have closely non-linear relationships. According to the second consensus report about TDM of tacrolimus-personalized therapy issued by the Immunosuppressive Drugs Scientific Committee of the International Association of Therapeutic Drug Monitoring and Clinical Toxicity in 2017, R 2 values have great fluctuation with a range from 0.27 to 0.99 in the prediction of tacrolimus concentration using pharmacokinetic models. (Brunet et al., 2019) Our prediction model using XGBoost algorithm achieved an R 2 of 0.54 and the lowest MAE and RMSE, which outperformed other models investigated herein. One study used the SVR model to predict tacrolimus blood concentration in liver transplantation patients; the authors found that SVR models show superior prediction accuracy than multiple linear regression model. (Van Looy et al., 2007) Compared with pharmacokinetic models and linear regression models, machine learning models can capture the complex relationships of variables and analyze high-dimension data in clinical practice. Considering that our data were derived from real-world settings with some missing values, machine learning models have better ability to deal with incomplete data.
In conclusion, on the basis of 193 cases of tacrolimus TDM data from real-world settings, the study identified nine important variables for tacrolimus concentration, including height, other immunosuppressants, HCT, WBC, tacrolimus daily dose, MCV, MCH, LDL-C, and DBIL, and their specific quantified impacts. Then, we established a prediction model of tacrolimus concentration in patients with autoimmune disease via XGBoost algorithm, which had better prediction abilities than other models compared in the study.
To our knowledge, this study is the first to predict the blood concentration of tacrolimus in patients with autoimmune diseases using machine learning techniques on the basis of RWE. The prediction model is of important clinical significance as regimen adjustment may be required to maintain blood concentrations within the therapeutic range in patients whose clinical data are changing. One limitation is the limited sample size, and all data were obtained from one center. In the future, larger data from multiple centers are needed for further analysis.