Characterizing the Chemical Space of ERK2 Kinase Inhibitors Using

This is an open access article published under an ACS AuthorChoice License, which permits copying and redistribution of the article or any adaptations for non-commercial purposes.

Article pubs.acs.org/jcim

Characterizing the Chemical Space of ERK2 Kinase Inhibitors Using Descriptors Computed from Molecular Dynamics Trajectories Jeremy Ash and Denis Fourches* Department of Chemistry, Bioinformatics Research Center, North Carolina State University, 322 Ricks Hall, 1 Lampe Drive, Raleigh, North Carolina 27695, United States S Supporting Information *

ABSTRACT: Quantitative Structure−Activity Relationship (QSAR) models typically rely on 2D and 3D molecular descriptors to characterize chemicals and forecast their experimental activities. Previously, we showed that even the most reliable 2D QSAR models and structure-based 3D molecular docking techniques were not capable of accurately ranking a set of known inhibitors for the ERK2 kinase, a key player in various types of cancer. Herein, we calculated and analyzed a series of chemical descriptors computed from the molecular dynamics (MD) trajectories of ERK2-ligand complexes. First, the docking of 87 ERK2 ligands with known binding affinities was accomplished using Schrodinger’s Glide software; then, solvent-explicit MD simulations (20 ns, NPT, 300 K, TIP3P, 1 fs) were performed using the GPU-accelerated Desmond program. Second, we calculated a series of MD descriptors based on the distributions of 3D descriptors computed for representative samples of the ligand’s conformations over the MD simulations. Third, we analyzed the data set of 87 inhibitors in the MD chemical descriptor space. We showed that MD descriptors (i) had little correlation with conventionally used 2D/3D descriptors, (ii) were able to distinguish the most active ERK2 inhibitors from the moderate/weak actives and inactives, and (iii) provided key and complementary information about the unique characteristics of active ligands. This study represents the largest attempt to utilize MD-extracted chemical descriptors to characterize and model a series of bioactive molecules. MD descriptors could enable the next generation of hyperpredictive MDQSAR models for computer-aided lead optimization and analogue prioritization. community challenge,8 we showed that even the most reliable quantitative structure−activity relationships (QSAR) models and sophisticated molecular docking techniques were actually not capable of accurately ranking a set of ERK2 kinase inhibitors.9 Interestingly, they were not even able to distinguish the most active ERK2 inhibitors from weakly active and inactive compounds. These negative results were disappointing, especially when considering that all participants in that CSAR community challenge failed to afford better results,8 thus making the ERK2 data set a very challenging but real case study for benchmarking novel predictive cheminformatics technologies. Computational models have become everyday tools for researchers in drug discovery. QSAR models typically rely on structural fingerprints and 2D/3D molecular descriptors to characterize chemicals and forecast their experimental activities. These descriptors are generally fast to compute, making QSAR

1. INTRODUCTION Kinase inhibitors are highly bioactive molecules of great interest for combating various types of cancer.1,2 With tens of compounds already on the market or currently in clinical trials,3 the research community is on the path to developing a large collection of potent kinase-mediated drugs and chemical probes. In particular, the identification of selective inhibitors with a finetuned kinase polypharmacology is critical. However, this research heavily relies on the expensive and time-consuming synthesis of large matching molecular series (MMS)4,5 and the exploration of structure−activity relationships within these local chemical spaces of structural analogues. In that context, the ERK2 (or MAP) kinase is considered as a key player in various cancers.6,7 Several ERK2 inhibitors have been discovered,3 but more potent and selective ERK2 inhibitors are needed. To do so, both ligand-based and structure-based molecular modeling techniques can be utilized to accelerate the discovery process and prioritize new active analogues to be synthesized and tested. Previously, in the context of the CSAR © 2017 American Chemical Society

Received: January 25, 2017 Published: May 4, 2017 1286

DOI: 10.1021/acs.jcim.7b00048 J. Chem. Inf. Model. 2017, 57, 1286−1299

Article

Journal of Chemical Information and Modeling

study. First, the structure was preprocessed using the Protein Preparation wizard14 in the Schrodinger Suite v2015-4. Bond orders were assigned to untemplated residues. All hydrogens were removed from the original structure and then added back. There was no missing side chain, and all water molecules were deleted. The PROPKA15 program was used to predict the protonation states for protein residues at pH 7. The hydrogen bonding network was optimized by the minimization of sampled hydrogens. Finally, a restrained minimization was performed using the OPLS316 force field. In this study, all small-molecule structures were docked into the active site of the target protein using Glide and two scoring functions: the standard docking precision (Glide SP) and extra docking precision (Glide XP).17,18 The ERK2 set of ligands was further preprocessed according to the following protocol. Explicit hydrogen atoms were added, and ionizable compounds were converted to their most probable charged forms at pH 7.0 ± 2.0 using the LigPrep19 software. Up to eight tautomeric forms were generated for each ligand. Chiralities were retained if specified. If undefined, at most 32 stereoisomers were generated per ligand. The binding region was defined by a 10 Å × 10 Å × 10 Å grid box. A scaling factor of 0.8 was applied to the van der Waals radii. The OPLS3 force field was used for grid generation. Default settings were used for all remaining parameters. Once Glide docking was performed, the top-five poses with the lowest Emodel score underwent postdocking minimization. The top ranked pose was selected for each ligand, and then ligands were sorted according to the docking scores (termed SP and XP scores in this manuscript depending on the precision mode used) of the top member. The poses used for downstream analysis were those predicted using the docking (Glide SP or XP) that resulted in the largest Spearman correlation with pKi values. Following molecular docking, independent 20 ns MD simulations were performed for each ERK2 ligand-protein complex using the GPU-accelerated Desmond software.20 Simulations were conducted with a TIP3P explicit solvent model. We chose a 1 ps recording interval, resulting in a MD trajectory with 20,000 unique conformations per ligand. The NPT ensemble was employed with a temperature fixed at 300 K and pressure at 1.01 bar. The integration time step was set at 2 fs. Default settings were used for all other parameters. 2.3. Computation of Molecular Descriptors. The following types of descriptors were computed using the KNIME software:21 1D-MACCS fingerprints,22 2D-RDKit descriptors,23 3D-D Moments descriptors,24 and 3D-WHIM descriptors.25 For the 3D-WHIM descriptors, all available weightings schemes were used: unit weights, atomic masses, van der Waals (VdW) volumes, Mulliken atomic electronegativities, and atomic polarizabilities. MD descriptors were computed as follows: 400 conformations were sampled at regular intervals (50 ps) for each ligand’s MD trajectory. 3D-D Moments and 3D-WHIM descriptors were computed for each ligand conformation. For each ligand, the mean and standard deviation of each 3D descriptor distributions were calculated as follows

