In Silico Prediction and Insights Into the Structural Basis of Drug Induced Nephrotoxicity

Drug induced nephrotoxicity is a major clinical challenge, and it is always associated with higher costs for the pharmaceutical industry and due to detection during the late stages of drug development. It is desirable for improving the health outcomes for patients to distinguish nephrotoxic structures at an early stage of drug development. In this study, we focused on in silico prediction and insights into the structural basis of drug induced nephrotoxicity, based on reliable data on human nephrotoxicity. We collected 565 diverse chemical structures, including 287 nephrotoxic drugs on humans in the real world, and 278 non-nephrotoxic approved drugs. Several different machine learning and deep learning algorithms were employed for in silico model building. Then, a consensus model was developed based on three best individual models (RFR_QNPR, XGBOOST_QNPR, and CNF). The consensus model performed much better than individual models on internal validation and it achieved prediction accuracy of 86.24% external validation. The results of analysis of molecular properties differences between nephrotoxic and non-nephrotoxic structures indicated that several key molecular properties differ significantly, including molecular weight (MW), molecular polar surface area (MPSA), AlogP, number of hydrogen bond acceptors (nHBA), molecular solubility (LogS), the number of rotatable bonds (nRotB), and the number of aromatic rings (nAR). These molecular properties may be able to play an important part in the identification of nephrotoxic chemicals. Finally, 87 structural alerts for chemical nephrotoxicity were mined with f-score and positive rate analysis of substructures from Klekota-Roth fingerprint (KRFP). These structural alerts can well identify nephrotoxic drug structures in the data set. The in silico models and the structural alerts could be freely accessed via https://ochem.eu/article/140251 and http://www.sapredictor.cn, respectively. We hope the results should provide useful tools for early nephrotoxicity estimation in drug development.


