The Development, Validation, and Use of Quantitative Structure Activity Relationship Models of 5-Hydroxytryptamine (2B) Receptor Ligands to Identify Novel Receptor Binders and Putative Valvulopathic Compounds among Common Drugs
Some antipsychotic drugs are known to cause valvular heart disease by activating serotonin 5-HT2B receptors. We have developed and validated binary classification QSAR models capable of predicting potential 5-HT2B binders. The classification accuracies of the models to discriminate 5-HT2B actives from the inactives were as high as 80% for the external test set. These models were used to screen in silico 59,000 compounds included in the World Drug Index and 122 compounds were predicted as actives with high confidence. Ten of them were tested in radioligand binding assays and nine were found active suggesting a success rate of 90%. All validated binders were then tested in functional assays and one compound was identified as a true 5-HT2B agonist. We suggest that the QSAR models developed in this study could be used as reliable predictors to flag drug candidates that are likely to cause valvulopathy.
Introduction<!>Dataset<!>Preprocessing of the Dataset<!>Dataset Division for Model Building and Validation<!>Computational Methods<!>Molecular Descriptors<!>DRAGON Descriptors<!>MolConnZ Descriptors<!>MOE Descriptors<!>Subgraph Descriptors (SG)<!>Balancing the Dataset Using Similarity Searching<!>Fingerprints (FP)<!>Dragon Descriptors<!>MolConnZ Descriptors (MZ)<!>k Nearest Neighbors (kNN) QSAR<!>Classification Based on Association (CBA)<!>Distance Weighted Discrimination (DWD)<!>Robustness of QSAR Models<!>Applicability Domain of kNN QSAR Models<!>Consensus Prediction<!>Virtual Screening and Compound Selection for Experimental Validation<!>Balancing the Dataset<!>k NN Classification<!>Classification Based on Association (CBA)<!>DWD Modeling<!>Comparison of Binary QSAR Approaches for Classifying 5-HT2B Actives vs. Inactives<!>Model Validation by Predicting Drugs Known to be 5-HT2B Actives and Valulopathogens<!>Model Validation by Predicting an External Set<!>The Importance of Variable Selection<!>Virtual Screening of the World Drug Index Database to Identify Putative 5-HT2B Ligands<!>Experimental Validation<!>Conclusions<!>Radioligand Binding Assays<!>Chemistry
<p>During the last decade, several drugs have been shown to cause cardiac valvulopathy in humans. The initial discovery of drug-induced valvulopathy occurred when the anorectic drug fenfluramine (approved by the FDA in 1973), one of the active ingredients of the anorectic drug combination fen-phen, was found to increase the risk of developing two potentially serious conditions, pulmonary hypertension and valvular heart disease (VHD), in individuals receiving these medications to treat obesity.1 More recently, a group at the Mayo Clinic reported VHD in patients taking the anti-Parkinson drug pergolide.2 After the initial 2002 report, other cases of VHD associated with pergolide or other dopamine agonists such as cabergoline used as anti-parkinsonian drugs were identified.3-5 In January of 2007, the New England Journal of Medicine published two large European studies that independently verified the association of VHD with pergolide and cabergoline.6;7 Finally, on March 29, 2007, the Food and Drug Administration issued a Public Health Advisory for the voluntary market withdrawal of pergolide. These stunning withdrawals of drugs from the market stressed the importance of elucidating the mechanism by which these drugs induce valvulopathy and of determining the valvulopathic risk that may be associated with new drug candidates or even existing drugs.</p><p>To date, all but two of the VHD-associated drugs are ergoline derivatives (dihydroergotamine, methysergide, pergolide and carbergoline) (see Table 1). The two non-ergoline VHD-associated drugs are fenfluramine1 and 3,4-methylenedioxymethamphetamine (MDMA, ecstasy),8;9 both of which are amphetamine analogues (see Table 1). Thus, it appears that compounds from both the ergoline and phenylisopropylamine families can produce VHD.10</p><p>There is increasing evidence that activation of serotonin 2B receptors (5-HT2B) may play a significant role in the pathogenesis of drug-induced valvulopathy.11-13 For instance, VHD-associated drugs such as fenfluramine,14 ergotamine,14 pergolide9;15 and cabergoline, and/or selected active metabolites (such as norfenfluramine and methylergonovine),14 all potently activate 5-HT2B receptors. Chemically similar medications that do not activate 5-HT2B receptors (e.g., lisuride) seemingly do not cause valvular heart disease, further implicating the 5-HT2B receptor—but not other receptors that bind ergopeptines/ergolines and phenylisoproylamines with high affinity—in the pathogenesis of heart-valve disease.13</p><p>Additionally, valvulopathy-associated drugs have been shown to induce DNA synthesis in cultured interstitial cells from human cardiac valves via 5-HT2B receptor activation.9 It has been suggested that the valvulopathy induced by 5-HT2B receptor agonists is caused by the inappropriate mitogenic stimulation of normally quiescent valve cells, resulting in an overgrowth valvulopathy.9;13 Although the precise signaling pathways underlying drug-induced valvulopathy remain elusive, 5-HT2B receptors are known to activate mitogenic pathways through the phosphorylation of Src kinase and extracellular regulated kinases (ERK), as well as through receptor tyrosine kinase transactivation,16;17 consistent with a role in regulating heart valve interstitial cell proliferation.</p><p>The discoveries that 5-HT2B receptors were (1) abundantly expressed in heart valves,18 (2) activated by fenfluramine and its metabolite, norfenfluramine,11;18 and (3) activated by other valvulopathy-inducing drugs9;11 suggested that 5-HT2B receptors were involved in the etiology of valvulopathy.11;18 Subsequently, several other 5-HT2B agonists were also found to be valvulopathogenic.9 Since 5-HT2B agonists have the potential of causing valvulopathic side-effects, it has been suggested that all pharmaceuticals should be screened for activity at 5-HT2B receptors prior to further commercial development.13;19</p><p>Similar to experimental high throughput screening (HTS), virtual screening (VS) is typically employed as a 'hit' identification tool.20 The experimental screening of all molecules against all biological targets is generally cost- and time-prohibitive. Therefore, pre-selection of compounds by VS that have a reasonable probability to act against a given biological target is highly attractive. Typically, VS approaches imply the use of structure based methodologies; nevertheless, we have repeatedly advocated for the use of ligand based cheminformatics approaches such as QSAR models in virtual screening (reviewed in a recent monograph21).</p><p>Herein, we report on the development of in silico screening tools for identifying compounds with potentially serious valvulopathic side effects. These tools can be employed as filters to flag and de-select the potentially harmful compounds at the preclinical stage of drug development, thereby potentially avoiding significant economic and human health consequences incurred at later stages of drug discovery. To achieve this goal, validated and externally predictive, binary QSAR models were generated for 5-HT2B active vs. inactive compounds as defined in 5-HT2B functional assays. Similar studies to develop QSAR models for 5-HT2B actives vs. inactives were reported recently by Chekmarev et al.22 However, in our investigations we considered a larger dataset that contained the most complete set of all known valvulopathogens reported by Huang et al,23 and we validated our predictions experimentally in binding assays.</p><p>To obtain the most statistically robust and predictive models, we have employed the combinatorial QSAR strategy24;25 implemented as part of our predictive QSAR modeling workflow (reviewed in Tropsha and Golbraikh26). All models were subjected to rigorous internal and external validation. The results confirmed the high external prediction accuracy of our computational models, which led us to conclude that these models can be used reliably to screen chemical databases to identify putative 5-HT2B actives. Screening the World Drug Index (WDI) database using these models led to the identification of 122 possible 5-HT2B actives; 10 of these computational hit compounds were experimentally tested in 5-HT2B radioligand binding assays at the NIMH Psychoactive Drug Screening Program (PDSP), UNC Chapel Hill (http://pdsp.med.unc.edu/). Experiments confirmed that 9 out of 10 compounds were true actives implying a hit rate of 90%. These results indicate the reliability of our computational models as efficient predictors of compounds' affinity towards 5-HT2B receptors. We suggest that the computational models developed in this study could be used as drug liability predictors similar to commonly used predictors27;28 of other undesired side effects such as carcinogenicity,29-31 mutagenicity,29;32;33 PGP binding,24 or hERG binding.34-37 Our models can be used to flag compounds that are expected to bind to 5-HT2B receptors but they cannot distinguish agonists from antagonists. Nevertheless, as demonstrated in this study, these putative 5-HT2B binders can be tested in functional assays for their potential to activate 5-HT2B receptors to further assess their valvulopathic potential.</p><!><p>The PDSP recently screened roughly 2,200 FDA-approved drugs and investigational, drug-like molecules against 5-HT2B receptors.23 However, this modeling study was initiated prior to the completion of the screening of the entire compound library. At the time this study began, screening against 5-HT2B receptors had been completed for 800 compounds. This set became the basis for our model development. After preprocessing of the 800-compound dataset and deleting duplicates, the final dataset consisted of a class of 146 'actives', and another class of 608 'inactives'. Detailed PDSP protocols are available online (http://pdsp.med.unc.edu/) and in Huang et al.23 All chemical structures were obtained from PubChem38 as SDF files. By the time our modeling studies were completed, functional data for the remainder of the 2,200 compounds (1,400 compounds) had become available. These 'new' data became a source for additional, independent validation sets.</p><!><p>For the purposes of this work, the data were curated as follows: First, all molecules were "washed" using the Wash Molecules tool in MOE39 (v.2007.09). Using this tool, we processed chemical structures by carrying out several standard operations including 2D depiction layout, hydrogen correction, salt and solvent removal, chirality and bond type normalization (all details are found in the MOE manual39). Second, we used ChemAxon Standardizer40 to harmonize the representation of aromatic rings. Finally, the analysis of the normalized molecular structures resulted in detection of 46 duplicate compounds (i.e., different salts or isomeric states). The functional data for duplicated compounds were found to be identical, so in each case a single example was removed. The curated subset of the original 5-HT2B dataset used in this work contains 754 unique organic compounds (146 actives and 608 inactives). All details about the dataset are available in Supporting Information.</p><!><p>All QSAR models generated in this study to classify actives vs. inactives were validated by predicting two external validation sets. Each dataset employed in QSAR studies was first randomly divided into a modeling and a validation sets. Additionally, as described above, an independent validation set became available after we completed our modeling studies. Details about this external set are available in Supporting Information, and in Huang et al.23</p><p>Another level of internal validation was achieved by comparing model performance for training and test sets. This approach is always employed as a part of our predictive QSAR modeling workflow26 to emphasize the fact that training-set-only modeling is not sufficient to obtain reliable models that are externally predictive.41 Thus, for each collection of descriptors, the modeling sets were further partitioned into multiple chemically diverse training and test sets of different sizes using the Sphere Exclusion method implemented in our laboratory.42 Only models that were highly predictive on the test sets were retained for the consensus prediction of the external validation sets. Finally, only those models that were shown to be highly predictive on both external sets were used in consensus fashion for virtual screening of external compound libraries.</p><!><p>A combinatorial QSAR approach (Combi-QSAR)24;25 was used to generate classification models for actives vs. inactives (Fig. 1). In this study, four types of descriptors were applied in combination with three types of statistical methods.</p><!><p>Four sets of molecular descriptors were considered in our modeling studies: Dragon,43 MolConnZ (MZ),44 MOE,39 and subgraph descriptors (SG)45 developed in this laboratory. Each type of descriptors was used separately with each of the classification methods in the context of our Combi-QSAR strategy.</p><!><p>The Dragon Professional version 5.4 software43 was used to calculate 2D descriptors. These included topological descriptors, constitutional descriptors, walk and path counts, connectivity indices, information indices, 2D autocorrelations, edge adjacency indices, Burden eigenvalues, topological charge indices, eigenvalue-based indices, functional group counts, atom-centered fragments and molecular properties. The initial descriptor set was reduced by eliminating the constant and near-constant variables using built-in functions within the software. The pairwise correlations for all descriptors were examined and one of the two descriptors with the correlation coefficient R2 of 0.95 or higher was excluded. The calculation procedures for these descriptors, with related literature references, are reported by Todeschini and Consonni.46 Finally, the remaining descriptors were normalized by range-scaling so that their values were distributed within the interval 0-1.</p><!><p>The MolConnZ software44 available from EduSoft affords the computation of a wide range of topological indices of molecular structure. These indices include, but are not limited to, the following descriptors: valence, path, cluster, path/cluster and chain molecular connectivity indices,47-49 kappa molecular shape indices,50;51 topological52 and electrotopological state indices,53-56 differential connectivity indices,47;57 graph's radius and diameter,58 Wiener59 and Platt60 indices, Shannon 61 and Bonchev-Trinajstić62 information indices, counts of different vertices, counts of paths and edges between different types of vertices (http://www.edusoftlc.com/molconn/manuals/400). Descriptors with zero values or zero variance were removed; the remaining descriptors were normalized by range-scaling so that their values were distributed within the interval 0-1.</p><!><p>MOE 2007.09 software39 was used to generate 2D MOE descriptors. These included physical properties, subdivided surface areas, atom and bond counts, Kier and Hall connectivity47-49 and kappa shape indices,50;51 adjacency and distance matrix descriptors,58;59;63;64 pharmacophore feature descriptors, and partial charge descriptors. 39 Descriptors with zero values or zero variance were removed; the remaining descriptors were normalized by range-scaling so that their values were distributed within the interval 0-1.</p><!><p>Frequent subgraph mining of chemical structures is a novel approach to generating fragment descriptors that was developed recently in our group.45 SG descriptors are derived from each dataset, i.e., not pre-defined which gives the advantage of finding important chemical fragments that may have not been defined a priori by other fragment descriptor generating methods. The fragments are derived based on recurring substructures found in at least a subset of molecules (defined by a support value σ) in the dataset. These recurring substructures can implicate chemical features responsible for compounds' biological activities. First, chemical structures were converted into labeled, undirected graph representations where nodes were labeled by atom types and edges corresponded to chemical bonds. Fast frequent subgraph mining (FFSM) algorithm65 was then used to find common frequent subgraphs for a given support value (σ), which is one of the variables defined by the user that determines the size of the set of subgraphs generated using FFSM. Obviously, the larger is the value of the support, the smaller is the number of subgraphs descriptors. As the support value decreases, the number of subgraphs increases dramatically. Redundant subgraphs were identified and removed leaving only the so called "closed subgraphs". A subgraph SGi is closed in a database if there exists no supergraph SGj such that SGi ⊆ SGj and σSGi = σSGj. However, subgraph SGi would not be deleted if it also occurs by itself (not as part of the SGj) in the graph database. Removing redundant subgraphs (fragments) reduces the number of subgraphs descriptors drastically and therefore makes the subsequent calculations more efficient. The frequency of individual 'closed subgraphs' in each molecule of the dataset is calculated and used as the descriptor value for each molecule. In this study, a support value of 12 % was used, and the upper size limit of the generated subgraphs was 7 atoms.</p><!><p>The dataset used for model building was imbalanced, consisting of 146 actives vs. 608 inactives. Therefore, only a subset of the larger class of inactives of approximately the same size as the actives was used in model building unless otherwise indicated. This subset was selected to include inactives that were most similar to the actives. Given the vast array of available chemical descriptors and the large number of similarity measures, it is always difficult to decide a priori which combination of descriptors/similarity metrics to use. This problem has been highlighted in several recent publications.66;67 Therefore, similarity searching studies were performed using three types of molecular descriptors: fingerprints (FP), Dragon, and MZ, and applying two similarity metrics, i.e., Euclidean distance and Tanimoto coefficient (Tc). The similarity cutoff was chosen to obtain the most balanced (with roughly equal number of compounds from each class) subset of compounds.</p><!><p>166 MACCS68 structural keys implemented in MOE 2007.09 software39 were calculated for all compounds. The similarity searching was performed using an in-house written script applying Tanimoto coefficients for similarity measures.</p><!><p>Normalized Dragon descriptors of the original dataset were employed to calculate similarities between all compounds in the dataset using Euclidean distance as similarity metric; variable similarity thresholds were used to down-sample the larger class (inactives). Although many schemes could be considered for down-sampling the larger classes, we used the similarity threshold based approach since it restricts the larger class to compounds most similar to the smaller class molecules. This approach makes it more challenging to develop statistically significant models capable of discriminating smaller class compounds from most chemically similar molecules in the larger class. Therefore, it increases the robustness of the binary QSAR models.</p><!><p>Similar procedures to those described above for Dragon descriptors were used.</p><!><p>The kNN QSAR method69 is based on the k nearest neighbors principle and the variable selection procedure. It employs the leave-one-out (LOO) cross-validation (CV) procedure and a simulated-annealing algorithm70;71 to optimize variable selection. The procedure starts with the random selection of a predefined number of descriptors from all descriptors. If the number of nearest neighbors k is higher than one, the estimated activities ŷi of compounds excluded by the LOO procedure are calculated using the following formula: (1)y^i=∑j=1kyjwij∑j=1kwij where yj is the activity of the j-th compound. Weights wij are defined as: (2)wij=(1+dij∑j′=1kdij′)−1 and dij is Euclidean distances between compound i and its j-th nearest neighbor. However, if the number of nearest neighbors k is equal to one, then the estimated activity ŷi of the compound will be equal to the activity of this one nearest neighbor.</p><p>For classification kNN, the predicted ŷi values (see expression (1)) are rounded to the closest whole numbers (which are, in fact, the class numbers), and the prediction accuracy (correct classification rate, CCRtrain) is calculated as follows: (3)CCR=0.5(N1corrN1total+N2corrN2total) where Njcorr and Njtotal are the number of correctly classified and total number of compounds of class j (j=1, 2). Then, a predefined small number of descriptors are randomly replaced by other descriptors from the original pool, and the new value of CCRtrain is obtained. If CCRtrain (new) > CCRtrain (old), the new set of descriptors is accepted. If CCRtrain (new) ≤ CCRtrain (old), the new set of descriptors is accepted with probability p = exp (CCR (new) - CCR (old))/T, or rejected with probability (1-p), where T is a simulated annealing (SA) "temperature" parameter. During this process, T is decreasing until a predefined threshold. Thus, the optimal (highest) CCRtrain is achieved. For the prediction, the final set of selected descriptors is used, and expressions (1) and (2) are applied to predict activities of compounds of the test sets. Then the activities are rounded to the closest whole numbers, and the correct classification rate for the test set is calculated using formula (3).</p><p>In the case when compounds belong to two classes (e.g., active and inactive compounds), a 2 × 2 confusion matrix can be defined, where N(1) and N(0) are the number of compounds in the data set that belongs to classes (1) and (0) respectively. TP, TN, FP, and FN are the number of true positives (actives predicted as actives), true negatives (inactives predicted as inactives), false positives (inactives predicted as actives), and false negatives (actives predicted as inactives), respectively. The following classification accuracy characteristics associated with confusion matrices are widely used in QSAR studies: sensitivity (SE=TP/N(1)), specificity (SP=TN/N(0)), and enrichment E = TP*N/[(TP+FP)*N(1)]. In this study, we have employed normalized confusion matrices. A normalized confusion matrix can be obtained from the non-normalized one by dividing the first column by N(1) and the second column by N(0). Normalized enrichment is defined in the same way as E but is calculated using a normalized confusion matrix: En = 2TP*N(0)/[TP*N(0)+FP*N(1)]. En takes values within the interval of [0, 2].25;72</p><!><p>This method integrates both classification rule mining,73;74 which aims to discover a small set of rules in the database that forms an accurate classifier, and association rule mining,75 which finds all the rules existing in the database that satisfy some minimum support and minimum confidence constraints. For association rule mining, the target of discovery is not pre-determined, while for classification rule mining there is one and only one predetermined target. The integration is done by focusing on mining a special subset of association rules, called class association rules (CARs). An efficient algorithm is also used for building a classifier based on the set of discovered CARs.</p><p>The CBA algorithm76;77 consists of two parts, a rule generator, which is based on the a priori algorithm for finding association rules, and a classifier builder. The candidate rule generator is similar to the a priori one. The difference is that CBA updates the support value in each step while the a priori algorithm only updates this value once. This allows us to compute the confidence of the ruleitem. A ruleitem is of the form: <condset, y> where condset is a set of items, y ∈ Y is a class label. The support count of the condset (called condsupCount) is the number of cases in the dataset (D) that contain the condset.</p><p>Next, a classifier is built from CARs. To produce the best classifier out of the whole set of rules would involve evaluating all the possible subsets of it on the training data and selecting the subset with the right rule sequence that gives the least number of errors. There are 2m such subsets, where m is the number of rules. It is a heuristic algorithm. Given two rules, ri and rj, ri precedes rj if (1) the confidence of ri is greater than that of rj, or (2) their confidences are the same, but the support of ri is greater than that of rj, or (3) both the confidences and the supports of ri and rj are the same, but ri is generated earlier than rj. If R is a set of generated rules (i.e. CARs) and D the training data, the basic idea of the algorithm is to choose a set of high precedence rules in R to cover D. The classifier follows this format: <r1, r2, . . ., rn, default_class>, where ri ∈ R. In classifying an unseen case, the first rule that satisfies the case will classify it. If there is no rule that applies to the case, it takes on the default class.</p><p>The descriptors used with CBA need to be discrete in nature76 as is the case with SG descriptors but not Dragon, MolConnZ or MOE. Hence, this method was only used with SG descriptors using CBA (v2.1) software.78</p><!><p>This method was initially proposed by Marron and Todd79 with the goal of improving the performance of SVM80;81 in high dimensional low sample size (HDLSS) contexts. The main idea is to improve upon the criterion used for "separation of classes" in SVM. SVM has data piling problems along the margin, because it is maximizing the minimum distance to the separating plane, and there are many data points that achieve the minimum. A natural improvement is to replace the minimum distance by a criterion that allows all the data to have an influence on the result. DWD does this by maximizing the sum of the inverse distances. This results in directions that are less adversely affected by spurious sampling artifacts. The major contribution of this new discrimination method is that it avoids the data piling problem, to give the anticipated improved generality. Like SVM, DWD is based on computationally intensive optimization; however, while SVM uses well known quadratic programming algorithms, DWD uses interior-point methods for so-called Second-Order Cone Programming (SOCP) problems.82 Detailed discussion of these issues may be found in Marron and Todd (2007),79 which is available with the supporting information at https://genome.unc.edu/pubsup/dwd/. All DWD computations were performed using the DWD software83 written in Matlab84 and kindly provided by Dr. Marron.</p><!><p>Y-randomization test is a widely used validation technique to ensure the robustness of a QSAR model.85 This test includes (i) randomly shuffling the dependent-variable vector, Y-vector of training sets (class labels in this study) and (ii) rebuilding models with the randomized activities (class labels) of the training set. All calculations are repeated several times using the original independent-variable matrix. It is expected that the resulting QSAR classification models, built with randomized activities for the training set, should generally have low CCRs for training, test, and external validation sets. It is likely that sometimes, though infrequently, high CCR values may be obtained due to a chance correlation or structural redundancy of the training set. However, if some QSAR classification models obtained in the Y-randomization test have relatively high CCR it implies that an acceptable QSAR classification model cannot be obtained for the given dataset by the particular modeling method. Y-randomization test was applied to all datasets considered in this study, and the test was repeated five times in each case.</p><!><p>Formally, a QSAR model can predict the target property for any compound for which chemical descriptors can be calculated. However, since the training set models are developed in kNN QSAR modeling by interpolating activities of the nearest neighbor compounds, a special applicability domain (i.e., similarity threshold) should be introduced to avoid making predictions for compounds that differ substantially from the training set molecules.86</p><p>The similarity was estimated using Euclidean distances in high-dimensional descriptors space. Compounds with the smallest distance between them have the highest similarity. The distribution of distances (pairwise similarities) of compounds in our training set is computed to produce an applicability domain threshold, DT, calculated as follows: (4)DT=y¯+Zσ Here, ȳ is the average Euclidean distance of the k nearest neighbors of each compound within the training set, σ is the standard deviation of these Euclidean distances, and Z is an arbitrary parameter to control the significance level. Based on previous studies, we set the default value of this parameter as 0.5, which formally places the boundary for the applicability domain at one-half of the standard deviation. Thus, if the distance of the external compound from at least one of its nearest neighbors in the training set exceeds this threshold, no prediction is made.86 In this study two types of applicability domains were employed. First, a global applicability domain that ensures some level of global similarity (using all descriptors for similarity calculations) between the predicted compounds and the compounds in the modeling set. The second is a local domain which is the applicability domain of each of the individual models using only descriptors used for the model building.</p><!><p>Our experience suggests that consensus prediction of the target property for external compounds, i.e., when the compound activity is calculated by averaging values predicted by all individual models that satisfy our acceptability criteria always provides the most stable and accurate solution87. In general, consensus prediction implies averaging the predictions for each compound by majority voting for classification QSAR models, using all models passing the validation criteria (e.g., CCRtrain ≥ 0.70 and CCRtest ≥ 0.70). In order to determine the confidence in the obtained predictions we need to define a consensus score. The consensus scores employed in this study take into account the total number of models used to predict the compound's activity, and the number of models that predicted the compound to belong to a specific class. Since we define two classes of compounds, i.e., class 1 (actives) and class 0 (inactives), some models may predict a compound to belong to class 0 and others may predict it to belong to class 1. As a result, a consensus score between 0 and 1 will be obtained for each of the predicted compounds. As an additional measure of confidence (and an additional applicability domain criterion) we only accepted those predictions that had an average predicted value (consensus score) above 0.7 (for actives) or below 0.3 (for inactives).</p><!><p>To identify putative actives, validated consensus models generated for 5-HT2B ligands were used for virtual screening of ca. 59,000 molecules within the WDI chemical library; the selection of hits was limited by the applicability domains of each models.88 122 compounds were identified as VS hits (by consensus agreement between all accepted models, see Table S1 of Supporting Information for details) and 10 structurally diverse and commercially available hits were purchased from different suppliers and tested at PDSP in 5-HT2B radioligand binding assays.</p><!><p>The original dataset of 146 actives and 608 inactives was first balanced by downsizing the class of inactives. Similarity searching between active and inactive compounds using Tc cutoff of 0.7 resulted in 195 inactives (that were similar to at least one active compound with Tc above 0.7), which were combined with the 146 actives to form the modeling set of 342 compounds. Dragon and MZ descriptors were generated for this 342-compound modeling set to be used separately with kNN. However, similarity searching using Dragon and MZ descriptors and applying Euclidean distance-based threshold resulted in a 304- (146 actives and 158 inactives) and 325-compound (146 actives and 179 inactives) modeling sets respectively. Thus, slightly different modeling sets were used depending on the type of descriptors.</p><!><p>kNN method was used with each of the following descriptor types: DRAGON, MZ, MOE, and SG descriptors. Models were built for the three datasets resulting from the down-sampling of the original dataset. First, a validation set (15-20% of the dataset) was excluded from each of the resulting datasets randomly. The compounds in the remaining modeling set (85-80% of the original dataset) were divided into multiple training and test sets (28-40 divisions). Multiple QSAR models were generated independently for all training sets and applied to the test sets. Generally, we accepted models with CCR values for both the training and test set greater than 0.70. kNN combined with subgraphs and Dragon descriptors were the two best performing methods based on validation set statistics (Table 2). kNN-subgraphs (kNN-SG) had a CCRevs = 0.80, while kNN-Dragon gave a CCRevs = 0.72.</p><p>Results of the Y-randomization test (Table 2) confirmed that kNN classification models with CCRtrain and CCRtest values ≥ 0.70 were robust. None of the models with randomized class labels of the training set compounds had CCRrand > 0.54 for any dataset.</p><!><p>The CBA method was applied to classify the dataset using SG descriptors. A dataset of 342 compounds (146 actives and 196 inactives), resulting from the downsizing process with FP and Tanimoto distances, was used. The dataset was split randomly into training (267 compounds) and validation sets (75 compounds). A total of 1371 closed frequent subgraphs were generated with FFSM (see Methods) from the training set using a support value of 12% and a maximum size limit of the fragments of 7. The training set consisting of 267 compounds (111 actives and 156 inactives) was then used to build the classifier in CBA. The classifier gave a CCRtrain of 0.79. Then the validation set consisting of 75 compounds (35 actives and 40 inactives) was used to assess the robustness of the classifier. The CCRevs was 0.65 which is not as high as the CCR value for the training set.</p><!><p>The DWD method was applied to classify the dataset using Dragon descriptors. A dataset of 304 compounds (146 actives and 158 inactives), resulting from the downsizing process with Dragon descriptors and Euclidean distances, was used. The dataset was split randomly into training (244 compounds) and validation sets (60 compounds). A total of 387 Dragon descriptors were generated for the training set. The training set consisting of 244 compounds (120 actives and 124 inactives) was then used to build the DWD model. The DWD model was able to group compounds in this dataset based on their biological classes with a CCRevs = 0.70 (TP=18, TN=24, FP=10, FN=8), setting the cutoff at "0.15". DWD was further used to rank Dragon descriptors according to their importance for discriminating the two classes of compounds (actives vs. inactives). DWD uses class label information where positive (for actives) and negative (for inactives) signs are assigned to each descriptor value to indicate its importance to the corresponding class. The top 20 highly weighted descriptors (based only on weights' values and ignoring the signs) are presented in Table S2 of Supporting Information.</p><!><p>The performance of different binary QSAR approaches employed as part of combinatorial QSAR strategy for 5-HT2B, and based on validation set statistics, is summarized in Figure 2. kNNSG, and kNN-Dragon were the best performing methods for classifying 5-HT2B actives vs. inactives based on validation set statistics (Table 2), yielding the highest CCRevs of 0.80 in case of kNN-SG. On the contrary, kNN-MZ was the worst performing method with a CCRevs of 0.57 which was very close to random. It was also interesting to see that kNN-SG performed much better than CBA-SG with CCRevs = 0.80 in the former case and 0.65 in the latter. These results confirm the importance of employing the combinatorial QSAR approach to find the most predictive QSAR method/descriptor combination for each specific dataset.</p><p>Our models also indicated that the nature of the descriptors used has a dramatic effect on the performance of the modeling methods. It was clear that MOE and MolConnZ descriptors did not perform very well in all tested cases irrespective of the applied modeling techniques. On the contrary, Dragon descriptors afforded most significant models with all methods and in all tests, for both validation and external sets.</p><!><p>Both fenfluramine and dexfenfluramine (known to be 5-HT2B actives and agonists, which were not included in our modeling sets) were predicted as 5-HT2B actives using consensus models to classify actives vs. inactives. The consensus scores using kNN-Dragon were 0.79 for both compounds. Our previous studies suggest that consensus prediction that is based on the results obtained by all validated predictive models always provides the most stable solution.87 A 5-HT2B active compound can have consensus scores in the interval [0.5-1.0]. The closer value to 1.0 the greater is the confidence in the prediction. Therefore, we can claim that both compounds were predicted as actives with statistically significant consensus scores.</p><p>These results highlight the predictive power of our validated models that could have predicted the possible dangerous side effects of these two drugs by suggesting that they may be 5-HT2B actives. This prediction would have suggested that these compounds should be tested experimentally in 5-HT2B functional assays and prevented from further development as potentially unsafe medicines. This example illustrates the potential use of models developed in this study as computational drug safety alerts.</p><!><p>An additional 16-compound set was obtained from PDSP after we finished out modeling studies. This external set was used to further assess the robustness and the predictive power of our models. All 16 compounds were 5-HT2B actives including 4 agonists and 12 antagonists.</p><p>The 16 external compounds were predicted using all consensus models built to classify actives vs. inactives. kNN-Dragon was the best performing method on this external set with a CCRex of 0.81. Predictions were made by applying local model applicability domains with Z = 0.5 (see Applicability Domain of kNN QSAR Models). It was interesting to find that kNN-Dragon had CCR ≥ 0.72 with both the validation (CCRevs = 0.72) and the external (CCRex = 0.81) sets. However, kNN-SG (the best performing method on validation sets) was not as good with the external set (CCRex = 0.65) as it was with the validation set (CCRevs = 0.80). CBA-SG gave a CCRex = 0.65, which was consistent with its performance with the validation set (CCRevs = 0.65) but less than CCRtrain (0.79). The latter results using SG descriptors with kNN and CBA might be due to the limitation that frequent subgraphs are derived from the training set compounds; therefore, it is possible that fragments that are frequent in the external set are not represented in the frequent subgraphs used for prediction. Our current applicability domain filter, which is calculated using the fragments in the training set, does not account for this possibility. It is clear that a more stringent applicability domain filter could be applied in this case, which uses the distribution of subgraphs counts between the training and test set, but this has not been implemented yet within our current method.</p><!><p>Since kNN-Dragon was the best performing method to classify actives vs. inactives based on the results for all validation sets, we thought it would be interesting to check the performance of kNN using all 387 Dragon descriptors, generated for the actives vs. inactives modeling set, without variable selection. The results of this test are shown in Table 3. Comparison of modeling results for kNN-variable-selection (CCRevs = 0.72) vs. kNN-without-variable-selection (CCRevs = 0.52) clearly indicates that variable selection is a vital part of modeling. Furthermore, the top 20 most frequent descriptors (MFD) selected by kNN models (Table S3 of Supporting Information) and top 20 highly weighted descriptors by DWD based only on weights and ignoring the sign (Table S2 of Supporting Information) were used independently with the kNN method (with no variable selection) to predict actives vs. inactives (Table 3). Models built with either the top 20 DWD-selected Dragon descriptors or MFD from Dragon-kNN and using 1-5 nearest neighbors gave CCRevs ~ 0.5 (Table 3). These results illustrated again that SA-based variable selection procedures implemented in our kNN QSAR method69 lead to models with the highest external predictive power as compared to any other approach not relying on variable selection for model optimization.</p><p>Mechanistic interpretability is frequently regarded as very important feature of QSAR models. We generally argue that only models that have been extensively validated on external datasets and identified experimentally-confirmed hits should be subjected to interpretation. Furthermore, very few classes of models, specifically, those based on (multiple) linear regression and small number of descriptors can afford a relatively straightforward interpretation. The interpretation of multi-parametric statistical models developed with non-linear optimization algorithms (as in this study) should be attempted with great care because of strong and often poorly understood interplay between descriptors. Furthermore, although we could foresee that in some cases medicinal chemists may want to modify their candidate compounds to prevent 5HT2B binding, the tools developed in this study are predominantly intended for virtual screening of libraries of drug candidates to flag and possibly eliminate compounds that are likely to bind 5HT2B receptor, not to design new compounds; and any compound designed by chemists could be passed through our models. Therefore, we only restricted the discussion in this paper to the most frequent descriptors found by all acceptable kNN models and the most highly weighted descriptors selected by DWD to stress that the process of variable selection employed as part of model optimization has indeed converged on a small number of descriptors.</p><!><p>Since our models proved to be reasonably accurate based on two external validation sets, we used the best models to mine a large external database of approved and potential drugs for putative 5-HT2B actives. An important condition that assures reliable predictions by the model is the use of AD. Therefore, two types of AD were employed in the virtual screening of compound databases. The first is a global AD that acts as a filter and ensures some level of global similarity between the predicted compounds and the compounds in the modeling set. The second is a local AD which is defined for each of the individual classification models.</p><p>The WDI database of ca. 59,000 compounds (approved or investigational drugs) was used for virtual screening (Fig. 3). This original collection had many duplicates (i.e., many salt forms for the same chemical entity). The duplicates were removed using MOE: keeping unique structures and deleting duplicates. We also removed all compounds included in our modeling and external validation sets. Dragon descriptors were generated for the remaining 46,859 unique compounds in the database; 926 compounds were excluded because Dragon was unable to calculate at least one of the descriptors generated for the modeling set. The remaining 45,933 compounds were then subjected to a global AD filter for the actives vs. inactives modeling set using a strict Z cutoff of 0.5 (which formally places the allowed pairwise distance threshold at the mean of all pairwise distance distribution for the training set plus one-half of the standard deviation). Obviously, increasing the AD would increase the number of computational hits identified by virtual screening. However, our experience suggests that such increase is typically accompanied by the decrease in prediction accuracy. Additionally, we required that the nearest neighbor in the modeling set of a compound from the virtual library be an active. The resulting 7,286 compounds were then classified into actives vs. inactives using DWD-Dragon classifier resulting in 891 actives. Next, all kNN-Dragon models with CCRtrain and CCRtest ≥ 0.70 were employed in consensus fashion to predict these 891 compounds resulting in a selection of the 500 active hits. At this point, SG descriptors were generated for these 500 molecules. CBA-SG classifier followed by kNN-SG consensus models were used as final filters for the determination of 122 compounds regarded as putative 5-HT2B actives.</p><!><p>Ten structurally diverse hits (1-10, see Table 4) were selected from the final consensus virtual screening hits for further experimental validation taking into account both their commercial availability and cost (see Table 4). To our satisfaction, nine compounds were confirmed to inhibit 5-HT2B radioligand binding, which implies a hit rate of 90 %. Ki values were in the range 0.8 – 3,127 nM, with 4 compounds having Ki values < 100 nM. The four highest affinity compounds were: 4 (Ki=33 nM, see Fig. 4 (A)), 7 (Ki=0.8 nM, see Setola et al, 20039), 9 (Ki=70 nM, see Fig. 4 (B)), and 10 (Ki=69 nM, see Fig. 4 (C)).. It should be noted that methylergometrine, though not included initially in our dataset, was known to be a valvulopathic compound and had been tested against 5-HT2B receptors in both binding (K =0.8 nM)9 and functional assays (pEC50 for 5-HT2B-Mediated calcium flux = 7.67)23. In order to determine the activity of the remaining eight 5-HT2B ligands, all compounds were tested at the PDSP in 5-HT2B functional assays. Results indicated that methylergometrine was the only compound among the 9 5-HT2B ligands that possessed strong agonist activity.</p><p>This low hit rate of 11.1% for identifying validated agonists is in fact not surprising in light of Huang et al23 major finding that potent 5-HT2B receptor agonism is a relatively rare occurrence among drugs and drug-like compounds. However, to arrive at such conclusions, Huang et al screened a composite library containing three publicly available collections of FDA-approved and investigational medications and one internally compiled library. Of the approximately 2200 compounds screened, 27 5-HT2B receptor agonists were identified; thus, the validated hit rate was 1.2%.</p><p>These results illustrate that the validated QSAR workflow, as employed in this paper, could be used as a general tool for identifying 5-HT2B ligands by the means of virtual screening of chemical libraries using rigorously built QSAR models. As we demonstrated in this study, our models identify a relatively small number of VS hits making it feasible to employ experimental tools to validate predictions in 5-HT2B binding and functional assays. Ten compounds selected from a large external library have been tested experimentally in this proof-of-concept study resulting in very high experimentally confirmed hit rate. The list of all compounds predicted to be 5-HT2B actives is available in the Supporting Information (Table S1).</p><p>To verify the diversity of the experimentally validated hits, we have compared the results of QSAR-based virtual screening with simple similarity searches. Similarity calculations were done using two different descriptor-metric combinations: (1) MACCS structural keys and Tanimoto coefficients (as a standard similarity searching approach, see Table S9 and Figure S1 in Supporting Information) and (2) Dragon descriptors and Euclidean distances (to compare directly with our best performing QSAR models of kNN-Dragon, see Table S10 and Figure S2 in Supporting Information). The nearest neighbor compounds (based on Tanimoto similarities and MACCS keys) from the active compounds in the dataset and the 10 experimentally validated VS hits are reported in Table 5. Results of similarity analyses indicated that neither technique would be able to efficiently identify the diverse hits obtained with our methods (see Supporting Information for details). Hence, our studies illustrated the power of combi-QSAR-based VS in prioritizing compounds (which are not just close analogs of the modeling set compounds) from screening libraries to achieve high success rates when experimentally validated.</p><p>We also think that agonist vs. antagonist models will be highly useful as more data about agonist compounds become available. The small number of known 5-HT2B agonists made it impossible at this stage to develop statistically significant models that could distinguish agonists from antagonists. Thus, the current study was limited to building binder vs. non-binder models. We will continue with our efforts to develop quantitative 5-HT2B agonist predictors as we accumulate more experimental data.</p><!><p>QSAR models are becoming increasingly attractive as robust computational tools for virtual screening due to both their computational efficiency and success rates [reviewed in26 as well as in a recent monograph21]. In this study, we have applied a combinatorial QSAR approach to a dataset of 800 compounds experimentally annotated as 5-HT2B receptor agonists, antagonists and inactives resulting in statistically validated and externally predictive models. Specifically, we have applied a combi-QSAR approach utilizing three different classification methods (kNN, CBA and DWD) and four different descriptor types (Dragon, MZ, MOE and SGs) to generate classification QSAR models to discriminate between 5-HT2B actives (agonists and antagonists) from inactives. Predictive models with classification accuracies as high as 0.80 for actives vs. inactives, as estimated on external validation sets, were obtained.</p><p>Classification models for actives vs. inactives were further validated by predicting an external validation set obtained after we completed the modeling studies. The high accuracy of prediction for the second external validation set proved that our models were indeed rigorous. Therefore, we posited that our studies afforded a robust computational tool to predict potential 5-HT2B activity and consequently prioritize hits for testing in functional 5-HT2B assays to predict valvulopathic side effects of drugs and drug candidates that act as 5-HT2B agonists. We suggested that this computational predictor could be used to eliminate high risk compounds at the early stages of the drug development process. To illustrate this point, we have used this predictor retrospectively to evaluate the valvulopathic potential of two drugs withdrawn from the U.S. market for this reason, i.e., fenfluramine and dextrofenfluramine. Both drugs were not included in our modeling set and both were indeed predicted with high confidence as actives for binding to 5-HT2B receptors.</p><p>Encouraged by our model validation results, we have applied these models for virtual screening of the 59,000 compounds in WDI database. Our classification strategies identified 122 potential 5-HT2B ligands. Ten structurally diverse VS hits were experimentally tested at PDSP. Nine compounds were experimentally confirmed as 5-HT2B ligands thereby demonstrating a very high success rate of 90%.</p><p>The predictor developed in this report is similar in its potential use to other predictors of drug liability such as carcinogenicity and mutagenicity that are widely used in pharmaceutical industry. For instance, the TOPKAT program available in the Discovery Studio,89 is a QSAR-based system that generates and validates accurate, rapid assessments of various types of chemical toxicity solely from a chemical's molecular structure. In contrast, our predictor is a unique specialized tool for the prediction of 5-HT2B activity and therefore prioritizing compounds for functional testing against 5-HT2B receptors to assess their valvulopathic potential. Therefore, this predictor can be used, along with other computational chemical health risk assessment tools, to evaluate compounds' safety at early stages of the drug development. It can be used as well to verify that all drugs available on the market are free from possibly fatal valvulopathic risk. This predictor will be made publicly available at the ChemBench server established in the Laboratory for Molecular Modeling (chembench.mml.unc.edu). We will also gladly apply this predictor to any compound library that may be of interest to any researcher.</p><!><p>This screen was performed by the National Institute of Mental Health Psychoactive Drug Screening Program (PDSP). Radioligands were purchased by PDSP from Perkin-Elmer or GE Healthcare. Competition binding assays were performed using transfected or stably expressing cell membrane preparations as previously described (Shapiro et al. 2003;90 Roth et al. 200291) and are available online (http://pdsp.med.unc.edu). All experimental details are available online (http://pdsp.med.unc.edu/UNC-CH%20Protocol%20Book.pdf).</p><!><p>Chemical compounds predicted as hits from the virtual screening were obtained from commercial suppliers according to their availability. All compounds were ordered to have ≥ 95% purity. Additionally, all compounds were subjected to purity assessment using LC/MS by the Center for Integrative Chemical Biology and Drug Discovery at UNC-Chapel Hill. LC/MS spectra of all compounds were acquired from an Agilent 6110 Series system with UV detector set to 220 nm. Samples were injected (5 uL) onto an Agilent Eclipse Plus 4.6 × 50 mm, 1.8 uM, C18 column at room temperature. A linear gradient from 10% to 100% B (MeOH + 0.1% Acetic Acid) in 5.0 min was followed by pumping 100% B for another 2 minutes with A being H2O + 0.1% acetic acid. The flow rate was 1.0 mL/min.</p>