models suitable for virtual screening of large sets of chemicals. Although very high reliability is not required per se when screening a set of 10 million of compounds (as compared to hit enrichment), it is an absolute requirement when prioritizing the top-10 analogues to be synthesized among a set of 200 possible ones. The latter task necessitates a cheminformatics technique capable of evaluating the pKi or IC50 value with very high reliability. However, as shown by the aforementioned community-wide benchmarking on the ERK2 data set,8 the prediction accuracy of the current models is far from being excellent. Thus, there is a need for next-generation, hyperpredictive QSAR models capable of facilitating the reliable screening of small focused chemical libraries of analogues in order to identify and prioritize the most promising chemicals. Interestingly, the recent development of GPU computing has dramatically decreased the resources needed to run molecular dynamics (MD) simulations of protein−ligand complexes (e.g., up to one microsecond of simulation per day on a GPU workstation10). We hypothesize that chemical descriptors computed from MD trajectories could be utilized to better characterize dynamic noncovalent protein−ligand interactions and thus build targetspecific QSAR models with enhanced prediction performances. This research hypothesis is not new and has been explored by Hopfinger and co-workers11 at a lower scale (see Discussion). Herein, we present the calculation, characterization, visualization, and benchmarking of the MD chemical space of a set of ERK2 inhibitors. To do so, we computed several pools of chemical descriptors at different levels of molecular representation. Traditional MACCS keys and 2D RDKit descriptors were considered first. Then we docked 87 ERK2 ligands with known binding affinities using Schrodinger’s Glide software. Solventexplicit MD simulations were performed for each ERK2-ligand complex. MD descriptors were computed for representative samples of the ligands’ conformations over the MD simulations. In this paper, we show that MD descriptors can provide key information about the unique dynamic characteristics of ERK2 active ligands and help discriminate the most active ones. This study represents the largest attempt of utilizing MD-extracted descriptors to characterize and model a series of bioactive molecules. This technology could enable the next generation of hyperpredictive MD-QSAR models for computer-aided lead optimization and analogue prioritization.

2. MATERIALS AND METHODS 2.1. ERK2 Data Set. In this study, we considered the ERK2 data set we previously analyzed.9 At that time, we mined ChEMBL1212 for all known ligands associated with ERK2 in order to construct the modeling set. Only 91 Ki values were retrieved for ERK2 (CHEMBL4040). Chemical curation13 encompassing the deletion of stereoisomers and compounds with uncertain and approximate Ki values was conducted. In the end, pKi values for the 48 remaining unique compounds ranged from 4.6 to 8.7. Meanwhile, the external CSAR set comprised 39 ligands with experimental pKi values provided by the CSAR organizers after the end of the community benchmark. Their activities ranged from 4.8 to 9.0. Overall, ligands with pKi above 7.5 were considered active, while those below 7.5 pKi were considered inactive. The full data set is available in Supporting Information. 2.2. Protein Preparation, Glide Docking, and Molecular Dynamics. The X-ray crystal structure of ERK2 (PDB: 3I60) was used as input for the molecular docking we conducted in this

n

xi̅ =

∑ j = 1 xij n n

si = 1287

∑ j = 1 (xij − x ̅ )2 n−1 DOI: 10.1021/acs.jcim.7b00048 J. Chem. Inf. Model. 2017, 57, 1286−1299

Article


measured using Cohen’s D value.30 The classical definitions of moderate effect size (absolute Cohen’s D greater than 0.5) and large effect size (absolute Cohen’s D greater than 0.8) were used.30 2.5. Cluster Analysis and Identification of Activity Cliffs. For activity cliff identification, we defined an activity cliff as a pair of ligands that have an absolute difference in pKi (|ΔpKi|) greater than 1.5 and a Tanimoto distance less than 0.2 in the MACCS fingerprint descriptor space or a Euclidean distance less than 6 in all other descriptor spaces. Hierarchical clustering was performed using fastcluster according to an average linkage.31 All descriptors were z-score normalized prior to the analysis so that they were on the same scale. Circular dendrograms were constructed using the ggtree package.32 The cophenetic correlation is a measure of the agreement between hierarchical clustering dendrograms.33,34 The cophenetic distance between two observations (i.e., compounds here) in a dendrogram is the height of the node at which they first join. To compute the cophenetic correlation between two dendrograms, all (N(N−1)/2) pairwise cophenetic distances were computed for each dendrogram. The cophenetic correlation is the Pearson correlation of cophenetic distances. The cophenetic correlation was computed using the dendextend package.35 The cophenetic correlation coefficient (CCC) of a dendrogram is the Pearson’s correlation of all pairwise cophenetic distances for a dendrogram versus the pairwise distances in the corresponding descriptor space. This is a measure of whether the distance between observations in a dendrogram is a good representation of the “real” distance between compounds in the true descriptor space.33,34 Sokal and Rohlf33 considered a CCC > 0.8 to indicate that there are only small distortions going from “real” distances to cophenetic distances. To assess how well ERK2 active and inactive ligands clustered together, external validation was used to select the optimal cut in the dendrogram or, equivalently, the optimal number of clusters. The Adjusted Rand Index36 (ARI) was computed for every possible number of clusters that could be generated by cutting at a height in the cluster dendrogram. The Rand Index37 (RI) is a commonly used measure for external validation of clustering

where xi̅ and si are the mean and standard deviation for the ith 3D descriptor, xij is the value of the ith 3D descriptor for the ligand conformation at time step j, and n is 400. Sampling 400 conformations for each ligand gave us an approximately normal distribution for most 3D descriptors, making the mean and standard deviation suitable summary measures for the distributions. The initial numbers of descriptors are given in Table 1. For each descriptor set, all descriptors whose variance was in the Table 1. Types and Counts of Descriptors Used in This Study descriptor set

total

post low variance filter

post correlation filter

|r| > 0.4 with pKi

active and inactive sig diff

MACCS 2D-RDKit 3D-Whim, DMoments MD

166 117 87

121 87 54

97 65 20

5 3 5

18 27 6

174

108

44

2

12

lower quartile of descriptor variances were considered low variance descriptors and thus discarded. Table 1 shows the number of descriptors that remained. All data analysis was performed using R 3.3.1.26 2.4. Hypothesis Testing and Analysis of Correlation Structure. The caret package27,28 was used to identify and filter highly correlated pairs of descriptors. If the absolute correlation (|r|) between any pair of descriptors exceeded 0.9, the descriptor with the largest mean |r| was removed. The remaining numbers of descriptors are shown in Table 1. This filtering avoided overcounting the number of significant descriptors and increased the statistical power by decreasing the number of multiple comparisons we had to correct for (see below). Paired t tests were used to test for the significant differences between the mean of active and inactive ligand descriptors. To correct for multiple comparisons, p-values were adjusted to qvalues using a Benjamini-Hocheberg (BH)29 control of the false discovery rate (FDR). A significance level of 0.05 was used, ensuring that the FDR was no larger than 0.05. Effect size was

Figure 1. Binding mode of ERK2 inhibitor CHEMBL570366. Hydrogen bond donor regions (purple) and hydrogen bonding accepting regions (green) of the binding pocket are shown. 1288


Article


Figure 2. Distribution of pKi values and Glide scores for all Model and CSAR set compounds. Stacked histograms for the (A) pKi, (B) SP GlideScore, and (C) XP GlideScore distributions. Boxplots for (D) pKi, (E) SP GlideScore, and (F) XP GlideScore distributions for the Model and CSAR set compounds.