INTRODUCTION
Drug induced nephrotoxicity (DIN) can be defined as any renal injury caused directly or indirectly by medications (Sales and Foresto, 2020), which has been a major issue for patients and the pharmaceutical industry. The real world data (RWD) showed that incidence of drug induced nephrotoxicity to be approximately 14-26% in adult populations (Mehta et al., 2004;Uchino et al., 2005;Hoste et al., 2015).
The syndrome always contributes to considerable morbidity, mortality, and high costs, and it can lead to the development of chronic kidney disease (CKD) or end-stage renal disease (ESRD). Nephrotoxicity has become an important concern in modern drug development. It is one of the major reasons for safety-related failures at all phases of drug development and even marketed drugs being restricted or withdrawn .
The mechanisms of drug induced nephrotoxicity were very complex and may be different between various drug classes. Based on the histological component of the affected kidney, the drug induced nephrotoxicity should be categorized as three mechanisms, including proximal tubular injury and acute tubular necrosis (ATN), tubular obstruction by crystals or casts containing drugs and their metabolites, and interstitial nephritis induced by drugs and their metabolites (Nolin and Himmelfarb, 2010;Kwiatkowska et al., 2021). On the other hand, the mechanisms of drug induced nephrotoxicity can also be classified according to the mode of action of the drugs, including cytotoxicity (necrosis or apoptosis), immune injury, and ischemic injury. With the help of transdermal transport system, especially organic anion transporter 1 (OAT1), drugs were accumulated in proximal convoluted tubule epithelial cells. The cell necrosis or apoptosis would be caused when high concentration was reached, then cell necrosis or apoptosis would be caused (Sekine and Endou, 2009). It was reported that drugs can act on mitochondria and block production of adenosine triphosphate (ATP), resulting in cell necrosis or apoptosis (Gai et al., 2020). Nephrotoxic drugs can also increase superoxide free radical production and decrease antioxidant free radical production in epithelial cells, which can also lead to cell necrosis or apoptosis (Paller et al., 1984). Ferroptosis was a type of cell death characterized by iron overload and accumulation of toxic lipid peroxides . In recent years, studies have shown that ferroptosis was closely related to drug induced nephrotoxicity, especially acute kidney injury, but the exact molecular biological mechanism has not been clarified, and this needs more research (Hu et al., 2019). The immune response caused by some drugs acting as antigens or haptens can result in inflammation of the blood vessels and tubules of the kidney, such as penicillin-induced interstitial nephritis (Spanou et al., 2006). This kind of immune damage was related to individual hypersensitivity to drugs, so there were significant individual differences. Besides, nephrotoxic drugs can cause renal ischemia and injury by reducing blood perfusion in renal tissues (Raza and Naureen, 2020). For example, nonsteroidal anti-inflammatory drugs (NSAIDs) were able to reduce renal blood flow by changing the resistance of glomerular entry and exit arteries, and then induce renal ischemic injury (Hörl, 2010). Unfortunately, there are few studies on the structural characteristics of nephrotoxic drugs.
The evaluation of the nephrotoxic potential of chemicals in early stage is quite important and useful for reducing the failure of drug development. The in vivo test for drug induced nephrotoxicity evaluation is always very complicated, costly, and time-consuming, and it is not suitable for screening a large number of chemicals, and especially for virtual structures. Besides, the experimental results are easily affected by various factors such as model animals and technology and environment. Compared with biological experimental methods, the use of computational toxicology for nephrotoxicity estimation of compounds has obvious advantages: (1) large quantities of compounds can be rapidly processed and predicted; (2) the toxicity of the compound can be predicted by computational models as long as the structure is known, even if the compound has not been synthesized; and (3) computational toxicological methods can also contribute to the study of the mechanisms. Consequently, it should make a lot of sense to develop fast and accurate computational tools to estimate the risk of nephrotoxicity. Over the past decades, many computational models have been developed for toxicity prediction, but only a few models reported related to kidney injury or urinary tract toxicity, and due to the variety and complexity of the symptoms and mechanisms. Lei et al. summarized the reported models related to urinary tract toxicity until 2017 (Lei et al., 2017). Most of them were established based on biomarker descriptors, and only five models were based on theoretical descriptors with the datasets of drugs in various development phases. In Lei et al.'s study, they developed a series of qualitative and quantitative structure activity relationship (QSAR) models for urinary tract toxicity prediction using 258 compounds, the best regression model reached q 2 ext of 0.845 for the test set, and the best classification model gave global accuracy of 90.77% for the test set. Zhang et al. developed an in silico prediction model for chemical induced urinary tract toxicity using naïve Bayes classifier based on mouse intraperitoneal data set . The model provided 84.2% overall accuracy for the external test set. They also obtained several important molecular descriptors and fragments. More recently, Sun et al. developed QSAR models for screening nephrotoxicity of the ingredients in TCMs based on natural product or mixed dataset (Sun et al., 2019). The models performed well on external validation with 30 ingredients in the TCMs. The published models related to nephrotoxicity always provided high statistical performance. However, the structural characteristics of nephrotoxic and non-nephrotoxic drugs were rarely analyzed in these studies, and the usefulness of most published models was restricted because of poor availability. Besides, it should be more useful to develop the models based on real world data with human nephrotoxicity, and due to the species specificity in drug toxicity between rodent animals and human beings.
In the present study, we focused on the in silico prediction and insights into the structural basis of drug induced nephrotoxicity based on the medications causing human nephrotoxicity in the real world.

Data Source and Preparation
In this study, only medications with human nephrotoxicity data were included. The nephrotoxic structures were extracted from the Side Effect Resource (SIDER) database (Kuhn et al., 2016). SIDER is a widely used database of adverse drug reactions (ADRs), which contained the information on approved drugs and the ADRs on humans. Herein, we retrieved the entire SIDER database and collected those drugs with nephrotoxicity related ADRs with frequency ≥0.1% in the real world. The corresponding structures of included chemicals were downloaded in smiles format from the PubChem database (Kim et al., 2016). The non-nephrotoxic structures were extracted from Zhang's work . They built the negative drug dataset with the drugs without nephrotoxicity from the SIDER database. The included structures were carefully prepared as follows: (1) removing duplicate substances; (2) keep only the main ingredients in mixtures; and (3) salts were converted to their parent forms.

