In Silico Prediction of Hemolytic Toxicity on the ... - ACS Publications

Article pubs.acs.org/jmc

Cite This: J. Med. Chem. XXXX, XXX, XXX−XXX

In Silico Prediction of Hemolytic Toxicity on the Human Erythrocytes for Small Molecules by Machine-Learning and Genetic Algorithm Suqing Zheng,*,†,‡ Yibing Wang,§,∥ Wenxin Liu,† Wenping Chang,† Guang Liang,†,‡ Yong Xu,⊥ and Fu Lin*,† †

School of Pharmaceutical Sciences, Wenzhou Medical University, Wenzhou 325035, Zhejiang P. R. China Chemical Biology Research Center, Wenzhou Medical University, Wenzhou 325035, Zhejiang, P. R. China § Genetic Screening Center, National Institute of Biological Sciences, Beijing 102206, P. R. China ∥ Tsinghua Institute of Multidisciplinary Biomedical Research, Tsinghua University, Beijing 100084, P. R. China ⊥ Center of Chemical Biology, Guangzhou Institute of Biomedicine and Health, Chinese Academy of Sciences, Guangzhou 510530, Guangdong P. R. China

Downloaded via KEAN UNIV on July 18, 2019 at 06:41:38 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

S Supporting Information *

ABSTRACT: Hemolytic toxicity of small molecules, as one of the important ADMET end points, can cause the lysis of erythrocytes membrane and leaking of hemoglobin into the blood plasma, which leads to various side effects. Thus, it is very crucial to assess the hemolytic potential of small molecules during the early stage of drug development process. However, so far there is no computational model to predict the human hemolytic toxicity of small molecules. To this end, we manually curate the hemolytic toxicity data set for the small molecules experimentally evaluated on the human erythrocytes, develop the first machine-learning (ML) based models to predict the human hemolytic toxicity of small molecules, harness the genetic algorithm (GA) and ML based model to optimize human hemolytic toxicity based on the molecular fingerprint to derive “optimal virtual fingerprints (OVFs)” with the desired hemolytic/ nonhemolytic property, and finally implement a free software for the users to predict/ optimize the human hemolytic toxicity with ML and GA in the automatic manner.

■

saponin molecule named “dioscin” from the molecular level. By use of MD simulation, we can directly observe that how dioscin molecules penetrate into the membrane, interact with the cholesterols in the membrane, accumulate in the lipid raft, and induce the dramatic curvature change of lipid raft.14 In spite of this intuitive merits of MD simulation technique, it is still computationally demanding to perform the virtual screening against the huge compound database to identify the possible hemolytic compounds by the large-scale MD simulations with the membrane. Therefore, the fast machinelearning method could be a better option to predict the hemolytic toxicity of small molecules. Hitherto there are two works about the hemolytic peptide prediction based on the peptide sequence15,16 and one recent work about the prediction of hemolytic saponins, which are a special type of compounds consisting of a hydrophobic steroid/triterpenoid moiety and hydrophilic carbohydrate branches. However, all the prediction models are derived from the hemolytic/ nonhemolytic peptides (linear peptides with the natural amino acids) or saponins, which are experimentally assayed

INTRODUCTION ADMET prediction plays more and more critical role in the identification and optimization of lead compounds during the drug discovery process.1−9 Hemolytic toxicity, as one of the crucial ADMET end points, can cause the disruption of erythrocyte membrane and subsequently trigger the release of hemoglobin into the blood plasma,10 which consequently can incur various side effects such as acute renal failure.11,12 Therefore, it is imperative to evaluate the hemolytic potential of lead compounds during the research and development (R&D) of novel drugs, especially for the parenteral drugs. However, experimental assessment of hemolytic toxicity on the human erythrocytes is a slow, expensive and tedious process, because in vitro hemolysis assay for each compound requires the onerous exploration of cosolvent by trial-and-error to find an optimal one with the sufficient solubilizing capacity, and minimum hemolytic toxicity introduced by cosolvent.13 Obviously, this time-consuming experimental hemolysis assay will dramatically restrict the rapid high-throughput screening of chemical database for the candidate hemolytic/nonhemolytic compounds of interest. Hence the computational prediction could be a good alternative prior to the laborious experiments. In our previous work, the long-time coarse-grained molecular dynamics (MD) simulation method is used to investigate the hemolytic toxicity of a typical hemolytic © XXXX American Chemical Society

Special Issue: Drug Metabolism and Toxicology Received: May 25, 2019 Published: June 24, 2019 A

DOI: 10.1021/acs.jmedchem.9b00853 J. Med. Chem. XXXX, XXX, XXX−XXX

Journal of Medicinal Chemistry

Article

on the erythrocytes from different species (e.g., rat, mouse, and rabbit). Hence the more attractive human hemolytic toxicity prediction was never probed before. Furthermore, all these three works focus on only one specific series of compounds (linear peptides with the natural residues or saponins) and thereby the corresponding prediction models are not applicable for more diverse small molecules. Therefore, it is quite necessary to develop a machine-learning based model to predict the human hemolytic toxicity for the diverse small molecules. In this work, we curate a data set of 337 hemolytic and 424 nonhemolytic small molecules, which are all experimentally measured on the human erythrocytes and manually compiled from various literature. On the basis of this new data set, we employ K-nearest neighbors (KNN),17 support vector machine (SVM),18 random forest (RF),19 gradient boosting machine (GBM),20 and extended-connectivity fingerprint (ECFP)21 to build the classification models for the human hemolytic toxicity prediction of small molecules, which are further deployed in our user-friendly program “e-Hemolytic-Classification” for the convenience of experimental medicinal chemists or toxicologists. Furthermore, we make full use of genetic algorithm (GA) and machine-learning based model to optimize the human hemolytic toxicity of small molecule to derive “optimal virtual fingerprints (OVFs)” with the desired hemolytic/nonhemolytic property, which is also implemented in our “e-Hemolytic-Classification” program.

Figure 1. Scatter plot of cLogP vs MW for our curated human hemolytic toxicity data set including 337 hemolytic (red spheres) and 424 nonhemolytic (blue spheres) compounds.