quality and was selected because it equally rewards true positives (compounds of the same activity class are clustered together) and true negatives (compounds that are of different class are in separate clusters) RI =

Purity =

m

∑ Pi i=1

where m is the number of clusters, and n is the total number of observations. Therefore, the purity is a measure of how homogeneous the class labels are in each cluster. 2.6. Principal Components Analysis. A Yeo-Johnson transformation of the data to correct the heavy skewness in some variables was performed using the caret package prior to PCA. PCA plots were constructed using the ggbiplot39 and plotly40 packages.

TP + TN TP + FP + FN + TN

where TP, TN, FP, and FN stand for true positive, true negative, false positive, and false negative, respectively. The Adjusted Rand Index was adjusted for the random clustering of elements. An ARI equal to 0 indicates that ligands of similar class are clustering together as much as would be expected by chance, while an ARI equal to 1 indicates a perfect clustering of active and inactive ligands. Purity38 was also computed for the optimal clustering to assess the homogeneity of class labels within each cluster. We computed the purity for each cluster according to Pi =

ni n

3. RESULTS 3.1. Glide Docking of ERK2 Ligands. We used both Glide SP and XP methods to dock all 87 ligands (both Model and CSAR sets) into the binding site of ERK2. Figure 1 shows the highly active Model51 compound (CHEMBL570366, pKi = 8.4) in its docked conformation predicted by Glide SP. The extensive hydrogen bonding formed between this nanomolar inhibitor and the ERK2 pocket resulted in very low docking scores (Gscore,SP = −13.17, Gscore,XP = −12.34 kcal/mol). However, when Glide docking scores were used to rank the whole set of 87 ligands, we found moderate Spearman rank correlations (ρ) with experimental pKi values (Figure S1, Tables S1 and S2). Interestingly, Glide docking afforded reasonable

1 Maxj(nij) ni

where Pi is the purity for cluster i, nji is the number of observations in cluster i that have class j, and ni is the total number of observations in cluster j. The purity of the cluster is the fraction of the observations in cluster j that have the majority class. The overall purity is a weighted sum of each cluster’s purity 1289


Article


Figure 3. Heatmaps of (A) 116 MACCS, (B) 87 2D, (C) 54 3D, and (D) 108 MD descriptors computed for each ligand. Hierarchical clustering was performed on the descriptors (columns), whereas the ligands (rows) have been ordered from most to least active. Descriptors have been 0 to 1 normalized. Cells are colored by descriptor values (0: blue, 1: yellow).

prediction accuracy for the 48 Model Set compounds (SP: ρ = 0.74, XP: ρ = 0.65). However, the prediction performance on the CSAR Set was considerably worse (SP: ρ = 0.43, XP: ρ = 0.51). These new results obtained with the most recent version of Glide confirmed our previous conclusions on the high difficulty to (i) accurately rank these ERK2 inhibitors using molecular docking, (ii) discriminate the most active ERK2 inhibitors, and (iii) discard the most inactive compounds.9 As illustrated in Figure S1, many CSAR compounds have their binding affinities dramatically overestimated by both SP and XP scoring functions. For example, CSAR6 obtained SP and XP docking scores as good as −12.39 and −12.73 kcal/mol, respectively, meaning this compound would be predicted to be a strong active. These docking scores even put that compound at predicted ranks in the CSAR set of sixth and eighth, respectively. However, the experimental pKi of CSAR6 is only 6.1 (inactive if one takes 1 μM as the threshold), which corresponds to rank 28th.

Figure 2 shows the distribution of the computed docking scores and the experimental pKi values for the CSAR and Model Set ligands. The distribution of pKi values for the CSAR set is slightly shifted up relative to the Model Set (i.e., the CSAR set compounds are, on average, more active), though the difference in the mean pKi values for CSAR and Model Set compounds is not statistically significant. These observations on the pKi distributions translate to the docking score distributions. There is a shift downward in the distribution of CSAR docking scores relative to the Model Set, though again the differences in the means were not statistically significant. One should note there are no large outliers in the pKi distribution, so the presence of outliers with extremely low or high binding ERK2 affinity is not the reason for the low ranking reliability obtained by structurebased docking. 3.2. Analysis of Descriptor Profiles. Next, we computed the classical 166 MACCS structural keys as well as two additional sets of 2D and 3D molecular descriptors for each of the 87 ERK2 1290


Article


Figure 4. Correlation between descriptors and their relationship with ligand activity. (A) Boxplots of Pearson’s correlation between descriptors in different descriptor sets. (B) Stacked histograms of the correlation between each MACCS, 2D, 3D, and MD descriptors and ligand pKi. (C) Manhattan plot showing the BH corrected p-values for tests of association between ligand activity and each descriptor. The blue line shows significance threshold for a FDR of 0.05.