Principal Component Analysis for the Definition of the Chemical Space
The sufficiently structural diversity was a key issue for global QSAR models to ensure a reasonable predictive accuracy (Ancuceanu et al., 2019). Principal component analysis (PCA) is a well-known technique for reducing the dimensionality and increasing interpretability, which can solve an eigenvector problem by creating new uncorrelated variables that successively maximize variance (Jolliffe and Cadima, 2016). In this study, the chemical space of data sets was analyzed with the first two principal components of CDK (Chemistry Development Kit) Descriptors. The PCA was performed using SPSS Statistics 26.

Algorithms For Model Building
The model building was performed on the online chemical database and modeling environment (OCHEM), which is a user friendly web-based platform for automatic and simple QSAR modeling (Sushko et al., 2011). OCHEM supports the typical steps of QSAR modeling, and the models can be published and publicly used on the web (Oprisiu et al., 2013;Cui et al., 2019;Pawar et al., 2019;Cui et al., 2021;Huang et al., 2021;Ta et al., 2021). Among the many state-of-the-art modeling methods available on OCHEM, we applied five widely used traditional machine learning (ML) approaches and five different deep learning (DL) algorithms. As an application of artificial intelligence (AI), ML has been an effective tool for modeling in computational toxicology. Five highly effective and robust ML approaches were used in this study, including associative neural network (ASNN) (Tetko, 2009), support vector machine (SVM) (Chang and Lin, 2011), C4.5 decision tree (WEKA J48) (Hall et al., 2009), random forest (RFR) (Breiman, 2001), and extreme gradient boosting (XGBoost) (Chen and Guestrin, 2016). DL is an extension of machine learning, and its concept came from the research of artificial neural networks. In this study, we used five different DL approaches, including convolutional neural network fingerprint (CNF) , transformer convolutional neural network (TRANSNN) (Karpov et al., 2020), TEXTCNN algorithm available from DeepChem (TEXTCNN) (Wu et al., 2017), Graph Isomorphism Network (GIN) (Capela et al., 2019), and edge attention based multirelational graph convolutional networks (EAGCNG) (Shang et al., 2018). The detailed descriptions of these ML and DL approaches can be found in the corresponding literature.
The individual parameters for each model algorithm were optimized automatically by the method itself in an inner loop of cross-validation. For instance, the SVM algorithm used libSVM (Chang and Lin, 2011) on OCHEM. There were two important configurable options for this method, including SVM type (ε-SVR and μ-SVR, etc.) and the kernel type (linear, polynomial, radial basis function, sigmoid, etc.). In the OCHEM workflow, classic ε-SVR and radial basis function kernels were used. The other SVM parameters, namely cost C and width of the RBF kernel (γ, g), were optimized using default grid search, and this was performed according to the LibSVM manual .
For the deep learning models, the SMILES string of each compound was served as the input without descriptors.

