A k-nearest neighbor classification of hERG K+ channel blockers

A series of 172 molecular structures that block the hERG K+ channel were used to develop a classification model where, initially, eight types of PaDEL fingerprints were used for k-nearest neighbor model development. A consensus model constructed using Extended-CDK, PubChem and Substructure count fingerprint-based models was found to be a robust predictor of hERG activity. This consensus model demonstrated sensitivity and specificity values of 0.78 and 0.61 for the internal dataset compounds and 0.63 and 0.54 for the external (PubChem) dataset compounds, respectively. This model has identified the highest number of true positives (i.e. 140) from the PubChem dataset so far, as compared to other published models, and can potentially serve as a basis for the prediction of hERG active compounds. Validating this model against FDA-withdrawn substances indicated that it may even be useful for differentiating between mechanisms underlying QT prolongation. Electronic supplementary material The online version of this article (doi:10.1007/s10822-016-9898-z) contains supplementary material, which is available to authorized users.


Introduction
The human ether-a-go-go related gene (hERG, KCNH2) encodes for a voltage dependent K ? ion channel (Kv11.1). Blocking of this channel has been associated with potential severe heart arrhythmia, and because of this, several drugs have been withdrawn from the market [1][2][3][4][5][6]. Further, the drug-induced long QT syndrome may cause avoidable sudden cardiac arrest [3,4]. With the intention of protecting clinical trial participants and patients, the International Conference of Harmonization published a guideline (S7B) recommending that ''all new drugs'' should be tested pre-clinically for hERG sensitivity and cardiac safety before submitting an application to regulatory reviews [7]. Accordingly, the early assessment of hERG-related cardiotoxicity has become a common practice in drug discovery.
Many in vitro assays exist for the pre-clinical evaluation of hERG-related cardiotoxicity [8], examples include rubidium-flux assays, radioligand binding assays, in vitro electrophysiology measurements, and fluorescence-based assays [9]. In addition, in silico models have been proposed for identifying potential hERG blockers in drug discovery processes [10,11].
Efforts to use computational methods for the prediction of hERG blocking effects have ranged from the use of simple rules based on structural and functional features, through to more complex quantitative structure-activity relationship (QSAR) models [12][13][14][15][16]. A number of QSAR models have been developed for the hERG toxicity endpoint using different machine learning algorithms, such as multiple linear regressions [17], partial least squares (PLS) [18], k-nearest neighbor algorithms (k-NN) [19], artificial neural networks [20], support vector machines (SVM) [21], random forest [22] and naive Bayesian classifications [23]. Despite these efforts there is significant scope for development of more powerful and more easily deployed predictive models.
The recent development of open source fingerprints, such as PaDEL fingerprints, which are libraries of descriptors [24], allows for ready access to tools for predicting biological endpoints. A recent report on the use of PaDEL fingerprints in conjunction with a k-NN strategy aimed at the prediction of chronic toxicity [25] prompted us to apply this approach to hERG-channel blockers, a far more focused system. It was envisaged that publicly available data on a series of hERG-channel blockers could function as a starting point for model construction, and a series of 1953 PubChem compounds could act as basis for validation.

Methodology
Description of dataset IC 50 data for 172 I kr ('rapid' delayed rectifier current) channel blockers were retrieved from the webservers OCHEM [26] and Fenichel [27]. These 172 compounds are structurally diverse and belong to different therapeutic classes. The compounds were authenticated with respect to structure and IUPAC name. After authentication, the SMILES notations for all the 172 compounds were verified using ChemSpider [28], SigmaAldrich [29] and PubChem [30]. A PubChem dataset comprised of 1953 entries was chosen for the external validation [31]. Dataset entries that were mixtures or salts were discarded, leading to a final PubChem validation set of 1795 compounds. More details about the training and test set compounds are provided in the Online Resources 1 and 2, respectively.

Descriptor calculation
The descriptor calculation was a primary requirement for the construction of the classification model. Eight types of PaDEL fingerprints were calculated for both the training and test set compounds using PaDEL software [24]. These consisted of the CDK, Extended CDK, CDK Graph, Estate, MACCS, PubChem, Sub-structure and Sub-structure count fingerprints. Each of the eight types of fingerprints was then used, separately, to develop a classification model.