exceeded 0.9, the descriptor with the largest mean absolute correlation was removed. The remaining numbers of descriptors in each descriptor set are given in Table 1. Interestingly, several pairs of 3D descriptors with |r| > 0.9 afforded |r| < 0.9 when those same descriptors were used to construct the corresponding MD descriptors. For example, the “Ftf.Mean” 3D D-Moments descriptor has |r| > 0.9 with three other 3D descriptors (“Ctd.Mean”, “Fct.Mean”, “Atomic.Masses.WT”). However, the standard deviation of the corresponding MD descriptor (“Ftf.Mean.sd”) afforded a max |r| of 0.56 with all other descriptors. In fact, this scenario occurred for 11 3D descriptors. In each case, a descriptor afforded |r| > 0.9 with another 3D descriptor, whereas its corresponding MD descriptor (mean or SD) was characterized with |r| < 0.9 for all other MD descriptors. This result means that, at the 3D level, many pairs of descriptors encoded the same information about ligands’ molecular structures, whereas MD descriptors were more distinct and unique. Figure 4A shows the Pearson pairwise correlation between all descriptors after redundant descriptors were removed within each descriptor set. Only a few descriptors between descriptor sets are highly correlated postfiltering, suggesting that each descriptor set contains distinct information. In particular, these results again emphasized that MD descriptors contained new and unique information about the ERK2 ligands, information that is not captured by the three other descriptor sets. Postfiltering, only three pairs of descriptors from different descriptor sets obtained |r| > 0.9. The first was MACCS fingerprint bit91 (presence/absence of the substructure “QHAAACH2A”, Q being an heteroatom and A being any

ligands. Note that the 3D descriptors are computed for the toppose of each ligand obtained from molecular docking. In parallel, we performed MD simulations for each ERK2-ligand complex (see Methods). For each ligand, we computed MD descriptors as the mean and standard deviation of the distributions of a ligand’s 3D descriptors over the course of its MD trajectory. These MD descriptors describe the distribution of the 3D conformations sampled over each ligand’s MD simulation. Figure 3 shows the distribution profiles for each descriptor set. ERK2 inhibitors with a pKi above 7.5 pKi were classified as actives, whereas compounds with a pKi below 7.5 were classified as inactives. For the 3D and MD descriptors, there is a subset of descriptors where there is a clear difference in the profiles of active and inactive compounds. However, for MACCS and 2D descriptors, the difference in active and inactive profiles is less apparent. To better quantify these observations, paired t tests were used to test for significant association between ligand activity and interval valued (continuous) descriptors. For binary valued descriptors, logistic regression was used. To correct for multiple testing, the Benjamini-Hochberg procedure for false discovery rate (FDR) control was employed with a significance level of 0.05 (see Methods). Prior to testing, the correlation structure of the descriptors was examined, and descriptors with redundant information were eliminated. To examine the correlation structure of the sets of descriptors, the Pearson’s correlation was computed between all pairs of descriptors within the same descriptor set (for two binary variables, this is equivalent to the phi correlation coefficient, and for one interval and one binary variable, this is the point-biserial correlation coefficient). Again, if the |r| value between any pair of descriptors 1291


Article

Journal of Chemical Information and Modeling atom) and SMR-VSA4, an RDKit descriptor related to the polarizability of the molecule. Next was MACCS fingerprint bit110 (presence/absence of a cyanate functional group) and MQN12 (heavy atom count) in the RDKit descriptor set. The final pair of descriptors was the atomic mass weighted 3D Whim descriptor (WV) and the mean of the same 3D descriptor’s distribution over ligand’s MD conformations (i.e., the WV MD descriptor). Figure 4B shows the correlation between the descriptors in each descriptor set with the ligands’ activity toward ERK2. No descriptor afforded a strong correlation with pKi. The max |r| between pKi and the continuous descriptors was a 0.57 correlation with a 3D Whim descriptor (atomic masses weighted WV). The max of the |r| between pKi and the binary fingerprints was found to be 0.59 with MACCS fingerprint bit113 (presence/ absence of “Onot%A%A”, with %A being an aromatic element and not%A being a nonaromatic element). It is important to note that a high Pearson’s correlation between a binary fingerprint and pKi values indicates a large difference in the mean pKi values for compounds with and without that particular substructure. Meanwhile, there are several descriptors within each descriptor set that showed moderate correlation with pKi. The number of descriptors with |r| > 0.4 is given in Table 1. Figure 4C shows the BH corrected p-values for the paired t tests and logistic regressions along with the significance threshold for a 0.05 FDR (see Methods). The number of significant descriptors in each descriptor set is given in Table 1 (also see Table S3 for all descriptor names). Figure S4 shows that all of these descriptors have a moderate effect size (absolute Cohen’s D greater than 0.5) and many have a large effect size (absolute Cohen’s D greater than 0.8). For instance, the MD descriptors “Atomic.Masses.WV.mean” and “Atomic.Masses.WV.sd” have a Cohen’s D of 1.10 and 0.93, respectively. This means these two descriptors can discriminate ERK2 actives from ERK2 inactives reasonably well. Several MD descriptors were indeed found to be significantly associated with ERK2 binding affinity. Interestingly, all significant MD descriptors defined as standard deviations of 3D descriptor distributions showed a larger standard deviation for active ligands (Figure S4). Table S3 identifies five MD descriptors (“Atomic.Electronegativities.wlambda3.mean”, “VdW.Volumes.WV.sd”, “Atomic.Electronegativities.Wlambda2.mean”, “Atomic.Masses.WA.sd”, and “Atomic.Masses.Wlambda1.sd”) that showed a significant association with ERK2 activity, while their corresponding 3D WHIM did not. Therefore, if one were to examine only the top-ranked docked conformation of the ligands, one would infer that these 3D descriptors were not important determinants of ligands’ activity. Meanwhile, MD simulations revealed that their corresponding MD descriptors actually are of importance. Furthermore, each descriptor set contained significant descriptors, suggesting that each descriptor set contains information that can help discriminate ERK2 actives from inactives. 3.3. Hierarchical Clustering of ERK2 Inhibitors. Next, we used hierarchical clustering to investigate whether active ERK2 ligands clustered together in each descriptor space. Clustering was conducted using an average linkage as the metric between clusters. Tanimoto distances were used for the MACCS fingerprints, and Euclidean distances were used for all other interval valued descriptor sets (see Methods). Ligands were labeled with their active or inactive classes, and then external validation with the ARI metric was used to select the optimal number of clusters for each cluster dendrogram (see Table 2).

Table 2. Measures of Clustering Accuracy for the Optimal Number of Clusters Found by Hierarchical Clustering of ERK2 Inhibitors for Each Descriptor Seta MACCS 2D 3D MD a

ARI

RI

purity

no. of clusters

0.05 0.12 0.25 0.29

0.50 0.54 0.64 0.65

0.78 0.80 0.78 0.78

15 14 4 4

The optimal number of clusters is reported for each descriptor space.

In each dendrogram (Figure 5), we found the optimal clustering to be a large cluster with a high proportion of actives that is distinct from all other clusters with a high proportion of inactives. Table 2 shows the exact measures of how well the optimal clustering grouped together active and inactive ligands. For each descriptor set, the cluster assignments have equally high purity values, indicating that clusters have similar levels of class homogeneity in each dendrogram. However, the MACCS and 2D dendrograms have much lower RI and ARI values compared to those obtained for the 3D and MD dendrograms (because the optimal number of clusters is much larger for the MACCS and 2D dendrograms). This results in a larger number of false positives (compounds from the same activity class in different clusters), which is obviously penalized by the RI and ARI metrics. One can note that MD descriptor-based clustering obtained the highest ARI and RI values, demonstrating that ligands with similar activity are clustering together best in the MD descriptor space. These results indicate that MD descriptors are able to distinguish active and inactive ligands very well and that the information they provide exceeds that of 3D descriptors computed for the top-docked conformation of a ligand only. They also suggest that active ligands clustering together were characterized by highly similar MD descriptor profiles (even though the MD descriptors were only computed for 400 conformations sampled from all generated MD conformations). Moreover, the actual cluster assignments differed considerably between descriptor sets. We quantified the similarity between hierarchical clustering dendrograms using the cophenetic correlation. Table S4 shows the pairwise cophenetic correlation between all pairs of dendrograms. The largest correlation (0.65) was obtained between 3D and MD descriptor dendrograms, which is sensible given that MD descriptors are calculated as summary measures of the distribution of 3D descriptors over MD simulations. There was also a moderate cophenetic correlation (0.62) between the MACCS and 2D descriptor dendrograms. All other pairs of descriptor set dendrograms afforded cophenetic correlations below 0.3. Clustering using different descriptor sets resulted in dendrograms with no more than low/moderate similarity. This supports the conclusion that each descriptor set is capturing different information that is relevant to characterize ligands’ structural information with regard to their ERK2 binding affinity. 3.4. Resolution of Activity Cliffs with MD Descriptors. For each dendrogram corresponding to each descriptor set, the cophenetic correlation coefficient (CCC) was computed as the Pearson correlation of the cophenetic distance of each pair of compounds and their true distances in the corresponding descriptor space (see Methods). The resulting cophenetic correlation coefficients were 0.89 for MACCS, 0.79 for 2D, 0.77 for 3D, and 0.74 for MD. These elevated CCC values demonstrated that the distance between two compounds in each dendrogram is an accurate representation of the true distance 1292