Consensus Modeling
Consensus modeling is to unify the prediction of unknown samples by multiple individual models to achieve a unified result, and then improve the prediction of the model. By averaging the prediction of individual models, noise can be reduced, thus, and the consensus model can provide better predictive power than most individual models alone (Tapia Garcı'a et al., 2012). Herein, the consensus model was built with simple average of predictions from the best performed individual models.

Applicability Domain Assessment
For QSAR models, it is important to estimate the applicability domain (AD) to determine whether the test compound is suitable for the specific model. In this study, we determined the AD with distance to model (DM) in OCHEM proposed by Sushko (Sushko et al., 2011). The DM assesses the distance from the target compound to the model. The larger DM means the lower applicability for the target compound.

Analysis of Molecular Properties Differences Between Nephrotoxic and Non-nephrotoxic Structures
The molecular properties of compounds always can make a significant difference in the toxicity. In this study, the analysis of differences of molecular properties between nephrotoxic and non-nephrotoxic structures was performed, in order to investigate the relevance of these molecular properties with drug induced nephrotoxicity. Several commonly used physicalchemical properties which have been widely adopted in the analysis for other endpoints were calculated and analyzed, including molecular weight (MW), molecular polar surface area (MPSA), AlogP, molecular solubility (LogS), the number of hydrogen bond acceptors (nHBA) and donors (nHBD), the number of rotatable bonds (nRotB), and the number of aromatic rings (nAR). These properties are relative to the molecular size, lipophilicity and solubility, and hydrogen bonding ability and complexity, respectively. The comparison between groups was tested by T-test, and the p value <0.05 was considered to indicate statistical significance. The molecular properties were calculated with PaDEL-Descriptor package (Yap, 2011).

Identification of Structural Alerts Responsible for Nephrotoxicity
Structural alert (SA), or privileged substructure, was defined as the substructure which can cause the chemicals to become toxic. It has been well accepted in toxicity assessment, because of the direct derivation from mechanistic knowledge. SAs have been commonly used for toxicity assessment of many different endpoints (Claesson and Minidis, 2018;Yang et al., 2018;Wedlake et al., 2020). In this study, we identified the structural alerts for nephrotoxicity by calculating f-score and positive rate of each fragment from Klekota-Roth fingerprint (KRFP, 4,860 bits). The specific substructure should be regarded as a SA if presented more frequently in nephrotoxic drugs than non-nephrotoxic drugs. The positive rate (PR) of a substructure was defined as below: where N fragment_positive was the number of nephrotoxic drugs containing the fragment, and N fragment was the total number of drugs containing the fragment.

Validation and Evaluation of Models
All the ML and DL models were first internally validated with 5-fold cross validation, and the best performed models were further validated with the external validation set. The structural alerts were assessed with the whole dataset. The classifiers and SAs were evaluated based on the counts of true positives (TP), false positives (FP), true negatives (TN), and false negatives (FN). Several statistical parameters were calculated to evaluate the classifiers, including the total accuracy (Q), sensitivity (SE), specificity (SP), enrichment factor (EF), and Matthews correlation coefficient (MCC), which were calculated with Eqs (2-6).
Additionally, the receiver operating characteristic (ROC) curve was also plotted for the QSAR models, and the values of area under the ROC curve (AUC) were provided, too.

Data Set Analysis
In the study, 565 diverse structures were kept after preparation, including 287 nephrotoxic drugs and 278 non-nephrotoxic drugs. As shown in Table 1, the whole data set was randomly divided into a training set with 456 chemicals (232 nephrotoxic and 224 non-nephrotoxic) and an external validation set with 109 chemicals (55 nephrotoxic and 54 non-nephrotoxic). The structures of the drugs can be found in Supplementary Table S1.
The diversity of structures is a crucial factor for the applicability of global models. Thus, we performed the principal component analysis (PCA) based on CDK descriptors to analyze the chemical space of the included compounds. PCA can simplify the complexity in high-dimensional data while retaining trends and patterns by transforming the data into fewer dimensions, which act as summaries of features (Ringnér, 2008). The first two principal components were kept and used for the definition of the chemical space of compound in the training and validation sets in this study. As shown in Figure 1, the distribution scatter diagram illustrated that the data sets shared a similar chemical space. In addition, the Tanimoto similarity index (TSI) was also calculated based on the ECFP-4 fingerprint to evaluate similarities among the structures. The average value of the entire data set was 0.13, which indicated an evidently chemical diversity of the entire data set.

Results of ML and DL Models
Combined with five different ML algorithms and eight molecule descriptor packages, 40 ML models were developed. Meanwhile, five DL models were developed with different algorithms using chemical SMILES-string as input. The performances of ML and DL models on 5-fold cross validation were shown in Table 2. The prediction accuracy (Q) of the models ranged from 58.54 to 73.90%. Among them, three models performed much better on 5-fold cross validation than others, including a DL model (CNF) and two ML models (XGBoost_QNPR and RFR_QNPR). These models all provide good predictive ability with Q value ≥ 70.00% and AUC value ≥0.80.
The consensus approach has been proved able to improve the accuracy of models. In this study, a consensus model was developed based on the three best models. The consensus model performed much better than the individual models. As shown in Table 2

External Validation of Models
Due to the complete independence from the model training, the external validation set could be well used for evaluating the predictive ability of models objectively. As shown in Table 3; Figure 2, the best performed individual models also achieved good predictive results on external validation. The results of ML model building and validation suggested that ML models developed with the QNPR descriptor performed much better than others. The QNPR descriptors were derived directly from the chemical SMILES strings. For each structure either canonical SMILES or IUPAC name would be split into fragments of a specified length determined by the configuration (Thormann et al., 2007). For the different ML algorithms, XGBoost and RFR performed better than others. XGBoost is a scalable end-to-end tree boosting system, which has been widely used in many application scenarios. It supports multi-thread computation and uses regularization enhancement technology to reduce over-fitting, so as to ensure the robustness of the model. Meanwhile, it has the advantages of flexibility, fast calculation speed, and good robustness. Therefore, XGBoost is not easy to be disturbed by outliers. It can achieve state-of-the-art results on  many machine learning challenges (Chen and Guestrin, 2016). RFR is an ensemble learning method by constructing a multitude of decision trees (DT) at training time (Breiman, 2001). It can help to improve the accuracy by reducing overfitting in decision trees and automating missing values present in the data. CNF contributed to the best DL model. CNF is one of the state-of-theart deep learning methods, which is based on ideas of text processing. It could achieve the high prediction accuracy due to the augmentation technique employed during both training and inference steps. The DL model in our study did not show significant better predictive power compared with the ML model. However, considering that the DL models did not require additional molecular description, as long as SMILES are provided, the DL algorithms still showed obvious advantages over ML algorithms.

The Differences of Molecular Properties Between Nephrotoxic and Non-nephrotoxic Drugs
The molecular physical-chemical properties could provide useful information for biological activities. Herein, we analyzed the differences of several commonly used physical-chemical properties between nephrotoxic and non-nephrotoxic drugs. The distributions of these descriptors for nephrotoxic and non-nephrotoxic drugs can be seen in Figure 3.
The characteristics MW and MPSA can simply assess the size and complexity of compounds. For the entire data set in the study, the values of MW were distributed between 59.04 and 4,491.88, with a mean of 416.86. The mean value was 459.71 for nephrotoxic drugs and 372.63 for non-nephrotoxic drugs. The difference between the mean MW of nephrotoxic and nonnephrotoxic drugs was not significantly different (p < 0.01). The values of MPSA were distributed from 0 to 1902.88, with a mean of 119.80 for the entire dataset. The mean value was 137.52 for nephrotoxic drugs and 101.51 for non-nephrotoxic drugs (p < 0.01). These results suggested that it may be significantly different in structure size and polar surface area between nephrotoxic and non-nephrotoxic compounds.
The AlogP is commonly used to represent the lipophilicity of compounds. For the entire data set, the values of AlogP ranged from −32.81 to 10.25, with a mean of 1.48. The mean value of AlogP was 1.34 for nephrotoxic chemicals and 1.64 for non-nephrotoxic drugs. As can be seen in Figure 3, the distributions were not significantly different between nephrotoxic and non-nephrotoxic drugs with p-value 0.34. The result indicated that chemical lipophilicity may be weakly correlated with drug induced nephrotoxicity.
Chemical hydrogen bonding ability also is an important character for its activity and toxicity, and it was usually in terms of nHBA and nHBD. In this data set, mean values of nHBA for nephrotoxic and non-nephrotoxic drugs were 7.06 and 5.41, respectively, and the difference was significant (p < 0.01). Meanwhile, the nephrotoxic and non-nephrotoxic drugs had mean values of nHBD with 3.29 and 2.47, respectively, and the difference is not significant (p 0.07). The data indicated that nHBA was obviously associated with drug induced nephrotoxicity while nHBD was not.
LogS is an estimation of molecular solubility in water. For the entire data set, the values of LogS were distributed from −35.86 to 2.45, and the mean value was −4.52. For nephrotoxic structures, the mean value of LogS was −4.81, and it is −4.21 for non-nephrotoxic drugs. The difference also proved statistical significance (p 0.03). This result demonstrated that there was difference of molecular solubility between nephrotoxic and non-nephrotoxic drugs.
As shown in Figure 3, drug induced nephrotoxicity was also obviously associated with nRotB (the mean values were 8.05 for nephrotoxic chemicals, and 5.78 for non-nephrotoxic drugs, and with p-value 0.01) and nAR (the mean values were 1.70 for nephrotoxic chemicals and 1.08 for non-nephrotoxic drugs, with significant difference (p < 0.01).
Through the analysis of molecular properties, several physicalchemical properties have obvious differentiating effect on drug induced nephrotoxicity. Nevertheless, nephrotoxicity is a complex endpoint yet. It's not easy to explain the mechanism of drug induced nephrotoxicity with individual simple chemical descriptors.

Structural Alerts Responsible for Nephrotoxicity
The structural alerts (SA) responsible for nephrotoxicity were identified using f-score and positive rate analysis of each fragment from KRFP fingerprint. Only the fragments existing in 6 or more drugs were kept. The fragments with f-score ≥0.005 and positive rate ≥0.75 were identified. Finally, 87 representative fragments were filtered and listed in Supplementary Table S2. Among them, 16 substructures presented in nephrotoxic active chemicals only, which    Table 4.
We analyzed the structural characteristics of drugs in the entire data set, and 75.04% structures (424/565) were correctly classified. Among the 290 drug structures contained at least one identified substructure, 218 structures were true nephrotoxic. The classification accuracy was 75.17%. For the drug structures that did not contain any identified substructure, 206 of 275 drugs (74.91%) were true non-nephrotoxic. The frequencies of these substructures were much higher in nephrotoxic drugs than non-nephrotoxic drugs, and showed a good ability to distinguish nephrotoxic drugs in the whole data set. To a certain degree, these fragments could be considered as the structural alerts responsible for nephrotoxicity. If a structure contains one or more structural alerts, it is more likely to be nephrotoxic than non-nephrotoxic.
Fluorine is widely used in medicinal chemistry to improve a molecule's potency and permeability. However, the metabolism of fluorinated compounds may produce fluoride and other toxic metabolites (Pan, 2019). As shown in Table 4, several selected structural alerts contained phenyl fluoride (Nos 3,8,9,11,12,and 14). Since the kidney is a main target organ of mammalian fluoride systemic exposure and renal toxicity can occur after acute and chronic fluoride intoxication (Quadri et al., 2016). The kidney plays a vital role in fluoride metabolism, since 50-80% of fluoride is removed via urinary excretion. Fluorinated compounds can be toxic to the kidneys in a number of ways. It has widely been reported that fluorinated compounds can increase the generation of reactive oxygen species (ROS) and free radicals, cause extensive oxidative stress and excessive lipid peroxidation, and reduce antioxidant enzyme activities in vivo or in vitro (Quadri et al., 2016). Thus, numerous renal structural, ultrastructural, and functional may be changed after receiving increased amounts of fluorinated compounds exposure. Certainly, this is not meant to sound the alarm on all fluorinated compounds, we suggest raising the awareness of common drug instability and metabolism issues leading to defluorination, as well as the resulting reactive/toxic metabolite (Pan, 2019). Polyamines and their derivatives (Nos 2 and 15) widely existed in aminoglycosides and the other nephrotoxic drugs. Polyamines catabolism can be stimulated through oxidation, which will lead to the generation of ROS and low abundance of free polyamines with antioxidant capacity (Pegg, 2013;Murray Stewart et al., 2018). Consequently, this process is associated with a number of pathologies, including cancer, neurological disorders, and kidney dysfunction. Recent studies have also suggested a role for polyamines derivatives in the p53mediated ferroptotic response to ROS stress. It is an irondependent and nonapoptotic mode of cell death and has been associated with drug-induced nephrotoxicity recently (Ou et al., 2016;Murray Stewart et al., 2018). Benzimidazoles derivatives were always proposed as new bioreductive prodrugs with the potential anticancer activity, due to their effect on the DNA destruction and growth inhibition into selected tumor cell lines (Błaszczak-Świątkiewicz et al., 2014). Besides, benzimidazoles derivatives also can induce concurrent apoptotic and pyroptotic cell death (Ren et al., 2021). However, when these processes take place in kidney cells, kidney damage can occur. Toluene-contained structures (Nos 1 and 7) were also selected as structural alerts for drug induced rhabdomyolysis (Cui et al., 2019), which always associated to oliguric renal failure. In the setting of toluene intoxication, electrolyte disturbances may play important roles on causing rhabdomyolysis (Camara-Lemarroy et al., 2015). Toluene has also been associated with direct induction of acute tubular necrosis and acute oliguric renal failure. Among the SAs we identified, there were several fragments that also could cause the formation of ROS (Nos 6 and 16, etc.), which may contribute to the nephrotoxicity.
In the present study, the SAs were identified with the f-score and frequency analysis of defined substructures. These methods have been widely used for the SA discovery for many other endpoints. In fact, these methods have some shortcomings. They are not able to characterize the spatial arrangement of identified substructures, and they cannot make a good distinction when more than one SA presented in the same structure. In spite of this, these structural alerts were able to well distinguish chemical structures with renal toxicity, and they can help to understand the specific fragments which lead to nephrotoxicity. Therefore, these SAs should be severed as a useful tool to visually evaluate the nephrotoxicity of chemicals.

Availability of QSAR Models and Structural Alerts
For ease of use, the QSAR models were made available at OCHEM. The consensus model could be accessed via https:// ochem.eu/article/140251. The three best individual models (RFR_QNPR, XGBOOST_QNPR, and CNF) were also available with the corresponding model IDs. Users can predict chemical nephrotoxicity by using the "Apply the model to new compounds" link. In addition, the data sets for modeling could be downloaded by using the "Export this basket" link.
The structural alerts responsible for nephrotoxicity have been integrated as part of our Web server SApredictor, which is a structural alert based expert system for drug toxicity prediction and freely available at http://www.sapredictor.cn. With the help of SApredictor, people can quickly evaluate whether the query chemicals are nephrotoxic, and the specific structural fragments that lead to the nephrotoxicity of the compounds will be intuitively shown to provide valuable reference for the modification of the structures.

CONCLUSION
In this study, we collected 287 drug structures which were proved nephrotoxic on humans in the real world. A comparable amount of non-nephrotoxic structures was also extracted from approved drugs. Then, in silcio models were developed using OCHEM tools. A total of 40 ML models were developed using 5 different machine learning algorithms along with 8 descriptor packages. Besides, 5 DL models were also developed using different deep learning methods. Among them, two ML models (RFR_QNPR and XGBoost_QNPR), and one DL models (CNF) provided best predictive ability. A consensus model was developed based on them, which performed much better on internal validation, and provided good predictive ability on external validation. The consensus model and the best individual models were freely available at https://ochem.eu/article/140251. Moreover, the differences of several commonly used physical-chemical properties between nephrotoxic and non-nephrotoxic drugs were investigated. The results indicated that several key molecular properties differ significantly between nephrotoxic and nonnephrotoxic structures, including molecular weight (MW), molecular polar surface area (MPSA), AlogP, number of hydrogen bond acceptors (nHBA), molecular solubility (LogS), the number of rotatable bonds (nRotB), and the number of aromatic rings (nAR). Thus, these molecular descriptors may be associated to drug-induced nephrotoxicity and could play an important role in the identification of nephrotoxic chemicals. Finally, we identified the structural alerts responsible for nephrotoxicity using f-score and positive rate analysis. There were 87 structural alerts identified from the fragments of KRFP fingerprint. A compound would be classified as nephrotoxic if it contains one or more such SAs. These structural alerts showed a good ability to distinguish nephrotoxic drugs in the entire data set. They have been integrated as part of our web server SApredictor, which is freely available at www.sapredictor.cn. The in silico models and the structural alerts could be useful tools for estimation of nephrotoxicity in drug discovery.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials, and further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
XL contributed to conception and design of the study. YS and YH collected the datasets and carried out the experiments. YS, BW, and RZ performed the analysis. YS, YH, and XL Frontiers in Pharmacology | www.frontiersin.org January 2022 | Volume 12 | Article 793332 interpreted the results and wrote the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This work was supported by the National Natural Science Foundation of China (Grant 81803433) and the Special Research project of Clinical Toxicology of Chinese Society of Toxicology (CST2020CT104).