Class assignment
The training set compounds were split into one of the two classes (active and inactive) using an IC 50 threshold value of 5 lM. The PubChem dataset derived test set compounds were similarly classified, i.e. as either active or inactive, here using a % inhibition threshold of 20 %. A summary of the numbers of the compounds and their classes is provided in Table 1.

Software and modules
The Matlab module ''classification_toolbox'' [32] was employed for the development of the k-NN classification model. The Matlab module is freely available at [33].

Classification model development
The k-nearest neighbor (k-NN) classification method employed used cross validation (CV) to identify optimal k values [34,35]. A series of k values (from 1 to 10) were assigned to construct the model, and by determining the lowest class error, optimal k values were identified.
A five-step cross validation was implemented by first dividing the training set into five equal groups, four of which were used for model construction and the remaining for validation. This procedure was repeated so that each of the five groups was used for validating the models constructed using the remaining four. After cross validation, the models were subjected to external validation using the 1795 PubChem compounds. The performance of each classification model was assessed by means of statistical parameters, such as non-error rate (NER), sensitivity, specificity, precision and error rate [36]. The models were then analysed and compared on the basis of these statistical parameters.

Results and discussion
Construction of eight k-NN classification models The k-nearest neighbor (k-NN) classification method was employed to construct classification models using each of the eight PaDEL fingerprints. Employing the k-NN algorithm requires that the optimal value of k is determined [34]. There are several ways to determine the k value, e.g. through application of a risk function or empirical rules, or through cross validation. Here, cross validation was used to determine the optimal k value. A series of eight k-NN classification models was constructed using each of the PaDEL fingerprints, and compared with respect to a series of statistical parameters, Table 2.
CDK fingerprints are one-dimensional 1024 bit long arrays that are arranged based upon the occurrence of particular structural elements. The Extended CDK fingerprints are extended versions of CDK fingerprints that include ring features. Graph fingerprints are specialized versions of the CDK fingerprints that exclude bond orders. Estate fingerprints represent the influence of substituent electronic effects in a given compound. PubChem  [37,38]. The Sub-structure fingerprints represent 307 SMARTS patterns for different functional groups, whereas the count of these SMARTS patterns is referred to as the Sub-structure count fingerprint [37]. The sensitivity expresses the prediction accuracy of hERG-active compounds, whereas specificity reflects the prediction accuracy for hERG-inactive compounds. The models performed similarly in terms of the statistical parameters examined. Thus, to further improve the predictive power of these models we developed a series of consensus models. Several methods have been reported for consensus model development [39]. For classification models, the majority principle [40] is commonly employed and we have used this strategy to develop consensus models based upon three, five and seven different fingerprint-based models. As it is more important to identify hERG-active compounds than hERG-inactive compounds, the eight models (from Table 2) were examined with respect to their sensitivity in the external prediction. The Estate-fingerprint-based model exhibited relatively poor sensitivity (0.49) and was discarded from the consensus model building procedure to provide an odd number (seven) of fingerprints. Six consensus models were built using different combinations of the seven remaining fingerprint-based models, Table 3.
Although consensus model 1 shows better overall accuracy of prediction (Q), consensus model 2 shows higher sensitivity for test set prediction, and was thus chosen for further studies.

Individual contribution of each model
With consensus model 2 in hand, we then examined how individual training set compounds were handled by the consensus model as well as the individual models, i.e. Extended CDK, PubChem and Substructure count fingerprint based, Fig. 1.
The consensus model correctly predicted 121 of the 172 training set compounds. 69 of these 121 compounds were predicted correctly by all three individual models, while the remaining 52 compounds were correctly predicted by any two of the three models. Conversely, the consensus  Sensitivity Â Specificity p Fig. 1 Venn diagram representing the number of training set compounds correctly predicted by all three models (yellow), by any two models (magenta), by only one model (blue) and by none of the models (green). The shaded area represents compounds correctly predicted by the consensus model model incorrectly predicted 51 training set compounds. Of these 51, 25 compounds were predicted correctly by any one of the three models, whereas the remaining 26 compounds were incorrectly predicted by all three models.
In the case of the Extended fingerprint based model, 113 of 172 compounds were correctly predicted, 65 of which were hERG actives. The PubChem fingerprint based model predicted 105 compounds correctly from the training set. Among the 105 correctly predicted compounds, 66 were from class 1 and 39 from class 2. The Substructure count fingerprint based model predicted 118 training set compounds correctly. These 118 compounds were comprised of 67 compounds from class 1 and 51 compounds from class 2.
Compounds for which activities were not correctly predicted by our models are of interest as awareness of factors contributing to the incorrect prediction of compounds can help in the refinement of models. In this case, the IC 50 value-based endpoints are derived from a range of studies so impact of inter-laboratory variation in the reported IC 50 data on model performance cannot be excluded.