Article


Figure 5. Hierarchical clustering of CSAR and Model set ligands using (A) MACCS, (B) 2D, (C) 3D, and (D) MD descriptors. Clustering was performed using Euclidean distances and average linkage. Ligands are colored by pKi (high: green, low: red). Blue circles demonstrate the cut in each dendrograms, whereas arrows demark the nodes at which CSAR 1 (purple) and CSAR 6 (cyan) merge with the nearest compound with |ΔpKi| > 1.5.

between compounds in the respective whole descriptor space. The fact that the MD descriptor space has the lowest CCC merely indicates that there is slightly more distortion going from true distances to cophenetic distances for this particular descriptor space. These results implied that cluster dendrograms can be reliably used to identify and study activity cliffs in each descriptor space. Moving from the tips of the tree inward, a pair of ligands that merge early on in the cluster dendrogram are indeed similar to each other in that descriptor space. If the pair of ligands have a large difference in pKi, these ligands are forming an activity cliff in that descriptor space. Figure 5 provides two examples of compounds (CSAR1 and CSAR6) that merge early in the MACCS, 2D, and 3D dendrograms with other compounds with |ΔpKi| > 1.5 but merge more deeply in the MD dendrograms tree with compounds with |ΔpKi| > 1.5. The compound CSAR1 (pKi = 4.8) forms an activity cliff with several ligands in MACCS, 2D, and 3D descriptor spaces but does not form an activity cliff with any ligand in the MD descriptor space (Table 3). Figure 6 shows the 3D descriptors

Table 3. Closest Compounds to CSAR1 (pKi = 4.8) with |ΔpKi| > 1.5 in MACCS, 2D, 3D, and MD Descriptor Spacesa MACCS 2D 3D MD

compound

distance

ΔpKi

CSAR 12 CSAR 15 Model 37 Model 26

0.13 4.0 4.3 6.8

3 3.6 2.6 1.8

Tanimoto distances are reported for MACCS fingerprint space, and Euclidean distances are reported for all other descriptor spaces.

a

computed over the entire MD trajectories of CSAR14 (pKi = 8.3) and CSAR1 (pKi = 4.8), whose Tanimoto similarity in MACCS fingerprint descriptor space is 0.83. There are only two slight structural differences between these two ligands: the propane1,3-diol on the terminal amine in CSAR14 (versus no substitution in CSAR1) and the α-ethanol versus α-propanol substitution of the amide group in CSAR14 compared to CSAR1. However, these slight structural modifications translate to a dramatic difference in protein−ligand interactions as predicted by MD simulations. In Figure 6, we represented all the dynamic 1293


Article


Figure 6. Example 1 of an activity cliff present in MACCS, 2D, and 3D descriptor spaces is resolved in the MD descriptor space. Heatmaps of 3D descriptor values for 20,000 conformations sampled over the course of the ligands’ MD trajectories. Cells are colored by descriptor values (0: blue, 1: yellow). Log2(fold change) between CSAR14 and CSAR1 descriptor means and variances above columns (low: red, high: green). Protein ligand interactions that occur in over 30% of the sampled conformations are shown.

with |ΔpKi| > 1.5 in all three MACCS, 2D, and 3D descriptor spaces. However, CSAR6 does not form an activity cliff with any ligand in the MD descriptor space (Table 4). CSAR6 (pKi = 6.1) and CSAR18 (pKi = 9) are structurally similar (Tanimoto similarity = 0.86). The only two differences between these compounds were a substitution of the N-isopropyl group in CSAR6 with an N-methylpyridine in CSAR18 and a methyl group on C3 of a benzene ring in CSAR6. However, there are major differences regarding the dynamic protein−ligand interactions these two ligands create. CSAR18 forms a water bridge with Asp165 via its hydroxyl group, while CSAR6 forms a hydrogen bond with Asn152 with that same hydroxyl group. Meanwhile, CSAR18 and CSAR6 MD descriptor profiles showed large fold changes in means and variances. Interestingly, many of the substantial changes were caused by a subset of MD descriptors with much lower variance - between −4 and −2 log2(fold change) - for the active CSAR18 compound (see Figure 7). 3.5. Distribution of Conformations in MD Space. Principal component analysis was employed to visualize the distribution of conformations that were sampled along the MD trajectories for each ligand (Figure 8 and Figure S5). Exactly 400 conformations were sampled at regular intervals (every 50 ps)

protein−ligand interactions for both compounds (as long as they were present in over 30% of all the 20,000 recorded frames of the MD simulations). Using this criterion, one can observe many differences in the dynamic noncovalent protein−ligand interactions formed by CSAR14 and CSAR1 with ERK2. For example, CSAR14 formed hydrogen bonds between its propanediol group and two ERK2 residues (Asp109 and Lys112), whereas CSAR1 cannot form these intermolecular interactions due to the absence of the propanediol group. Moreover, the lack of the propanediol in CSAR1 also results in two hydrogen bonds with Met106 that are not present in CSAR14. Finally, CSAR14 forms a strong hydrophobic interaction with Tyr34 and a hydrogen bond with Gln103 that are not present in CSAR1. Therefore, it is no surprise that CSAR14 and CSAR1 MD descriptor profiles showed dramatic dissimilarities with many MD descriptors having large differences in means and variances (Figure 6). That is why the two compounds are not activity cliffs in the MD descriptor space (unlike when using fingerprints, 2D, or 3D descriptors). Another interesting example of an activity cliff formed in MACCS fingerprint, 2D, and 3D descriptor spaces that is resolved in the MD descriptor space is shown in Figure 7. The compound CSAR6 (pKi = 6.1) is much closer to several ligands 1294


Article


Figure 7. Example 2 of activity cliffs present in MACCS, 2D, and 3D descriptor spaces are resolved in the MD descriptor space. The same caption as Figure 6.

deviation of the ligand atomic coordinates from their centroid (Ctd.sigma) and the standard deviation of the ligand atomic coordinates from their mediod (Cst.sigma). The one other descriptor was the 3D-Whim descriptor, Van der Waals volumes weighted WV. These three MD descriptors alone (each of which is a measure of global ligand size and shape) contain sufficient information to separate ERK2 actives and inactives. Figure S6 shows three scatterplots obtained for all pairwise combinations of Cst.Sigma, Ctd.Sigma, and Atomic.Masses.WV. Interestingly, Atomic.Masses.WV can be interpreted as an indirect and weighted measure of a ligand’s size, and thus it is not surprising to observe that VdW.Volumes.WV and Atomic.Masses.WV are correlated with r = 0.93. Both Ctd.sigma and Cst.sigma descriptors are more difficult to interpret but provide additional information about ligands’ three-dimensional shape. There is a clear separation between actives and inactives in each descriptor space (Figure S6). Figure S6C shows a particularly good separation between ERK2 actives and inactives. Actives are clearly characterized by larger values for both Atomic.Masses.WV and Cst.Sigma, and the variance for these values is much larger, confirming the highly significant differences between ERK2 actives and inactives in the corresponding MD descriptors shown earlier. More specifically, conformations of most active ligands have Atomic.Masses.WV above 200 and Cst.sigma above 4.