gives a similar conclusion. Furthermore, the distribution coefficient log D (PH = 7.4), a cell-permeability relevant descriptor, is computed by JCHEM (https://chemaxon.com) and used to tentatively evaluate the relation between the hemolytic toxicity and cell permeability for the small molecules because in one of the popular hemolytic mechanisms, it is assumed that the compound’s penetration into the membrane is very important,10 which leads to a hypothesis that the hemolytic compounds are probably prone to be cell permeable and the nonhemolytic compounds may tend to be cell nonpermeable. To test this hypothesis, the scatter plot of log D (PH = 7.4) vs MW is plotted. SI, Figure S3 illustrates that the majority of hemolytic/nonhemolytic compounds are cell permeable (log D > 0) from the computational perspective, which suggests that the human hemolytic toxicity of small molecules cannot be easily distinguished by the predicted cell permeability and MW, and our previous hypothesis is only partially correct. Moreover, three aforementioned basic chemotypes (“general small molecule”, “saponin”, “small (cyclic) peptide”) are mapped to the scatter plot of log D (PH = 7.4) vs MW to further inspect the effect of different chemotypes on the relation between the hemolytic toxicity and cell permeability. SI, Figure S4 displays that no sharp distinction is observed between the hemolytic and nonhemolytic compounds for each chemotype. Therefore, the simple descriptors such as log P, log D, and MW are not very intuitive for the classification of hemolytic/nonhemolytic compounds. However, ECFP-based similarity matrix (Figure 2) illustrates that the overall pairwise Tanimoto similarities between the hemolytic and nonhemolytic compounds are quite low and the mean of pairwise Tanimoto similarities is only 0.071 by averaging over the entire matrix. This similarity matrix indicates that ECFP fingerprint may be a good molecular descriptor for the classification of hemolytic and nonhemolytic compounds, and also implies that the compounds in our curated data set are structurally diverse. To further analyze which chemotype in our data set can be best distinguished by ECFP based descriptor, the whole similarity- matrix (Figure 2) has already been sorted according to three aforementioned chemotypes. Clearly, the submatrices highlighted in the green and purple dotted rectangles in Figure 2 demonstrate that the

■

RESULTS AND DISCUSSION Chemical Space of Our Curated Human Hemolytic Toxicity Data Set. Our curated data set contains 337 hemolytic and 424 nonhemolytic compounds, and it is the first publicly available data set on the hemolytic toxicity of small molecules that are fully evaluated on the human erythrocytes in the experimental hemolysis assays. In this human hemolytic toxicity data set, there are three basic chemotypes: “general small molecule”, “saponin”, and “small (cyclic) peptide”, and the distribution of each chemotype in our data set is displayed in Supporting Information (SI), Figure S1. Obviously, the predominant chemotype in our data set is “general small molecule” (187 hemolytic and 371 nonhemolytic compounds), thus the prediction model trained from our data set is more applicable to this chemotype. For the convenience of users, the SMILES string, category (hemolytic/nonhemolytic), chemotype, possible hemolytic mechanism, and corresponding reference for each compound in our data set are summarized in SI, Tables S1 and S2. It should be noted that most of original works only reported the hemolytic toxicity of the compounds without further exploration of their hemolytic mechanisms, albeit several possible hemolytic mechanisms were proposed in the literature.10,22 To roughly examine the difference of chemical spaces between the hemolytic and nonhemolytic compounds in our curated data set, octanol−water partition coefficient (log P) and molecular weight (MW) are calculated by both OpenBabel (v2.4)23 and XLogP (v2)24 programs. The scatter plot of cLogP vs MW (Figure 1) shows that the distributions for the hemolytic and nonhemolytic compounds are very similar, and it also displays that most of the compounds are located in the region (0 ≤ cLogP ≤ 10 and 0 < MW ≤ 1500). Hence, direct discrimination between the hemolytic and nonhemolytic compounds in terms of cLogP and MW does not work well. Additionally, the scatter plot of XLogP vs MW (SI, Figure S2) B



Article

Figure 3. Visualization of our curated human hemolytic toxicity data set by t-distributed stochastic neighbor embedding (t-SNE); t-SNE1 and t-SNE2 are the two dimensions reduced from 2048 dimensions of 2048bit-ECFP6. Most of the compounds are located in the large red rectangle. The minority compounds, which are narrowly concentrated in the two smaller rectangles, are saponins (blue rectangle) and small (cyclic) peptides (black rectangle). Green and orange spheres denote the hemolytic/nonhemolytic compounds, respectively.

Figure 2. Tanimoto similarity matrix for the hemolytic vs nonhemolytic compounds. Similarity is calculated based on 2048bitECFP6 fingerprint with “e-Hemolytic-Classification” program. For the nonhemolytic compounds, 1−371, 372−403, and 404−424 stand for the general small molecules, saponins, and small (cyclic) peptides, respectively (X-axis). For the hemolytic compounds, 1−187, 188− 271, and 272−337 refer to the general small molecules, saponins, and small (cyclic) peptides, respectively (Y-axis). Green dotted rectangle is the submatrix of similarities between the hemolytic and nonhemolytic saponins. Purple dotted rectangle is the submatrix of similarities between the hemolytic and nonhemolytic small (cyclic) peptides. The overall pairwise Tanimoto similarities between the hemolytic and nonhemolytic compounds are quite low, with the mean of 0.071 by averaging over the entire matrix.

of chemical space of compound database, are not able to well discriminate between the hemolytic and nonhemolytic compounds, and further analyses of the chemical space of our curated data set by using ECFP-based similarity matrix and t-SNE indicate that ECFP could be an informative descriptor to classify the hemolytic/nonhemolytic compounds, especially for those compounds with the main chemotype “general small molecule”. The Overall Performance of Human Hemolytic Toxicity Classification Model for the Small Molecules. The intensive model training and internal validation is conducted on the cross-validation data set (Dataset-CV defined in Experimental Section) with 5-fold cross-validation, and the systematic model evaluation is performed on the holdout test set (Dataset-Holdout-Test defined in Experimental Section) without any overlapping with Dataset-CV. Consequently, 1216 hemolytic/nonhemolytic classification models (M0001−M1216 in SI, Table S6) are harvested by taking account of all the combinations of four machine-learning methods (KNN, SVM, RF, and GBM), four ECFPs (1024bitECFP4, 2048bit-ECFP4, 1024bit-ECFP6, and 2048bitECFP6), four feature-number options (128, 256, 512, and full features), and 19 random data-partition schemes. The scatter plot of ΔF1-score vs F1-score (hold-out test set) for all 1216 models (Figure 4) manifests that F1-score predominantly varies from 0.82 to 0.94, and ΔF1-score (referring to |F1score(hold-out test set) − F1-score(cross-validation)|) for the majority of our models is lower than 0.06, revealing that most of our models exhibit the good performance without noticeable overfitting or underfitting from this perspective. Furthermore, Y-randomization tests for all these 1216 models are carried out in accordance with the procedure in SI, Text S6 to achieve “1216 models after applying Y-randomization test” (SI, Table S8 and Figure S10). The scatter plot of F1-score (hold-out test set) vs MCC (hold-out test set) showcases that our previous models without Y-randomization are quite robust, which can be supported by the substantial decline of the model performances after exerting Y-randomization (SI, Figure S10).

overall similarities between hemolytic and nonhemolytic compounds with the chemotype “general small molecule” are notably lower than the correspondence with the chemotype “saponin” or “small (cyclic) peptide”. Hence, ECFP is more informative for the classification of hemolytic and nonhemolytic compounds with the main chemotype “general small molecule”. To further visualize this human hemolytic toxicity data set, the unsupervised and nonlinear t-distributed Stochastic Neighbor Embedding (t-SNE) method is adopted to transform the similarities between data points to the joint probabilities and minimize the Kullback−Leibler (KL) divergence between the joint probabilities of the low-dimensional embedding and high-dimensional data, which is conducted based on 2048bitECFP6 by the Sklearn Python library. According to Figure 3, most of the compounds in our curated data set are widely scattered in the large red rectangle, albeit the minority compounds are narrowly concentrated in the small blue or black rectangles, which belong to the saponins (blue rectangle) or small (cyclic) peptides (black rectangle). In addition, the green and orange spheres in the red rectangle are largely separated, whereas they are significantly overlapped in the blue and black dotted rectangles, which also indicates that the hemolytic/nonhemolytic compounds with the chemotype “general small molecule” is less challenging to be distinguished relative to the counterparts with the chemotypes “saponin” and “small (cyclic) peptide”. Therefore, the results from the analyses of t-SNE and ECFP-based similarity matrix are quite consistent. Briefly speaking, the simple descriptors, such as log P, log D, and MW, which are commonly used for the rough estimation C



Article

ΔF1-score and F1-score (hold-out test set) and are robust according to Y-randomization test. Different Factors Impacting on the Model Performance of Human Hemolytic Toxicity Classification for Small Molecules. Observed from Figure 4, the apparently narrower distribution of orange spheres suggests that different random data-partition schemes have the dramatic effect on the model performance. Therefore, 64 average models (SI, Table S7), which are derived from 1216 individual models (SI, Table S6) by averaging over 19 data-splitting schemes, are more reasonable to fairly assess the different factors (ECFPs, machine-learning methods, and feature-number options) impacting on the model performance. Overall comparison reveals that the model performances for 64 average models do not differ a lot according to the comparable F1-scores of those average models (Table 1). So as to further gauge such difference from the statistical perspective, a series of two-sample t-Tests is conducted for 4032 pairs of different average models, and the corresponding results are summarized in SI, Table S9 and plotted in Figure 5, demonstrating that 3824 of 4032 pairs have no significant differences based on the criterion p-value < 0.001. Hence, different machine-learning methods, ECFPs, or feature numbers do not have the evident influence on the model performance in most circumstances. In a nutshell, the individual classification models are very sensitive to the different random data-splitting schemes, while most of our average classification models are insensitive to the choices of machine-learning methods, ECFPs, and featurenumber options.

Figure 4. Scatter plot of ΔF1-score vs F1-score (hold-out test set) for 1216 individual (green spheres) and 64 average (orange spheres) classification models; F1-score(hold-out test set) is used to evaluate the performance of the classification models on the hold-out test set, meanwhile ΔF1-score, referring to |F1-score(hold-out test set) − F1score(cross-validation)|, is harnessed to inspect the potential overfitting/underfitting of the classification models. The 64 average models (orange spheres) are obtained from 1216 individual models (green spheres) by averaging over 19 different random data-partition schemes.

In short, 1216 models are built to classify hemolytic/ nonhemolytic compounds by considering different machinelearning methods, ECFPs, feature numbers, and data-splitting schemes. The majority of our models are promising in terms of

Table 1. Different Feature Numbers, ECFPs, and Machine-Learning Methods Impacting on the Model Performance of 64 Average Models (SI, Table S3) with F1-Score (Hold-out Test Set) As the Performance Indicatora full features KNN

SVM

GBM

RF

1024bit-ECFP4 2048bit-ECFP4 1024bit-ECFP6 2048bit-ECFP6

0.877(0.023) 0.875(0.023) 0.870(0.030) 0.868(0.030)

0.877(0.028) 0.877(0.024) 0.878(0.029) 0.881(0.028) 512 features

0.870(0.026) 0.875(0.023) 0.884(0.026) 0.880(0.027)

0.882(0.025) 0.880(0.027) 0.891(0.023) 0.890(0.022)

KNN

SVM

GBM

RF


0.876(0.021) 0.881(0.023) 0.874(0.025) 0.874(0.028)

0.875(0.027) 0.876(0.024) 0.881(0.027) 0.879(0.022) 256 features

0.876(0.027) 0.874(0.025) 0.884(0.027) 0.886(0.024)

0.884(0.023) 0.881(0.024) 0.887(0.027) 0.890(0.020)

KNN

SVM

GBM

RF


0.877(0.020) 0.876(0.023) 0.876(0.028) 0.871(0.028)

0.871(0.029) 0.871(0.023) 0.884(0.026) 0.874(0.023) 128 features

0.875(0.026) 0.876(0.024) 0.883(0.024) 0.882(0.025)

0.880(0.022) 0.887(0.022) 0.890(0.022) 0.887(0.028)

KNN

SVM

GBM

RF

0.886(0.024) 0.878(0.027) 0.866(0.032) 0.865(0.031)

0.873(0.027) 0.867(0.024) 0.875(0.028) 0.871(0.029)

0.874(0.025) 0.870(0.025) 0.878(0.026) 0.871(0.027)

0.882(0.026) 0.878(0.027) 0.879(0.027) 0.877(0.029)

1024bit-ECFP4 2048bit-ECFP4 1024bit-ECFP6 2048bit-ECFP6 a

Notes: (1) The number in the table is the mean (standard deviation) by averaging over 19 data-splitting schemes. (2) F1-score (hold-out test set) refers to the F1-score on the hold-out test set. D



Article

set. Therefore, in this work, two representative models (a typical consensus model and a fast individual model) are proposed for the users with different applications. To further assess the performances of these two typical hemolytic toxicity classification models (CM and FIM), an external test set “Dataset-External-Test” (50 hemolytic and 32 nonhemolytic compounds evaluated on the human red cells, SI, Table S3) is newly curated from the literature (SI, Table S4) after April 2018, while our previous data curation was completed before April 2018 (SI, Table S2) and led to the human hemolytic toxicity data set (Dataset-Whole, 337 hemolytic and 424 nonhemolytic compounds in SI, Table S1), which is split into the cross-validation data set “DatasetCV” for model training and the hold-out test set “DatasetHoldout-Test” for model evaluation. It should be mentioned that the data curation procedure is the same for DatasetExternal-Test and the previous human hemolytic toxicity data set (Dataset-Whole including Dataset-CV and DatasetHoldout-Test), However, this Dataset-External-Test has no overlapped samples with Dataset-Whole. According to Table 3, the consensus model (CM) provides the accuracy (0.829), precision (0.909), specificity (0.875), sensitivity (0.800), F1score (0.851), and MCC (0.660), respectively, on DatasetExternal-Test, which is comparable with the model performance on Dataset-Holdout-Test in terms of F1-score (0.851 vs 0.909), precision (0.909 vs 0.886), and specificity (0.875 vs 0.902). Nevertheless, it should be pointed out that MCC (0.660) is lower than its counterpart (0.835) on DatasetHoldout-Test. Meanwhile, the fast individual model (FIM) gives the accuracy (0.659), precision (0.867), specificity (0.875), sensitivity (0.520), F1-score (0.650), and MCC (0.400), respectively, on Dataset-External-Test (Table 3), which is inferior to the model performance on DatasetHoldout-Test except the precision (0.867 vs 0.842) and specificity (0.875 vs 0.857). Although the overall performance on Dataset-External-Test is decreased, our models (CM and FIM) still retain the good capability of picking out the nonhemolytic compounds because 28 of 32 nonhemolytic compounds can be classified correctly by both models (Table 3). On the basis of these comparisons, our consensus model (CM) still exhibits the promising performance on DatasetExternal-Test and is more robust than the fast individual model (FIM). In short, a reliable consensus model (CM) and a fast individual model (FIM) are derived to predict the human hemolytic toxicity of small molecules. Both models demonstrate the good precision and specificity on two independent test sets (Dataset-Holdout-Test and Dataset-External-Test) and CM affords the robust performance on both test sets. Moreover, the computationally slower consensus model (CM) is more reliable than the fast individual model (FIM).

Figure 5. Model performance comparison matrix based on the F1scores (hold-out test set) of 64 average models (SI, Table S7) by twosample t-test analysis; 64 blue squares highlight that the same average models will not be considered in the two-sample t test. The 3824 gray squares manifest that these two different average models have no significant difference in the model performance according to p-value < 0.001, while 208 white squares show that these two different average models have significant difference in the model performance based on the criterion (p-value < 0.001). The original data of each element in this matrix is tabulated in SI, Table S9.

Performance Evaluation of the Typical Models on the Hold-Out Test Set and External Test set. As we know, harnessing all 1216 hemolytic toxicity prediction models is not very pragmatic for the users in the practical prediction, thus a typical consensus model is constructed by choosing 19 best individual models from SI, Table S6 according to the highest F1-score in each data-partition scheme. Seen from Table 2, this consensus model (CM) offers the accuracy (0.917 ± 0.019), precision (0.886 ± 0.035), specificity (0.902 ± 0.036), sensitivity (0.935 ± 0.027), F1-score (0.909 ± 0.019), and MCC (0.835 ± 0.036) on the hold-out test set by averaging over 19 data-splitting schemes to alleviate the performance bias resulted from the random data-partition (SI, Table S10). In addition, a fast individual model (FIM in Table 2), which is actually one constituent model (M0267 in SI, Table S10) of CM with the lowest ΔF1-score (0.007), is selected for the rapid screening of large compound database, because the screening of the huge database with the consensus model (CM) is still time-consuming due to the fact that 19 calculations based on 19 individual models are required in each prediction to achieve the final consensus result. This fast model is trained with GBM method and full features of 2048bit-ECFP6, and affords the accuracy (0.901), precision (0.842), specificity (0.857), sensitivity (0.955), F1-score (0.895), and MCC (0.807), respectively, on the hold-out test

Table 2. Performance of a Typical Consensus Model and a Fast Individual Model on the Hold-out Test Set and in the CrossValidationa model

accuracy (test set)

precision (test set)

specificity (test set)

sensitivity (test set)

F1-score (test set)

MCC (test set)

F1-score (CV)

ΔF1-score

CM FIM

0.917(0.019) 0.901

0.886(0.035) 0.842

0.902(0.036) 0.857

0.935(0.027) 0.955

0.909(0.019) 0.895

0.835(0.036) 0.807

0.883(0.012) 0.888

0.029(0.023) 0.007

a

Notes: (1) the number in each parentheses is the standard deviation (SD), which is obtained according to 19 random data-splitting schemes. (2) ΔF1-score referring to |F1-score(test set) − F1-score(CV)| is utilized to monitor the potential overfitting/underfitting of the classification model. (3) “CM”, “FIM”, “MCC”, and “CV” are short for “consensus model”, “fast individual model”, “Matthews correlation coefficient”, and “crossvalidation” respectively. E



Article

Table 3. Performance of a Typical Consensus Model and a Fast Individual Model on the External Test Set (50 Hemolytic and 32 Nonhemolytic Compounds Experimentally Evaluated on the Human Red Cells, SI, Table S3)a model

TP#

TN#

FP#

FN#

accuracy (test set)

precision (test set)

specificity (test set)

sensitivity (test set)

F1-score (test set)

MCC (test set)

CM FIM

40 26

28 28

4 4

10 24

0.829 0.659

0.909 0.867

0.875 0.875

0.800 0.520

0.851 0.650

0.660 0.400

Notes: (1) “TP#”, “TN#”, “FP#”, “FN#”, “CM”, “FIM”, and “MCC” are short for “the number of true hemolytic compounds”, “the number of true nonhemolytic compounds”, “the number of false hemolytic compounds”, “the number of false non-hemolytic compounds”, “consensus model”, “fast individual model”, and “Matthews correlation coefficient” respectively. (2) The external test set (50 hemolytic and 32 nonhemolytic compounds, SI, Table S3) is newly curated from the literature (SI, Table S4) after April 2018, while our previous data curation was completed before April 2018 (SI, Table S2), which leads to 337 hemolytic and 424 nonhemolytic compounds (SI, Table S1). It should be mentioned that this external test set has no overlapped samples with our previous human hemolytic toxicity data set.

a

Figure 6. Main features of “e-Hemolytic-Classification” program. (1) Prediction of hemolytic/nonhemolytic compound; (2) database screening; (3) visualization and enquiry of our curated human hemolytic toxicity data set; (4) automatic inspection of applicability domain of our model for the given compound; (5) visualization of important feature in the 3D viewer; (6) calculation and visualization of fingerprint in the 3D viewer; (7) model customization based on the user’s in-house data set; (8) genetic algorithm (GA) guided hemolytic toxicity optimization; (9) “optimal virtual fingerprint (OVF)” extraction from the GA evolution trajectory; (10) compound mapping of OVF in the chemical database.

Deployment of the Typical Models in the UserFriendly Software for the Practical Application. To enable the nonexpert users without the computational background to perform the hemolytic toxicity prediction for the small molecules, the aforementioned models (CM and FIM) and relevant nontrivial procedures such as feature generation and selection are fully encapsulated and deployed in an easy-to-use program “e-Hemolytic-Classification” (Figure 6) for the automatic prediction of hemolytic/nonhemolytic compounds, which is developed based on our previous program “e-Bitter” for the prediction of bitter compounds.25 Besides this automatic prediction, “e-Hemolytic-Classification” program also furnishes some useful auxiliary functions: (1) Similarity search in our curated data set, (2) virtual screening of the compound database, (3) generation and visualization of ECFP fingerprint, (4) model interpretation by the synchronous visualization of the feature importance and fingerprint-bit “1” in the context of 3D structure, (5) automatic inspection of the applicability domain of our model for the given molecule, and (6) model customization based on the user’s in-house data set. The full package including the program, its corresponding

manual, and tutorials can be freely accessible from the Dropbox public shared-folder with the link as follows (https://www.dropbox.com/sh/ewds3mvn7f39co5/ AABgGfgsVd7Z3w9wQn_0tBeua?dl=0). In sum, two typical models (a consensus model and a fast individual model) have been deployed in an easy-to-use program for the automatic prediction via the simple mouse clicks. Intuitive Interpretation of Human Hemolytic Toxicity Classification Model for the Small Molecules. Our human hemolytic toxicity classification model can be intuitively interpreted according to the feature importance (FI), which can be conveniently visualized in our “e-Hemolytic-Classification” program with the implementation detail in SI, Text S9. For an illustrative purpose, only three most important structural features (Figure 7 and SI, Figures S11−S12) with three largest FIs are taken as the examples. The fingerprint bits for these three structural features are 1613-bit (Figure 7), 173bit (SI, Figure S11), and 365-bit (SI, Figure S12), respectively, their corresponding FIs are 0.024966, 0.023243, and 0.018565, respectively, and all these three structural features are strongly F



Article

Figure 7. Visualization of the 3D structural feature for 1613-bit with the largest feature importance (FI = 0.024966) in “e-Hemolytic-Classification” program.

toxicity to produce “optimal virtual fingerprints (OVFs)” that can satisfy the specific condition (e.g., the target compounds should be hemolytic and within the applicability domain of the aforementioned prediction model), subsequently, we can retrieve the closest real compounds from the chemical database of interest via similarity search based on OVFs that have acquired the enough structural features (or fingerprint bits) to induce the target property (e.g., hemolytic toxicity) via GA evolution, although OVFs cannot be directly mapped back to the corresponding chemical structures. It is worthy to note that GA guided optimization based on the molecular fingerprints can explore the diverse chemical spaces due to the nature of the binary encoding of chromosomes, and the resultant OVFs are independent of chemical databases, albeit the huge and diverse compound databases will have more chance to match the OVFs to elicit the potential real hemolytic/nonhemolytic compounds with the higher similarity to OVFs. If the closest compound retrieved from the compound database bears the low Tanimoto similarity to the specific OVF of enquiry, it means that this compound only contains the partial structural features of OVF. For the further confirmation, the retrieved compound will be subjected to the human hemolytic toxicity prediction based on our previous model. Therefore, the huge compound database for the full match of OVF is not required, and the small compound database is also applicable, if it includes the compound sharing certain partial structural features with OVFs. To demonstrate whether OVFs are practical to recover some compounds that are highly similar to (or exactly the same as) the existing hemolytic compounds in our collected data set,

relevant to the aliphatic chain, Thus, the aliphatic chain contributes to the classification of hemolytic/nonhemolytic compounds. It is worth mentioning that if necessary, our program can highlight all the important structural features, the associated fingerprint-bits, and the pertinent FIs in the interactive and synchronous manner. Simply put, the important structural features based on FIs may be tentatively considered as the simple hemolytic toxicity alerts for the experimental chemists. However, these structural alerts should be used with caution due to the fact that the hemolytic toxicity of small molecule is usually induced by the synergic effect of multiple structural features, rather than the single effect of one important structural feature. Genetic Algorithm Guided Optimization of Human Hemolytic Toxicity for the Small Molecules. Genetic algorithm is a global search strategy inspired by the concept from Darwin’s theory of evolution and can drive a population of randomly generated chromosomes to be adaptively evolved to the best chromosomes (referring to the best solutions for the problem of interest) with the optimal fitness score by applying the genetic operators (e.g., selection, mutation, and crossover) in the iterative manner. Thus, GA has been widely used in the computer-aided drug design such as molecule docking.26 However, the optimization of human hemolytic toxicity by GA was never attempted before, and is carried out in this work by the seamless fusion of the ECFP fingerprint, machine-learning based prediction model, and GA. The whole GA optimization protocol is concisely illustrated in Figure 8. More specifically, based on the machine-learning based prediction model as mentioned above, we implement the genetic algorithm (GA)27 guided optimization of hemolytic G



Article

Figure 8. Basic flowchart for the genetic algorithm (GA) guided hemolytic toxicity optimization to demonstrate how to derive the optimal virtual fingerprint (OVF) and retrieve the real compound from the external compound database based on the OVF. The initial chromosomes with the length of 2048 can be automatically initialized by the randomly generated “1/0” bit strings or seeded from 2048bit-ECFP6 fingerprint(s) of the given molecule(s). It should be noted that in our test the initial 100 chromosomes are randomly generated to make sure that these chromosomes are highly dissimilar to the hemolytic compounds in our data set because the fitness function is set to “hemolytic/AD” (eq 3 in the Experimental Section) in this test. Moreover, “selection criterion” is dependent on the users’ applications (e.g., if the users would like to retrieve the OVFs standing for the novel compounds, the OVFs with low similarity with our human hemolytic toxicity data set are retained). “Convergent criterion” refers to the preset number of evolution generation. In this test, the evolution generation number is set to 1500, and the “average fitness” of each generation along the whole GA evolution is recorded for the convergence analysis of GA evolution.

GA optimization is performed with the fitness function “hemolytic/AD” (eq 3) that is defined in Experimental Section. In detail, a population of 100 chromosomes are used in the GA optimization, and each chromosome with the length of 2048 is randomly initialized to ensure that each initial structure denoted by each chromosome is highly dissimilar to the hemolytic compounds in our curated data set, which can be observed from Figure 9 with the low average similarity (0.103) for these 100 chromosomes. The number of evolution generation, mutation rate, and crossover rate are set to 1500, 0.3, and 0.7, respectively, which are the default values in our program. The fast individual model (FIM in Table 2) is adopted to compute the predicted probability of a compound being hemolytic, which is required in this fitness function (“hemolytic/AD” defined in eq 3). The entire GA evolution trajectory including all the chromosomes in each generation will be logged for the subsequent analysis. The average fitness of each generation along the whole GA evolution is monitored and recorded in Figure 10, which unambiguously illustrates that our GA evolution is convergent after about 900 generations. Therefore, our default parameters for the GA optimization are reasonable. Because of the intrinsic randomness of GA method, the whole GA run is repeated five times. Once completed, five GA evolution trajectories are combined to form 750000

chromosomes in total and are subjected to the elimination of the chromosomes that are duplicated, or nonhemolytic, or outside the applicability domain of our model, which refers to “selection criterion” in Figure 8. Finally, the remaining 249032 chromosomes are considered as OVFs. To vividly manifest whether some OVFs are similar to (or exactly the same as) the molecular fingerprints of our curated hemolytic compounds and dissimilar to the molecular fingerprints of our collected nonhemolytic compounds, the similarities between OVFs and their respective closest hemolytic/nonhemolytic compounds in our data set are calculated and plotted in Figure 11. Seen from Figure 11, it is encouraging that 12269 of 249032 OVFs located in the blue dotted rectangle are very close to the fingerprints of our curated hemolytic compounds and are remarkably different from our collected nonhemolytic compounds. It is worth mentioning that the OVF highlighted by blue sphere in Figure 11 is exactly the same as the fingerprint of one curated hemolytic compound with the similarity of 1.000. Therefore, our GA guided optimization is quite effective to well reproduce the molecular fingerprint with the hemolytic potential. More importantly, previous five runs of GA guided optimization can also generate the OVFs that represent the novel hemolytic compounds, which are highly dissimilar to the compounds in our curated data set and are of more interests to H



Article

Figure 9. Histogram of the similarities between 100 initial randomized chromosomes and their respective closest hemolytic compounds in our curated data set. The similarity refers to the Tanimoto similarity based on 2048bit-ECFP6. The average similarity for these 100 chromosomes is 0.103, which indicates that the initial chromosomes prepared for the GA optimization are highly dissimilar to our curated hemolytic compounds. It should be mentioned that “Similarity to Hemolytic Compound” in the X-axis refers to the Tanimoto similarity between the randomized chromosome and the closest single hemolytic compound in our curated data set.

Figure 11. Chemical space of 249032 OVFs from five runs of GA guided optimization with the fitness function (“hemolytic/AD” defined by eq 3 in the Experimental Section). X-axis refers to the similarity between OVF and its closest hemolytic compound in our curated data set, and Y-axis stands for the similarity between OVF and its closest nonhemolytic compound in our curated data set. The similarity is short for “Tanimoto similarity” computed based on 2048bit-ECFP6 with our “e-Hemolytic-Classification” program. In this scatter plot, 12269 OVFs in the blue dotted rectangle are highly similar to our curated hemolytic compounds, and one hemolytic OVF (blue sphere) is exactly the same as the fingerprint of one curated hemolytic compound with the similarity of 1.000, which suggests that our GA guided optimization can recover the existing hemolytic compound in our curated data set. 111217 OVFs in the red dotted rectangle are highly dissimilar to the compounds in our curated data set and represent the potential novel hemolytic compounds. One typical hemolytic OVF (purple sphere) in the red dotted rectangle is selected for the demonstration purpose and mapped to the ZINC compound (ZINC ID: 72326049 with the chemical structure shown above) via similarity search in the ZINC compound database.

Figure 10. Change of average fitness of each generation along the whole GA evolution (1500 generations). The curve levels off after about 900 generations, indicating that our GA evolution is convergent.

with 2048bit-ECFP6 fingerprint in our program, which returns a candidate compound (ZINC ID: 72326049 in Figure 11) with the highest Tanimoto similarity of 0.120. Although this retrieved compound has low similarity, it still shares 85 structural features with the hemolytic OVF and thereby still tends to be hemolytic. As expected, this retrieved compound is within the applicability domain of our hemolytic toxicity prediction model and predicted to be hemolytic (the predicted probability of being hemolytic: 0.962) according to the FIM model in our “e-Hemolytic-Classification” program (SI, Figure S13). Finally, this retrieved compound is also mapped to Figure 11 (orange sphere: X-coordinate = 0.344 and Ycoordinate = 0.230), indicating that this ZINC compound (ZINC ID: 72326049) is still very dissimilar to our curated hemolytic/nonhemolytic compounds and could be a novel hemolytic compound from the computational perspective. Therefore, the hemolytic OVF learned from GA evolution inclines to accumulate the hemolytic structural features as many as possible, and the partial match of structural features of hemolytic OVF in the common chemical database can help us discover the potential novel hemolytic compound, albeit

the experimental chemists. As illustrated in Figure 11, 111217 of 249032 OVFs highlighted in the red dotted rectangle are significantly different from the fingerprints of our curated hemolytic and nonhemolytic compounds. Similarity search of these OVFs in the chemical database can retrieve the potential novel hemolytic compounds whose fingerprints are the closest to these OVFs. For the demonstration purpose, one typical hemolytic OVF highlighted by the purple sphere (Xcoordinate = 0.300 and Y-coordinate = 0.105) in Figure 11 has been selected because this OVF has the lowest similarity to the closest nonhemolytic compound in our curated data set when it possesses the similarity of 0.300 relative to the closest hemolytic compound in our collected data set, which suggests that this typical hemolytic OVF is highly different from our curated nonhemolytic and hemolytic compounds. Subsequently, similarity search of this hemolytic OVF (purple sphere) in the extensively used ZINC database28 is performed I



Article

In brief, we adopt GA and machine-learning based prediction model to optimize the hemolytic toxicity to derive the hemolytic OVFs, some of which are shown to enable the recovery of some hemolytic compounds that are very similar to (or the same as) the compounds in our curated data set, and some of which are manifested to capture the potential novel hemolytic compounds that are quite different from our curated compounds. The Limitation and Prospect of This GA Guided Optimization Protocol for Human Hemolytic Toxicity. The pros and cons of this GA optimization protocol are also discussed to help the users determine whether this protocol is applicable to their own projects. The advantages of this protocol are described as below: (1) This protocol can explore the broad chemical space based on the binary molecular fingerprint and consequently can yield the diverse and novel OVFs, which are independent of the chemical database and can be used to search any given chemical database specified by the users. (2) This protocol is quite general and can be applicable to other AMDET end points (Figure 13) in our “eHemolytic-Classification” program. Meanwhile, the disadvantages of this protocol are also listed as follows: (1) The OVFs cannot be directly transformed back to the corresponding chemical structures. Hence, similarity search in the chemical database of the users’ interests is required to retrieve the real compounds whose fingerprints are the closest to the OVFs. (2) The inherent randomness in the GA optimization will give rise to different OVFs for each GA run, thus the multiple GA runs are highly recommended and all the GA evolution trajectories from different GA runs are suggested to be combined for the subsequent analysis, but undoubtedly it will introduce some extra computational burden. Therefore, this general protocol should be used with caution. To further enhance the reliability of this GA guided optimization protocol, two major aspects could be strengthened in the future. First, the improvement of hemolytic toxicity prediction model by training on more extended and diverse high-quality human hemolytic toxicity data set may further increase the accuracy of fitness calculation in the GA optimization, because the derived optimal virtual fingerprints (OVFs) could be affected by the curated data set in the step of “Fitness calculation” (Figure 8) with the following details: (1) The cross-validation data set (Dataset-CV) has very important influence on the prediction model that will be used to compute the predicted probability of being hemolytic/nonhemolytic for the fitness calculation of each chromosome. (2) The whole hemolytic toxicity data set is utilized to define the applicability domain of our prediction model by “average similarity” (defined in SI, Text S8), which is included in the fitness calculation if the fitness function “hemolytic/AD” or “nonhemolytic/AD” is specified. As such, we will continue to expand our human hemolytic toxicity data set in the future. Second, the larger external compound database may benefit the compound mapping of OVF to retrieve the novel real compound with the higher similarity. Overall, the merits and limitations of the GA guided optimization protocol are addressed to inform the users of its applicability, and the future endeavors on this protocol are discussed as well.

similarity search of OVF in the huge chemical space may afford the compounds containing more hemolytic structural features from the hemolytic OVF. To make this GA guided optimization protocol more practical to the users, an automatic function has been developed in our “e-Hemolytic-Classification” program to derive the OVFs with the specific property, and the molecular fingerprints (2048bit-ECFP6 fingerprints) of ZINC database28 (around 20 million compounds) with the vendor information (http://zinc.docking.org/subsets/all-purchasable) are also integrated into our program to retrieve the closest real compounds for the OVFs via similarity search. In our current implementation, there are two main functions: (1) GA optimization with the fitness function (“hemolytic”, “hemolytic/AD”, “non-hemolytic”, or “non-hemolytic/AD” defined in the Experimental Section) generates the diverse and new OVFs that denote the novel hemolytic/nonhemolytic compounds. (2) GA optimization with the fitness function (“hemolytic/ SIM” or “non-hemolytic/SIM” defined in Experimental Section) produces the OVFs that are similar to 2048bitECFP6 fingerprint(s) of the given compound(s) with certain similarity threshold. For the purpose of concise demonstration, the snapshot of these functions is shown in Figure 12, while the detailed usage and tutorial are provided in the program manual.

Figure 12. Parameter settings for the genetic algorithm (GA) guided hemolytic toxicity optimization. The default fitness function is “toxic/ AD”, which is equivalent to “hemolytic/AD” defined by eq 3 in the Experimental Section, and the default parameters for the evolution generation number, population size, mutation rate, and crossover rate are 1500, 100, 0.3, and 0.7, respectively. It should be noted that “toxic” actually refers to “hemolytic” in this test, because the “Model source” is set to “internally integrated hemolytic toxicity model” in this panel; however, the option “toxic” is more general considering that this panel can be also applicable to the user-customized models for other toxicity end points.

■

CONCLUSIONS In this work, we manually curate the hemolytic toxicity data set including 337 hemolytic and 424 nonhemolytic small J



Article

Figure 13. Main interfaces for “Model Builder” (A) and “Model Loader” (B) that furnish the general functions to build and load the usercustomized model for other end points. The detailed usage about the model customization and loading has been provided in the program manual.

(or exactly the same as) the existing hemolytic compounds in our curated data set. More importantly, this GA guided optimization has also been manifested to produce the diverse “OVFs” representing the potential novel hemolytic compounds, which are dissimilar to the compounds in our curated data set and can be retrieved via similarity search in the chemical database of interest. To our knowledge, we present a promising computational model and a practical software to predict/optimize the small molecules’ hemolytic toxicity on the human red blood cells in light of our protocol by fusing ML and GA, which is general and applicable to other end points.

molecules experimentally assayed on the human red blood cells. Based on this new data set, we propose two typical machine-learning based classification models: a robust consensus model (CM) providing the accuracy, precision, specificity, sensitivity, F1-score, and MCC of 0.917 ± 0.019, 0.886 ± 0.035, 0.902 ± 0.036, 0.935 ± 0.027, 0.909 ± 0.019, and 0.835 ± 0.036, respectively, on the hold-out test set, and a fast individual model (FIM) affording the accuracy, precision, specificity, sensitivity, F1-score, and MCC of 0.901, 0.842, 0.857, 0.955, 0.895, 0.895, and 0.807, respectively, on the holdout test set. Moreover, the evaluation of CM and FIM on the external test set indicates that CM still exhibits the promising performance on the external test set, and both models (CM and FIM) well retain the good precision and specificity on the external test set. For the convenience of users, both models have been implicitly integrated into a stand-alone software “eHemolytic-Classification” with the versatile functions. Furthermore, we make full use of genetic algorithm (GA) and machine-learning (ML) based prediction model to optimize the hemolytic toxicity based on the molecular fingerprint to derive “optimal virtual fingerprints (OVFs)” that have been demonstrated to capture some compounds highly similar to

■

EXPERIMENTAL SECTION

Data Curation of Hemolytic/Nonhemolytic Compounds Assayed on the Human Erythrocytes. Hemolytic/nonhemolytic compounds evaluated on the human red cells were manually retrieved from a large volume of literature (published before April 2018) that was collected by searching the keywords “hemolytic” and “haemolytic” in the “web of science” and “google scholar”. In the literature, the reported hemolytic or nonhemolytic compounds are generally defined by HD50 (dose causing 50% of maximum hemolysis) threshold with the concentration of 200 μM.29−33 If HD50 of the K



Article

compound is larger than 200 μM, this compound is usually considered to be nonhemolytic, otherwise, it is deemed hemolytic. According to this criterion, 337 hemolytic and 424 nonhemolytic compounds were compiled from 109 articles (SI, Table S2) containing the hemolytic toxicity measured on the human red blood cells after the redundant/conflicted samples were eliminated. Subsequently, this data set (337 hemolytic and 424 nonhemolytic compounds) was fully integrated into “e-Hemolytic-Classification” program, which provided the flexible options for the intuitive visualization and full export of this data set. Moreover, the whole data set including the SMILES strings, the categories (hemolytic/ nonhemolytic), the chemotypes (“general small molecule”, “saponin”, or “small (cyclic) peptide”), the proposed hemolytic mechanism, and reference is also given in the Supporting Information (Tables S1−S2 and “human-hemolytic-toxicity-dataset.csv”). Finally, all the compounds except the quaternary ammonium compounds (the positive charge on the quaternary nitrogen atom is explicitly assigned in the original literature) are neutralized for the subsequent calculation of molecular descriptors, which will be used in the following model training and testing. Data-Splitting Schemes for the Model Training and Testing. To build and test the human hemolytic toxicity classification model for the small molecules, the entire data set (337 hemolytic and 424 nonhemolytic compounds) is randomly split into two parts: the data set for cross-validation (Dataset-CV) and the hold-out test set for model validation (Dataset-Holdout-Test). The detailed data-splitting scheme is elaborated as follows: 80% of hemolytic compounds and 80% of nonhemolytic compounds, which are randomly selected from the whole data set (Dataset-Whole), are used to train the model with 5-fold cross-validation (CV) and denoted as Dataset-CV, while the rest of them are treated as the hold-out test set (Dataset-HoldoutTest). Furthermore, Dataset-CV is further randomly split into five chunks, while each chunk (except the last chunk) keeps the original ratio between the hemolytic and nonhemolytic compounds. Subsequently, one chunk is employed as the internal validation set (Dataset-Internal-Validation), and the remaining four chunks are combined to form the training set (Dataset-Training). This procedure is repeated five times, which leads to five sets of Dataset-Training and Dataset-Internal-Validation for the 5-fold cross-validation. At last, the whole data-splitting procedure as mentioned above will be repeated 19 times in order to reduce the possible bias from different random data-splitting schemes. For the further validation of human hemolytic toxicity classification models, an external validation set (Dataset-External-Test in SI, Table S3) is newly curated from the literature (SI, Table S4) after April 2018 and has no duplicated samples with our previously collected hemolytic toxicity data set (Dataset-Whole in SI, Table S1), which was compiled from the literature (Table S2) before April 2018. Thus, this Dataset-External-Test will be more challenging for the model testing than the previous Dataset-Holdout-Test. The Protocol for the Building, Evaluation, and Deployment of Model. In the present study, K-nearest neighbors (KNN),17 support vector machine (SVM),18 random forest (RF),19 and gradient boosting machine (GBM)20 are combined with the molecular fingerprints to build the hemolytic/nonhemolytic classification model for the small molecules. The full protocol was developed in our previous work regarding to the prediction of bitter compounds,25 which can be straightforwardly adapted for this work. Therefore, we make a brief summary of this procedure as follows and address each detail in the Supporting Information (SI, Texts S1−S9). The whole protocol can be basically divided into nine steps. (1) Generating the extended-connectivity fingerprint (ECFP)21 based molecular descriptor, which is natively implemented in “e-HemolyticClassification” program for the intuitive visualization of ECFP bit in the context of 3D structure viewer. (2) Conducting the featureselection according to the feature importance (FI) obtained from the random forest (RF) approach. (3) Model-training by considering all the combinations of four machine-learning methods (KNN, SVM, RF, and GBM), four ECFPs (1024bit-ECFP4, 2048bit-ECFP4, 1024bitECFP6, and 2048bit-ECFP6), four feature-number options (128, 256,

512, and full features), and 19 data-splitting schemes to achieve 1216 models. (4) Model-evaluation with the standard metrics such as the accuracy, precision, specificity, sensitivity, F1-score, and Matthews correlation coefficient (MCC) measured on the hold-out test set (Dataset-Holdout-Test) and external test set (Dataset-External-Test). (5) Y-Randomization test34 of the models to check the modelreliability. (6) Construction of a representative consensus model (CM) and selection of a fast individual model (FIM). (7) Inspecting the applicability domain (AD) of our model to ensure its confident prediction.35 (8) Model interpretation with the feature importance (FI), which can be intuitively displayed in “e-Hemolytic-Classification” program. (9) Model-deployment in the graphic program “eHemolytic-Classification” for the users to forecast the human hemolytic toxicity of small molecules in the automatic fashion. Genetic Algorithm Guided Optimization of Human Hemolytic Toxicity for the Small Molecules. To discover more nonhemolytic compounds for the therapeutic purpose and hemolytic compounds used as the chemical probes to further cast light on the hemolytic mechanism, genetic algorithm (GA)27 is implemented to optimize the hemolytic toxicity based on our machine-learning based prediction model as mentioned above (Figure 8). The detail of GA optimization is given as follows: the chromosome with the length of 2048 denotes the molecular fingerprint of one compound. Each chromosome in the whole population can be automatically initialized in the random way or seeded from 2048bit-ECFP6 fingerprint(s) of the given molecule(s) of interest. The mutation operator is to randomly flip one fingerprint bit from “0” to “1” (or from “1” to “0”). The crossover operator is to swap two fingerprints fragments based on the random crossover point. The selection operator is to randomly select the chromosome based on the roulette wheel algorithm. The fitness function is defined by one of six options (“hemolytic”, “nonhemolytic”, “hemolytic/AD”, “non-hemolytic/AD”, “hemolytic/SIM”, or “non-hemolytic/SIM”), which are elaborated as follows: (1) “Hemolytic” and “non-hemolytic” represent that the fitness functions are the predicted probabilities of the compounds being hemolytic and nonhemolytic respectively, which are expressed by eqs 1−2. (2) “Hemolytic/AD” and “non-hemolytic/AD” mean that the fitness functions consider the applicability-domain (AD) of the given molecule besides the predicted probability, which are formulated by eqs 3−4. (3) “Hemolytic/SIM” and “non-hemolytic/SIM” denote that the fitness functions include the predicted probability and the similarity to the given molecule(s), which are described by eqs 5−6. Moreover, the population size, the number of evolution generation, mutation rate, and crossover rate are set to 100, 1500, 0.3, and 0.7, respectively, by default but can be adjusted by the users. After GA evolution is completed, all the generated chromosomes (or molecular fingerprints) that satisfy the target property with the appropriate threshold (referring to the step “selection criterion” in Figure 8), are extracted from the whole GA evolution trajectory and denoted as “optimal virtual fingerprints (OVFs)”. Nevertheless, the direct one-toone mapping of OVF to the chemical structure of compound is quite challenging due to the intrinsic mapping-complexity and unavoidable bit-collision in the fingerprint with the fixed-length, thus we harness the similarity search to retrieve the most similar compounds by enquiring the given compound database(s) with OVFs. Fitness(hemolytic) = p

(1)

Fitness(nonhemolytic) = 1.0 − p

(2)

Fitness(hemolytic/AD) = p + βSAD

(3)

Fitness(nonhemolytic/AD) = (1.0 − p) + βSAD

(4)

Fitness(hemolytic/SIM) = p + βScompound

(5)

Fitness(nonhemolytic/SIM) = (1.0 − p) + βScompound

(6)

where p is predicted probability of a compound being hemolytic; β is weighting factor (default value: 5.0) that can be customized by the users and is employed to balance the convergence of GA evolution L



■

and diversity of OVFs. If the weighting factor is smaller, GA evolution will produce more diverse OVFs but will require much longer generations to get the convergent result. When the weighting factor is zero, eqs 3−4 and eqs 5−6 are degraded to eqs 1−2. SAD: the “average similarity” to our curated hemolytic/nonhemolytic data set, which is used to define the applicability domain of our prediction model (detailed in SI, Text S8) and is included in the fitness function (eqs 3−4) to ensure that the GA generated solutions will be within the applicability domain of our model. Scompound: the similarity to the given compound(s) that is utilized to drive the GA evolution toward the specific compound(s) if the users would like to harvest the various OVFs similar to the given compound(s) of interest.

■

REFERENCES

(1) van de Waterbeemd, H.; Gifford, E. ADMET in silico modelling: towards prediction paradise? Nat. Rev. Drug Discovery 2003, 2, 192. (2) Norinder, U.; Bergström, C. A. S. Prediction of ADMET properties. ChemMedChem 2006, 1, 920−937. (3) Dearden, J. C. In silico prediction of ADMET properties: how far have we come? Expert Opin. Drug Metab. Toxicol. 2007, 3, 635− 639. (4) Gleeson, M. P.; Hersey, A.; Montanari, D.; Overington, J. Probing the links between in vitro potency, ADMET and physicochemical parameters. Nat. Rev. Drug Discovery 2011, 10, 197. (5) Moroy, G.; Martiny, V. Y.; Vayer, P.; Villoutreix, B. O.; Miteva, M. A. Toward in silico structure-based ADMET prediction in drug discovery. Drug Discovery Today 2012, 17, 44−55. (6) Cumming, J. G.; Davis, A. M.; Muresan, S.; Haeberlein, M.; Chen, H. Chemical predictive modelling to improve compound quality. Nat. Rev. Drug Discovery 2013, 12, 948. (7) Cheng, F.; Li, W.; Liu, G.; Tang, Y. In silico ADMET prediction: recent advances, current challenges and future trends. Curr. Top. Med. Chem. 2013, 13, 1273−1289. (8) Kirchmair, J.; Göller, A. H.; Lang, D.; Kunze, J.; Testa, B.; Wilson, I. D.; Glen, R. C.; Schneider, G. Predicting drug metabolism: experiment and/or computation? Nat. Rev. Drug Discovery 2015, 14, 387. (9) Tian, S.; Wang, J.; Li, Y.; Li, D.; Xu, L.; Hou, T. The application of in silico drug-likeness predictions in pharmaceutical research. Adv. Drug Delivery Rev. 2015, 86, 2−10. (10) Bangham, A. D.; Horne, R. W. Action of saponin on biological cell membranes. Nature 1962, 196, 952−953. (11) Qian, Q.; Nath, K. A.; Wu, Y.; Daoud, T. M.; Sethi, S. Hemolysis and acute kidney failure. Am. J. Kidney Dis. 2010, 56, 780− 784. (12) Rapido, F. The potential adverse effects of haemolysis. Blood Transfus. 2017, 15, 218−221. (13) Amin, K.; Dannenfelser, R.-M. In vitro hemolysis: guidance for the pharmaceutical scientist. J. Pharm. Sci. 2006, 95, 1173−1176. (14) Lin, F.; Wang, R. Hemolytic mechanism of dioscin proposed by molecular dynamics simulations. J. Mol. Model. 2010, 16, 107−118. (15) Chaudhary, K.; Kumar, R.; Singh, S.; Tuknait, A.; Gautam, A.; Mathur, D.; Anand, P.; Varshney, G. C.; Raghava, G. P. S. A web server and mobile app for computing hemolytic potency of peptides. Sci. Rep. 2016, 6, 22843. (16) Win, T. S.; Malik, A. A.; Prachayasittikul, V.; Wikberg, J. E. S.; Nantasenamat, C.; Shoombuatong, W. HemoPred: a web server for predicting the hemolytic activity of peptides. Future Med. Chem. 2017, 9, 275−291. (17) Cover, T.; Hart, P. Nearest neighbor pattern classification. IEEE Trans. Inf. Theory 1967, 13, 21−27. (18) Vapnik, V. N. The Nature of Statistical Learning Theory; Springer-Verlag: Berlin, Heidelberg, 1995; p 188. (19) Breiman, L. Random forests. Mach. Learn. 2001, 45, 5−32. (20) Friedman, J. H. Stochastic gradient boosting. Comput. Stat. Data Anal. 2002, 38, 367−378. (21) Rogers, D.; Hahn, M. Extended-connectivity fingerprints. J. Chem. Inf. Model. 2010, 50, 742−754. (22) Manaargadoo-Catin, M.; Ali-Cherif, A.; Pougnas, J.-L.; Perrin, C. Hemolysis by surfactants - a review. Adv. Colloid Interface Sci. 2016, 228, 1−16. (23) O’Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchison, G. R. Open babel: an open chemical toolbox. J. Cheminform. 2011, 3, 33−46. (24) Cheng, T.; Zhao, Y.; Li, X.; Lin, F.; Xu, Y.; Zhang, X.; Li, Y.; Wang, R.; Lai, L. Computation of octanol−water partition coefficients by guiding an additive model with knowledge. J. Chem. Inf. Model. 2007, 47, 2140−2148. (25) Zheng, S.; Jiang, M.; Zhao, C.; Zhu, R.; Hu, Z.; Xu, Y.; Lin, F. e-Bitter: bitterant prediction by the consensus voting from the machine-learning methods. Front. Chem. 2018, 6, 82.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jmedchem.9b00853.

■

Article

Generation of molecular descriptors; feature selection based on feature importance; evaluation metrics to assess the performance of our classification model; model training with machine learning methods and full features; model training with machine learning methods and feature subsets; Y-randomization test of the classification models; construction of the typical consensus model and selection of fast individual model; define the applicability domain of the classification model; implementation of e-Hemolytic-Classification program (PDF) Molecular formula strings, human-hemolytic-toxicitydata set (CSV)

AUTHOR INFORMATION

Corresponding Authors

*For F.L.: phone, 86-577-8577-3060; E-mail, lin1449@126. com. *For S.Z.: E-mail, [email protected]. ORCID

Fu Lin: 0000-0002-9680-1536 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We acknowledge the support of Natural Science Foundation of Zhejiang Province (grant no. LY17B030007) and the startup funding of Wenzhou Medical University. We also thank the reviewers for their helpful suggestions to further improve this work.

■

ABBREVIATIONS USED ADMET,absorption, distribution, metabolism, excretion, and toxicity; HD50,dose causing 50% of maximum hemolysis; ECFP,extended-connectivity fingerprint; ML,machine-learning; KNN,K-nearest neighbors; SVM,support vector machine; RF,random forest; GBM,gradient boosting machine; MCC,Matthews correlation coefficient; FI,feature importance; MW,molecular weight; log P,octanol−water partition coefficient; log D,distribution coefficient; CM,consensus model; FIM,fast individual model; CV,cross-validation; GA,genetic algorithm; OVF,optimal virtual fingerprint; t-SNE,t-distributed stochastic neighbor embedding; KL,Kullback−Leibler M



Article

(26) Jones, G.; Willett, P.; Glen, R. C.; Leach, A. R.; Taylor, R. Development and validation of a genetic algorithm for flexible docking. J. Mol. Biol. 1997, 267, 727−748. (27) Mitchell, M. An Introduction to Genetic Algorithms; MIT Press: Cambridge, 1998; p 209. (28) Irwin, J. J.; Shoichet, B. K. ZINC-a free database of commercially available compounds for virtual screening. J. Chem. Inf. Model. 2005, 45, 177−182. (29) Riester, D.; Wirsching, F.; Salinas, G.; Keller, M.; Gebinoga, M.; Kamphausen, S.; Merkwirth, C.; Goetz, R.; Wiesenfeldt, M.; Stürzebecher, J.; Bode, W.; Friedrich, R.; Thürk, M.; Schwienhorst, A. Thrombin inhibitors identified by computer-assisted multiparameter design. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 8597− 8602. (30) Liu, H.; Zhao, D.; Chang, J.; Yan, L.; Zhao, F.; Wu, Y.; Xu, T.; Gong, T.; Chen, L.; He, N.; Wu, Y.; Han, S.; Qu, D. Efficacy of novel antibacterial compounds targeting histidine kinase YycG protein. Appl. Microbiol. Biotechnol. 2014, 98, 6003−6013. (31) Garrison, A. T.; Abouelhassan, Y.; Kallifidas, D.; Bai, F.; Ukhanova, M.; Mai, V.; Jin, S.; Luesch, H.; Huigens, R. W. Halogenated phenazines that potently eradicate biofilms, MRSA persister cells in non-biofilm cultures, and mycobacterium tuberculosis. Angew. Chem., Int. Ed. 2015, 54, 14819−14823. (32) Zhang, S.; Song, J.; Gong, F.; Li, S.; Chang, H.; Xie, H.; Gao, H.; Tan, Y.; Ji, S. Design of an α-helical antimicrobial peptide with improved cell-selective and potent anti-biofilm activity. Sci. Rep. 2016, 6, 27394. (33) Lee, J.; Kang, D.; Choi, J.; Huang, W.; Wadman, M.; Barron, A. E.; Seo, J. Effect of side chain hydrophobicity and cationic charge on antimicrobial activity and cytotoxicity of helical peptoids. Bioorg. Med. Chem. Lett. 2018, 28, 170−173. (34) Rücker, C.; Rücker, G.; Meringer, M. Y-Randomization and its variants in QSPR/QSAR. J. Chem. Inf. Model. 2007, 47, 2345−2357. (35) Tropsha, A. Best practices for QSAR model development, validation, and exploitation. Mol. Inf. 2010, 29, 476−488.

N


In Silico Prediction of Hemolytic Toxicity on the ... - ACS Publications

Recommend Documents