Comparison of our model with other models
External validation provides an assessment of the QSAR model's performance, and to compare models it is necessary that the external validations are performed on the same dataset. The PubChem dataset is comprised of 221 hERGactives and 1574 hERG-inactives. Sensitivity and specificity are generally used to assess classification performance in imbalanced binary class studies [41]. G-mean, which is a geometric mean of sensitivity and specificity, was also used to measure the performance of the classification method in predicting actives and inactives. In studies aimed at the effective detection of only one class, as in our case where the prediction of hERG-actives is a priority, sensitivity and F-measures are often adopted [41]. Accordingly, we have compared our model with previously published models that were externally validated with the PubChem dataset [18,[42][43][44], with respect to sensitivity, specificity, G-mean and F-measure, Table 4.
As presented in Table 4, three of the four previously described models demonstrate lower overall sensitivities than our model, though it should be pointed out that IC 50 thresholds used in the various studies varied between 5 and 40 lM. From a drug development perspective, it may be argued that it is of more interest to identify the potent hERG blockers (class 1) than hERG inactive compounds (class 2). Comparison on this point reveals that our model demonstrates better performance in predicting the hERG active compounds (True positives = 140, Sensitivity = 0.63) than the other models except that of Su et al. [18] in their model presented 2010. There, 163 hERG actives from the Pub-Chem dataset were used for the external validation, whereas in our study a somewhat more comprehensive external validation was performed using 221 hERG actives. From a practical perspective, ease of use is an issue of importance and an advantage of our model is that PaDEL fingerprints are fast and easy to calculate and do not involve complicated descriptor selection procedures. This is in contrast with all the other models presented in Table 4 that all employed 3D and 4D descriptors that require geometry optimization, a task necessitating significant computational resources. In addition, the application of different descriptor selection procedures makes these tasks more cumbersome. Therefore, in comparison to the other models, our model has the advantage of being fast, simple and relatively efficient in predicting hERG toxic compounds.
To further assess the potential of our consensus model, we turned our attention to the series of 47 substances withdrawn from use on account of QT-prolongation, which can be hERG-derived, as present in the WITHDRAWN database [45] (database last updated December 2015). Our training set had included 32 of these 47 drugs (shown in bold in Online Resource 1) of which our model had correctly predicted the IC 50 -based classes of 22. We interrogated the remaining 15 withdrawn substances (see Online Resource 3) using our model, which correctly predicted the IC 50 -based classes of 11 (73 %, see Online Resource 4). It is important to note that our model is solely based upon in vitro data (hERG IC 50 ), while the basis for withdrawal, QT prolongation, is in vivo data-derived. The interpretation of the QT prolongation endpoint is itself a major challenge as mechanisms other than hERG activity can also underlie QT prolongation [4,46,47]. This is reflected in the fact that substances were correctly classified as class 1 or class 2, five and six substances respectively, based on their hERG IC 50 . This observation suggests that the model may even be useful for differentiating between mechanisms underlying QT prolongation.
A general reflection upon examining the hERG active compounds predicted by our model was the prevalence of aromatic and basic functionalities in these compounds (for example, see Online Resource 2). These features have previously been identified as essential components in a pharmacophore for central nervous system activity [48,49] and we believe should be considered in future model development. Moreover, this may be considered indicative of a common evolutionary origin for the hERG voltage dependent K ? ion channel and CNS receptors [50,51].

Conclusion
In conclusion, PaDEL fingerprint-based k-NN classification models presented here show potential as tools for the prediction of the hERG toxicity endpoint, an important issue in modern drug development. In particular, the consensus model developed using the Extended CDK, Pub-Chem and Sub-structure count fingerprint-based models performed comparably with models employing more complicated descriptors in the validation with external datasets. Moreover, the model presented here, in terms of the prediction of hERG toxicity, compares most favorably with these previously published models. Moreover, validating this model against FDA-withdrawn substances indicates that the model may be useful for differentiating between hERG-derived QT prolongation and other QT prolongation mechanisms. Accordingly, we believe that this model may provide a basis for improved drug design.