Table 4. Closest Compounds to CSAR6 (pKi = 6.1) with |ΔpKi| > 1.5 in MACCS, 2D, 3D, and MD Descriptor Spacesa MACCS 2D 3D MD

compound

distance

ΔpKi

CSAR 15 Model 52 Model 43 CSAR 13

0.11 5.5 5.7 11

2.3 2.3 1.7 2

Tanimoto distances are reported for MACCS fingerprint space, and Euclidean distances are reported for all other descriptor spaces.

a

along each ligand’s MD trajectory. Again, the full set of 3D descriptors was computed for each conformation, and then PCA was performed on the resulting data. We found a clear separation between the active and inactive ERK2 ligands along the first principal component (PC1). Only a subset of the descriptors has large loadings for PC1 (Figure S2). When analyzing the distribution of the PC1 loadings (Figure S3), we found an apparent separation between the descriptors with an absolute PC1 loading > 0.15 and all other descriptors. Most of these descriptors were highly correlated. If a pair of these descriptors had |r| > 0.9, the descriptor with the largest mean absolute correlation was removed. After filtering, three descriptors with large PC1 loadings remained. Two of these descriptors were 3D-D Moments descriptors: the standard 1295


Article


Figure 8. MD descriptor space for conformations sampled from ligand MD trajectories reveals distinct regions predominately occupied by actives or inactives. 400 conformations were sampled at regular intervals for each ligand’s MD trajectory. The first two principal components from a PCA of the 3D descriptors computed for each conformation. On the z axis is the number of conformations that are present in each bin. A filled contour plot was mapped to the surface. Contour spaces are colored by pKi (high: blue, low: red). Figure 9. PCA of the conformations of (A) CSAR and (B) model set ligands sampled over their MD trajectories. The same PCA as shown in Figure 8, with Model and CSAR sets plotted separately. Convex hulls are drawn around the conformations that belong to each ligand, and their centroid is plotted. The hulls and centroids are colored by their corresponding ligand’s pKi (high: red, low: blue).

A logistic regression model was constructed to quantify more precisely the relationships between these two descriptors and ligands’ binding affinity at ERK2. To assess how well this model performed on held out data, conformations were ranked by the activity of their associated ligand, and every third observation was removed to construct a test set (n = 11,600 conformations). The logistic regression model with Atomic.Masses.WV and Cst.sigma as predictor variables and ligand activity as the response was fit to the training set (n = 23,200 conformations). Then, predicted probabilities of activity were obtained for each test observation. Conformations were predicted to be active if the probability exceeded 0.5. The resulting prediction performance showed this model performed much better than random chance (error rate: 0.25, specificity: 0.78, sensitivity: 0.66). The logistic regression model was then fit to the entire data set. The estimated coefficients for Atomic.Masses.WV and Cst.Sigma were 0.013 and 0.075, respectively. One can actually use this very simplistic MD-QSAR model to design new ERK2 inhibitors with enhanced binding affinity. Regarding the most active ligands, there is one main region on the PCA plot (Figure 8) where there are prominent peaks corresponding to actives. Indeed, active and very active ligands have most of their sampled conformations corresponding to these peaks (Figure 9). This could also provide strong support for the hypothesis of the existence of several possible and metastable binding modes leading to the inhibition of ERK2. Figure 9 further shows how the conformations sampled along ligand MD trajectories are distributed in the MD descriptor space. Clearly, there is much more separation between the active and inactive compounds for the Model set. However, in both Model and CSAR sets, the active compounds (pKi > 7.5) cluster closely together. Still, there are several inactive CSAR ligands, such as CSAR6 (pKi = 6.1), that have their MD conformational distributions that closely resemble active ligands in the model set.

For CSAR6, this confirms the result in the MD cluster dendrogram (Figure 5) that though CSAR6 may not be considered close in MD space to any ligands with large differences in pKi, many of its nearest neighbors having large differences in pKi. In Figure 9, the regions of the plot containing the majority of active compounds in the Model space contain several inactive compounds in the CSAR space. This not only presents a challenge for the prediction of a couple of CSAR set ligands using conformation samples from MD trajectories but also calls for the development of more complex MD descriptors computed on larger samples of conformations.

4. DISCUSSION This study relies on the hypothesis that characterizing dynamic noncovalent interactions of protein−ligand complexes is beneficial for differentiating the ligands with the highest binding affinities. To do so, we have developed a workflow that encompasses (i) the structure-based docking of a series of ligands in the binding site of the target of interest, (ii) the independent MD simulations of each protein−ligand complex, (iii) and the computation of MD descriptors as conformationdependent 3D descriptors of the ligands. We applied this modeling workflow on the challenging set of 87 ERK2 inhibitors. In a previous study9 dedicated to the benchmarking of QSAR and docking prediction performances, we showed how poorly the best models performed when assessing the ranking accuracy for 1296


Article


models and molecular docking alone did clearly not reach this level of segregation. MD descriptors were also shown to be useful in resolving activity cliffs present in all of the other descriptor spaces. Inactive CSAR1 and CSAR6 compounds showed substantial activity cliffs with highly active compounds in MACCS, 2D, and 3D descriptor spaces. Not surprisingly, we have previously shown that both of these compounds had their binding affinities dramatically overpredicted by the best performing QSAR models (ΔpKipred‑exp > 1.5).9 Also, both of these compounds were ranked much higher based on their docking score than by their experimental pKi. Examining these compounds in the MD descriptor space revealed that these ligands have MD trajectories that are not similar to other active compounds. Furthermore, CSAR1 was found to have a MD trajectory that is most similar to compounds of comparable activity. Therefore, ligands that appear to form activity cliffs in canonical descriptor spaces may not form an activity cliff when their MD trajectories are taken into account and examined. PCA of the MD simulation conformational space has also revealed the persistent challenge in using the MD descriptor set to predict the activity of the CSAR set using the Model set. Several inactive CSAR ligands - such as CSAR6 - took conformations that lie within the region predominately occupied by active Model set ligands. The MD cluster dendrogram also showed CSAR6 clustering with active ligands. This indicated that though CSAR6 may not be considered close in MD space to any ligands with large differences in pKi, many of its nearest neighbors have large differences in pKi. Therefore, QSAR models predicting the activity of this ligand may still struggle to accurately predict the activity of this compound. The difference in protein ligand interactions between CSAR6 and CSAR18 suggested that a new set of descriptors based on dynamic protein−ligand interactions may further improve QSAR prediction for this data set. We are currently exploring the definition of these new MD descriptors. Importantly, this project could not have been possible without GPU acceleration enabling the 87 independent, solvent-explicit MD simulations. These simulations were performed using a GPU workstation based on dual Intel Xeon 24-core processors and a 4-way Nvidia Titan X setup (∼12k CUDA cores). As we were able to run simultaneous MD simulations, the runs were completed in ∼8 days, whereas CPU Desmond would require ∼103 CPU days.

several data sets including ERK2. Herein, as shown by the updated results we obtained with Glide, we reconfirmed that current scoring functions were still not able to rank ERK2 inhibitors accurately (especially when it comes to several inactive compounds such as CSAR1 and CSAR6 that are dramatically mispredicted by docking). Our ultimate goal is the development of highly accurate QSAR models for the prediction of the binding affinity toward ERK kinases. In this context, analyzing the ability of each molecular descriptor set to distinguish ERK actives from inactives is critical. The results suggest that each descriptor set contains unique and complementary information to describe ligands’ abilities to bind to ERK2. Moreover, because no descriptor set was found to be strongly correlated with the experimental pKi values, our results underline the fact that each descriptor set might contain descriptors that are relevant, orthogonal, and synergistic in improving the prediction performance of QSAR models. Obviously, variable selection methods are needed to identify those top descriptors from each descriptor set. These observations support the perspective of Tzeng and Hopfinger41 that many different classes of chemical descriptors (e.g., 1D through nD) should be considered when constructing a QSAR model. Using only descriptors from one descriptor class (typically 2D) may result in a large loss of information about ligands’ characteristics. Several 3D descriptors we computed for the docked poses of each ligand seemed to be useless to accomplish the discrimination between ERK2 actives and inactives. However, when the distributions of these 3D descriptors over the course of MD simulations were examined, we have shown that some of these MD descriptors computed as conformation- and timedependent 3D descriptors could fully discriminate the most active from the less active ERK2 inhibitors. These results demonstrate how MD simulations can provide us with an additional layer of key information on ligands’ dynamic characteristics that are of importance for further structural optimization. In particular, we found that ligands with large values of Atomic.Masses.WV (related to ligands’ size) had higher probability of being active. Those size (Atomic.Masses.WV) and shape (Cst.sigma) related descriptor thresholds along with the simplistic MD-QSAR model based on logistic regression may help us to design and identify new analogs likely to be more potent at ERK2. Interestingly, some other MD descriptors computed as the standard deviations of the distributions of 3D descriptors over ligand MD trajectories (e.g., “Atomic.Masses.WV.sd”) were found to be significantly associated with ligands’ affinities to ERK2. For these descriptors, we surprisingly observed that active ligands tended to have larger standard deviations than inactive ligands. More work is needed to understand this result, but one could hypothesize that the most active ligands can adopt conformations that better and dynamically adapt to the flexible kinase pocket. Moreover, the hierarchical clustering results demonstrated that ligands with similar activities have similar MD descriptor profiles, much more so than any other descriptor set we examined. This, along with the results demonstrating a lack of correlation between the MD descriptors and the other descriptor sets, suggests that our MD descriptor set provides new and relevant information about ligands’ activity toward ERK2. In particular, the PCA-based visualization of the MD conformational space for the ERK2 inhibitors has revealed a clear segregation between active and inactive ERK2 ligands. QSAR

5. CONCLUSIONS This study demonstrates the yet untapped potential of exploiting MD trajectories of ensembles of protein−ligand complexes to inform QSAR models via the calculation of specific MD chemical descriptors. Our case study on ERK2 inhibitors shows that such MD-generated descriptors help in discriminating the most potent ERK2 binders, including some wrongly predicted by 2D QSAR models and 3D molecular docking. Therefore, MD descriptors computed for the characterization of conformational ensembles of molecules could boost the prediction performances of QSAR models. Importantly, these results further confirm the early pioneering studies of Hopfinger et al.11 with 4D-QSAR models. However, the astonishing technological advancements of modern computers, especially GPU acceleration, enable the conduction of much longer simulations (tens to hundreds of nanoseconds) for larger sets of protein−ligand complexes. This study paves the way to hyperpredictive MD-QSAR models involving classical 2D/3D molecular descriptors and MD 1297


Article


Using the 4D-QSAR Analysis Formalism. J. Am. Chem. Soc. 1997, 119 (43), 10509−10524. (12) ChEMBL Database. https://www.ebi.ac.uk/chembl/ (accessed Mar 13, 2013). (13) Fourches, D.; Muratov, E.; Tropsha, A. Trust, but Verify: On the Importance of Chemical Structure Curation in Cheminformatics and QSAR Modeling Research. J. Chem. Inf. Model. 2010, 50 (7), 1189− 1204. (14) Madhavi Sastry, G.; Adzhigirey, M.; Day, T.; Annabhimoju, R.; Sherman, W. Protein and Ligand Preparation: Parameters, Protocols, and Influence on Virtual Screening Enrichments. J. Comput.-Aided Mol. Des. 2013, 27 (3), 221−234. (15) Olsson, M. H. M.; Søndergaard, C. R.; Rostkowski, M.; Jensen, J. H. PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa Predictions. J. Chem. Theory Comput. 2011, 7 (2), 525−537. (16) Harder, E.; Damm, W.; Maple, J.; Wu, C.; Reboul, M.; Xiang, J. Y.; Wang, L.; Lupyan, D.; Dahlgren, M. K.; Knight.; et al. A. OPLS3: A Force Field Providing Broad Coverage of Drug-like Small Molecules and Proteins. J. Chem. Theory Comput. 2016, 12 (1), 281−296. (17) Friesner, R. A.; Banks, J. L.; Murphy, R. B.; Halgren, T. A.; Klicic, J. J.; Mainz, D. T.; Repasky, M. P.; Knoll, E. H.; Shelley, M.; Perry, J. K.; et al. Glide: A New Approach for Rapid, Accurate Docking and Scoring. 1. Method and Assessment of Docking Accuracy. J. Med. Chem. 2004, 47 (7), 1739−1749. (18) Friesner, R. A.; Murphy, R. B.; Repasky, M. P.; Frye, L. L.; Greenwood, J. R.; Halgren, T. A.; Sanschagrin, P. C.; Mainz, D. T. Extra Precision Glide: Docking and Scoring Incorporating a Model of Hydrophobic Enclosure for Protein−Ligand Complexes. J. Med. Chem. 2006, 49 (21), 6177−6196. (19) Schrödinger Release 2016-3: LigPrep; Schrödinger, LLC: New York, NY, 2016. (20) Bowers, K. J.; Chow, E.; Xu, H.; Dror, R. O.; Eastwood, M. P.; Gregersen, B. A.; Klepeis, J. L.; Kolossvary, I.; Moraes, M. A.; Sacerdoti, F. D. Scalable Algorithms for Molecular Dynamics Simulations on Commodity Clusters. 2006; No. October, 43. (21) Berthold, M. R.; Cebron, N.; Dill, F.; Gabriel, T. R.; Kötter, T.; Meinl, T.; Ohl, P.; Sieb, C.; Thiel, K.; Wiswedel, B. {KNIME}: The {K}onstanz {I}nformation {M}iner. In Studies in Classification, Data Analysis, and Knowledge Organization (GfKL 2007); Springer: 2007; DOI: 10.1007/978-3-540-78246-9_38. (22) Durant, J. L.; Leland, B. A.; Henry, D. R.; Nourse, J. G. Reoptimization of MDL Keys for Use in Drug Discovery. J. Chem. Inf. Comput. Sci. 2002, 42 (6), 1273−1280. (23) RDKit: Open-source cheminformatics. http://www.rdkit.org (accessed August 20, 2016). (24) Ballester, P. J.; Richards, W. G. Ultrafast Shape Recognition to Search Compound Databases for Similar Molecular Shapes. J. Comput. Chem. 2007, 28 (10), 1711−1723. (25) Todeschini, R.; Gramatica, P. The WHIM Theory: New 3D Molecular Descriptors for QSAR in Environmental Modelling. SAR QSAR Environ. Res. 1997, 7 (1−4), 89−115. (26) R Core Team. R: A Language and Environment for Statistical Computing; Vienna, Austria, 2016. (27) Kuhn, M. Building Predictive Models in R Using the Caret Package. J. Stat. Softw. 2008, 28 (5), 1−26. (28) Kuhn, M.; Wing, J.; Weston, S.; Williams, A.; Keefer, C.; Engelhardt, A.; Cooper, T.; Mayer, Z.; Kenkel, B.; The R Core Team. et al. Classification and Regression Training, R package version 6.0-73; 2016. https://Cran.R-Project.Org/Package=Caret. (29) Benjamini, Y.; Hochberg, Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J. R. Stat. Soc. 1995, 57 (1), 289−300. (30) Cohen, J. Statistical Power Analysis for the Behavioral Sciences; L. Erlbaum Associates: Hillsdale, NJ, 1988. (31) Müllner, D. {fastcluster}: Fast Hierarchical, Agglomerative Clustering Routines for {R} and {Python}. J. Stat. Softw. 2013, 53 (9), 1−18.

molecular descriptors together. These MD-QSAR models would not be intended for screening large libraries of chemicals but could become extremely useful for lead optimization and analogue prioritization, for which very high prediction reliability is required.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.7b00048. Information as mentioned in the text and chemical structures used in this study (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Phone: (+1) 919-515-8711. E-mail: [email protected]. ORCID

Denis Fourches: 0000-0001-5642-8303 Author Contributions

J.A. conducted all the calculations and wrote the manuscript. D.F. conceived the study and wrote the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors gratefully thank the NCSU Chancellor’s Faculty Excellence Program. J.A. would like to thank Dr. Jacqueline Hughes-Oliver (NCSU) for fruitful discussions and the NIEHS for support (training grant T32ES007329).

■

REFERENCES

(1) Zhang, J.; Yang, P. L.; Gray, N. S. Targeting Cancer with Small Molecule Kinase Inhibitors. Nat. Rev. Cancer 2009, 9 (1), 28−39. (2) Gross, S.; Rahal, R.; Stransky, N.; Lengauer, C.; Hoeflich, K. P. Targeting Cancer with Kinase Inhibitors. J. Clin. Invest. 2015, 125 (5), 1780−1789. (3) Wu, P.; Nielsen, T. E.; Clausen, M. H. Small-Molecule Kinase Inhibitors: An Analysis of FDA-Approved Drugs. Drug Discovery Today 2016, 21 (1), 5−10. (4) Hu, Y.; Bajorath, J. Compound Promiscuity: What Can We Learn from Current Data? Drug Discovery Today 2013, 18 (13-14), 644−650. (5) Griffen, E.; Leach, A. G.; Robb, G. R.; Warner, D. J. Matched Molecular Pairs as a Medicinal Chemistry Tool. J. Med. Chem. 2011, 54 (22), 7739−7750. (6) Sebolt-Leopold, J. S.; Herrera, R. Targeting the Mitogen-Activated Protein Kinase Cascade to Treat Cancer. Nat. Rev. Cancer 2004, 4 (12), 937−947. (7) Dhillon, A.; Hagan, S.; Rath, O.; Kolch, W. MAP Kinase Signalling Pathways in Cancer. Oncogene 2007, 26, 3279−3290. (8) Damm-Ganamet, K. L.; Smith, R. D.; Dunbar, J. B.; Stuckey, J. A.; Carlson, H. A. CSAR Benchmark Exercise 2011−2012: Evaluation of Results from Docking and Relative Ranking of Blinded Congeneric Series. J. Chem. Inf. Model. 2013, 53 (8), 1853−1870. (9) Fourches, D.; Muratov, E.; Ding, F.; Dokholyan, N. V.; Tropsha, A. Predicting Binding Affinity of CSAR Ligands Using Both StructureBased and Ligand-Based Approaches. J. Chem. Inf. Model. 2013, 53 (8), 1915−1922. (10) Götz, A. W.; Williamson, M. J.; Xu, D.; Poole, D.; Le Grand, S.; Walker, R. C. Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 1. Generalized Born. J. Chem. Theory Comput. 2012, 8 (5), 1542−1555. (11) Hopfinger, A. J.; Wang, S.; Tokarski, J. S.; Jin, B.; Albuquerque, M.; Madhav, P. J.; Duraiswami, C. Construction of 3D-QSAR Models 1298


Article

Journal of Chemical Information and Modeling (32) Yu, G.; Smith, D. K.; Zhu, H.; Guan, Y.; Lam, T. T.-Y. Ggtree: An R Package for Visualization and Annotation of Phylogenetic Trees with Their Covariates and Other Associated Data. Methods Ecol. Evol. 2017, 8 (1), 28−36. (33) Sokal, R. R.; Rohlf, F. J. The Comparison of Dendrograms by Objective Methods. Taxon 1962, 11 (2), 33−40. (34) Sokal, R. R. PHENETIC TAXONOMY: Theory and Methods. Annu. Rev. Ecol. Syst. 1986, 17 (1), 423−442. (35) Galili, T. Dendextend: An R Package for Visualizing, Adjusting, and Comparing Trees of Hierarchical Clustering. Bioinformatics 2015, 31 (22), 3718−3720. (36) Hubert, L.; Arabie, P. Comparing Partitions. J. Classif. 1985, 2, 193−218. (37) Rand, W. M. Objective Criteria for the Evaluation of Clustering Methods. J. Am. Stat. Assoc. 1971, 66 (66), 846−850. (38) Rendón, E.; Abundez, I.; Arizmendi, A.; Quiroz, E. M. Internal versus External Cluster Validation Indexes. Int. J. Comput. Commun. 2011, 5 (1), 27−34. (39) Vu, V. Q. ggbiplot: A ggplot2 Based Biplot, R package version 0.55; 2011. http://github.com/vqv/ggbiplot (accessed August 11, 2016). (40) Sievert, C.; Parmer, C.; Hocking, T.; Chamberlain, S.; Ram, K.; Corvellec, M.; Despouy, P. Plotly: Create Interactive Web Graphics via “Plotly.js”, R package version 4.5.6; 2016. https://CRAN.R-project.org/ package=plotly (accessed October 27, 2016). (41) Tseng, Y. J.; Hopfinger, A. J.; Esposito, E. X. The Great Descriptor Melting Pot: Mixing Descriptors for the Common Good of QSAR Models. J. Comput.-Aided Mol. Des. 2012, 26 (1), 39−43.

1299


Characterizing the Chemical Space of ERK2 Kinase Inhibitors Using

Recommend Documents