Protein-Ligand Empirical Interaction Components for Virtual Screening

Article pubs.acs.org/jcim

Protein−Ligand Empirical Interaction Components for Virtual Screening Yuna Yan,†,‡,§ Weijun Wang,†,‡,§ Zhaoxi Sun,†,‡,§ John Z. H. Zhang,†,‡,§ and Changge Ji*,†,‡,§ †

Shanghai Engineering Research Center for Molecular Therapeutics and New Drug Development, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai 200062, China ‡ State Key Laboratory of Precision Spectroscopy, East China Normal University, Shanghai 200062, China § NYU-ECNU Center for Computational Chemistry at NYU Shanghai, Shanghai 200062, China S Supporting Information *

ABSTRACT: A major shortcoming of empirical scoring functions is that they often fail to predict binding affinity properly. Removing false positives of docking results is one of the most challenging works in structure-based virtual screening. Postdocking filters, making use of all kinds of experimental structure and activity information, may help in solving the issue. We describe a new method based on detailed protein− ligand interaction decomposition and machine learning. Protein−ligand empirical interaction components (PLEIC) are used as descriptors for support vector machine learning to develop a classification model (PLEIC-SVM) to discriminate false positives from true positives. Experimentally derived activity information is used for model training. An extensive benchmark study on 36 diverse data sets from the DUD-E database has been performed to evaluate the performance of the new method. The results show that the new method performs much better than standard empirical scoring functions in structure-based virtual screening. The trained PLEIC-SVM model is able to capture important interaction patterns between ligand and protein residues for one specific target, which is helpful in discarding false positives in postdocking filtering.

1. INTRODUCTION

functions which are used to predict protein−ligand binding affinity.21−26 Empirical scoring functions are typically trained on a data set of complexes with known affinities. It is difficult to create a universal scoring function that meets all types of targets. A common phenomenon in many VS campaigns is the retrieval of false positives among true positives. Postdocking process is needed for true drug discovery when it comes to screening large chemical libraries against a novel target using the knowledge of known active binding modes. Detailed information on intermolecular interaction between proteins and their ligands is critical for a target specific virtual screening and can be used for constructing target specific scoring methods. Mooij et al. suggest that the target specific scoring methods perform better than the general scoring functions only if a sufficient amount of information is available.27 Several standalone programs (e.g., Silver,28 Post-Dock,29 and Viscana30) have been developed for postdocking process. Knowledge-based weighting for ligand-based similarity search31,32and molecular field analysis of docking poses (AFMoC)33 have used the structural information to create target specific scoring methods. Several methods decompose protein−ligand energy terms into a per-residue basis, including comparative binding

Virtual screening (VS) has made an increasing impact in lead identification and optimization, which make the drug discovery process more efficient.1−7 Compared to experimental screening methods, in silico virtual screening is much faster and more economical and productive. A typical virtual chemical library screening in silico can narrow down the number of ligands to be physically screened. Virtual screening can start either from the three-dimensional structure of a drug target or from known ligands of the target. For the ligand-based virtual screening methods (LBVS),8−12 the main computational technology focuses on developing fast and efficient fingerprints, descriptors, and similarity algorithms. For the structure-based virtual screening methods (SBVS),13−19 most of the scientific efforts are centered on docking algorithms and scoring functions development. A typical structure-based virtual screening process using docking programs usually includes two steps. In the first step, compounds are docked to the putative binding pocket on the premise of shape and interaction complementarity.15,20 In the second step, the binding free energy for a particular conformation of a compound interacting with the target is estimated by a scoring function and those with most favorable energy scores are selected for further biological study. A major challenge in docking programs is the accuracy of the scoring © 2017 American Chemical Society

Received: January 10, 2017 Published: July 5, 2017 1793

DOI: 10.1021/acs.jcim.7b00017 J. Chem. Inf. Model. 2017, 57, 1793−1806

Article

Journal of Chemical Information and Modeling Table 1. Protein Targets for PLEIC-SVM Benchmarking Collected from DUD-E

a

familya

protein

PDB

activesb

decoysc

family

protein

PDB

actives

decoys

protease protease protease protease protease protease protease protease nuclear nuclear nuclear nuclear nuclear nuclear nuclear nuclear kinase kinase

ada17 bace1 dpp4 fa10 hivpr mmp13 thrb try1 andr esr1 esr2 gcr ppara ppard pparg prgr akt1 cdk2

2oi0 3l5d 2i78 3kl6 1xl2 830c 1ype 2ayw 2am9 1sj0 2fsz 3bqd 2p54 2znp 2gtk 3kba 3cqw 1h00

532 283 533 537 536 572 461 449 269 383 367 258 373 240 484 293 293 474

35724 18100 40950 28325 35750 37200 27004 25980 14350 20685 20199 15000 19399 12250 25300 15650 16450 27850

kinase kinase kinase kinase kinase GPCR GPCR GPCR GPCR others others others others others others others others others

egfr lck mk14 src vgfr2 adrb1 adrb2 drd3 aa2ar aces cah2 dhi1 dyr fnta hivrt parp1 pde5a pgh2

2rgp 2of2 2qd9 3el8 2p2i 2vt4 3ny8 3pbl 3eml 1e66 1bcd 3frj 3nxo 3e37 3lan 3l3m 1udt 3ln1

542 420 578 524 409 247 231 480 482 453 492 330 231 592 338 508 398 435

35050 27400 35850 34500 24950 15850 15000 34050 31550 26250 31172 19350 17200 51500 18891 30050 27550 23150

Protein family classification of selected protein targets. bNumber of actives collected from DUD-E. cNumber of decoys collected from DUD-E.

energy analysis (COMBINE)34 and the tailored scoring function TScore.35 Both the COMBINE and TScore methods use decomposed protein−ligand interaction terms from structural data of ligand−macromolecule complexes to train a regression model for ligand activity prediction. Recent studies have focused on the use of interaction fingerprints (IFs) as a postdocking strategy.36−42 In 2003, Deng et al. introduced the concept of the structural interaction fingerprints (SIFt)37 to describe and analyze three-dimensional protein−ligand interactions. After that, many variations of SIFt were developed. In 2009, the Tingjun Hou group developed a method named MIEC-SVM43 using the force-based protein−ligand interaction matrix and SVM to discriminate the binding peptides from the nonbinders for protein domains, and the prediction results have been validated by various experiments.44,45 In 2010, Sato et al. developed pharmacophore-based interaction fingerprint (Pharm-IF)41 and examined its utility for virtual screening using machine learning techniques. In 2014, C. Da and D. Kireev introduced a new approach termed structural protein− ligand interaction fingerprints (SPLIF),36 which also exploits the general idea of quantifying and comparing all possible ligand−protein interactions that may occur. In most of the cases, IFs have been combined with similarity metrics or machine learning approaches to rank or classify compounds for in silico screening. Machine learning has already been widely used in the scope of drug design due to high accuracy of nonlinear fitting, such as target prediction,46 druglikeness,47,48molecular and biological activities,49−59 and ADME/Tox60,61 profiles of small compounds using different types of molecular fingerprints or descriptors. Widely used machine-learning approaches include support vector machines (SVM), decision trees (DT), k-nearest neighbors (k-NN), naive Bayesian methods, artificial neural networks (ANN), etc. Here we introduce an approach termed protein−ligand empirical interaction components (PLEIC) to map protein− ligand binding interactions on protein pocket residues. Typical empirical scoring function terms such as hydrogen bonds and van der Waals and hydrophobic interactions were accounted for in PLEIC. We combine the decomposed protein−ligand empirical interaction elements with machine learning techni-

ques to build a classification model (PLEIC-SVM) as a postdocking process for virtual screening. In order to quantitatively assess the performance of this new approach, 36 diverse protein targets are tested. The targets along with the sets of respective actives and decoys are collected from the Database of Useful Decoys: Enhanced (DUD-E).62 The enrichment performance of PLEIC-SVM is compared to that from Glide docking results. Detailed interaction patterns between ligands and protein are analyzed to illustrate possible reasons for better performance of the new method.

2. METHODS 2.1. Data Sets. The Database of Useful Decoys: Enhanced (DUD-E),62 a standard test set for virtual screening, is used to evaluate the performance of our approach. 36 targets are selected from the DUD-E database. All the selected targets have more than 200 active compounds. The initial number of actives and decoys and the PDBID used for docking and PLEIC generation are listed in Table 1. The 36 target sets cover a wide range of popular drug targets, which contain 4 GPCR, 7 kinases, 8 nuclear receptors, 8 proteases, and 9 other target families. To increase scaffold diversity and to make the decoy set smaller, the DUD-E database has clustered the raw ChEMBL ligands by their Bemis−Murcko atomic frameworks.62 For each target, half of the actives and decoys are randomly selected as the training set for the model construction, while the remaining molecules are used as the external test set for the model validation. 2.2. Data Preparation and Docking. Receptor structures are prepared with the Protein Preparation Wizard (PrepWizard) in Maestro prior to docking. LigPrep is used with default settings to prepare the actives and decoys. Ligands are docked into the active site of the target protein using the Glide program63,64 with standard docking precision (Glide SP). The Default Glide options for grid generation are employed including the use of a 10 Å cubic inner box. The outer box is defined by 30*30*30 centered on the native ligand contained in crystal structure. The top one pose for each ligand scored by Glide SP scoring function is retained for interaction component construction. 1794


Article

Journal of Chemical Information and Modeling

Figure 1. Illustration of PLEIC fingerprint generation. (a) Selecting compounds from libraries and docking compounds to the protein target. (b) Calculating the residue-based interaction. Different colors represent different interaction strengths. (c) Generating the PLEIC fingerprints. Values in the Excel grid are the scaled data, and the different colors represent different strengths.

2.3. Protein−Ligand Empirical Interaction Components (PLEIC). Up to now, a variety of scoring functions have been established. It can be roughly classified into three categories: force field-based, knowledge-based, and empirical scoring functions. The empirical scoring function is widely used because of its speed and accuracy. The empirical scoring functions estimate binding affinity by explicitly accounting for various physicochemical terms that can contribute to the binding free energy, such as van der Waals (VDW) interaction, hydrophobic contact, hydrogen bonding, electrostatic interaction, metal−ligand binding, etc. In this paper, we calculate three kinds of interactions including van der Waals interaction, hydrogen bonding, and hydrophobic contact, which are most widely used in several scoring functions.63−66 Although different ligands bind to a specific protein with different interaction modes, they share a same protein structure. Thus, we can do special training on a specific target protein by mapping each kind of protein−ligand interaction on protein residues in the binding pocket. The PLEIC-SVM model is based on the generation of the interaction matrix. The native ligand located in the complex structure is used as the reference ligand for pocket definition. The pocket is defined by the residues around 7 Å of the reference ligand. We use a large pocket to guarantee that most of the actives and decoys can bind in the pocket. A residue is included when it has hydrophobic contact or hydrogen bond with the active ligands. An interaction element is selected only if its frequency of the actives is more than 5%. Atomic radii for both protein and ligand atoms are set as follows: carbon, 1.9 Å; oxygen, 1.7 Å; nitrogen, 1.8 Å; phosphorus, 2.1 Å; sulfur, 2.0 Å; fluorine, 1.5 Å; chlorine, 1.8 Å; bromine, 2.0 Å; iodine, 2.2 Å; metals, 1.2 Å; hydrogen, 0 Å. Radii of other atoms not included in the above are set at 1.8 Å. The overall process is illustrated in Figure 1. Details of each interaction type are described as follows. 2.3.1. Van der Waals Interaction. The van der Waals (VDW) interaction is one of the most important noncovalent interactions. Here, VDW interaction is calculated using the Lennard-Jones 8-4 potential formula, which is a softer form of standard 12-6 equation.

ligand res

∑ ∑ VDWjk = ∑

VDWi =

j

k

⎡⎛ ⎞8 ⎛ ⎞4 ⎤ ⎢⎜ d0 ⎟ − 2⎜ d0 ⎟ ⎥ ∑ ⎢⎜ ⎟ ⎜d ⎟ ⎥ d ⎝ jk ⎠ ⎦ k ⎣⎝ jk ⎠

ligand res j

(1)

where VDWi denotes the van der Waals interaction energy between residue i and the ligand, which is calculated by considering all the atom pairs between the ligand and the protein residue i; djk denotes the distance between the ligand atom j and the atom k of protein residue i, and d0 denotes the sum of atomic radii of atoms j and k. 2.3.2. Hydrophobic Contact. The hydrophobic effect is calculated by summing up the hydrophobic atom pairs formed between the ligand and the protein residues. The algorithm adopted by ChemScore67,68 is used in our calculations. In our implementation, it is calculated as ligand res

HCi =

∑ ∑ f (djk) j

(2)

k

where f (d) = 1.0

d ≤ d0 + 0.5 Å

= (1/1.5) × (d0 + 2.0 − d) d0 + 0.5 Å < d ≤ d0 + 2.0 Å =0

d > d0 + 2.0 Å

where HCi represents the hydrophobic interaction energy between residue i and the ligand, which is calculated by considering all the hydrophobic contacts between the ligand the protein residue i; d0 is the sum of the atomic radii of atoms j and k. djk denotes the distance between the ligand atom j and atom k of protein residue i. 2.3.3. Hydrogen Bonding. A typical hydrogen bond is defined between a hydrogen bond donor and a hydrogen bond acceptor. The strength of the hydrogen bond is valued by the extent to which the angle between the N−H bond vector and the O···H vector differs from 180° and the distance between donor atom and acceptor atom. We define the distance cutoff at 3.5 Å and the angle cutoff at 120°. Hydrogen bond energy is calculated as 1795


Article


Figure 2. Workflow of the PLEIC-SVM model construction and activity predicting.

Table 2. Detailed Information of Descriptors and Parameters Used in PLEIC-SVM Model Construction systems

resa

desb

Cc

γc

systems

res

des

C

γ

aa2ar aces ada17 adrb1 adrb2 akt1 andr bace1 cah2 cdk2 dhi1 dpp4 drd3 Dyr egfr esr1 esr2 fa10

33 41 36 38 34 32 35 39 23 31 41 33 36 32 42 43 30 33

68 86 74 81 78 68 63 86 54 67 74 66 73 66 87 82 62 70

2.00 2.83 0.71 4.00 8.00 0.71 8.00 1.00 11.31 1.00 8.00 11.31 11.31 8.00 11.31 0.71 8.00 5.66

1.41 0.50 1.41 0.12 0.12 1.41 0.12 0.71 0.35 2.00 0.50 0.25 0.35 1.00 0.18 0.50 0.12 0.25

fnta gcr hivpr hivrt lck mk14 mmp13 parp1 pde5a pgh2 ppara ppard pparg prgr src thrb try1 vgfr2

44 46 44 35 36 43 32 41 36 43 45 45 40 35 52 36 41 48

91 82 96 68 79 82 68 80 71 85 93 92 83 64 96 77 84 88

11.31 4.00 2.00 11.31 5.66 2.00 2.00 4.00 11.31 11.31 5.66 2.83 0.71 11.31 5.66 11.31 8.00 2.00

0.12 0.35 0.25 0.25 0.50 1.00 0.71 0.50 0.12 0.18 0.12 0.12 2.00 0.18 0.25 0.18 0.12 0.50

a

Number of selected pocket residues surrounding the reference ligand. bNumber of descriptors used for PLEIC generation. cParameters used in PLEIC-SVM training. ligand res

HBi =

ligand res

⎛

1 6 ⎝ (1 + (djk /2.6)

∑ ∑ HBjk = ∑ ∑ ⎜⎜ j

k

j

k

⎞ 0.58⎟⎟ ⎠

and has a lower risk of model overfitting. The radial basis function (RBF) kernel transfers sample data into a higher dimensional space, and it is suitable for the case when the relation between the class labels and attributes is nonlinear. There are two parameters for an RBF kernel: C (the penalty parameter of the error term) and γ (the kernel parameter). Grid searching is employed to do cross-validation for the purpose of finding the best parameters. The default evaluating criterion for a model is accuracy, which is the degree of closeness of measurements of a quantity to that quantity’s true value. However, for some unbalanced data sets, accuracy may not be a good criterion. Here we use the AUC (the area under the ROC

(3)

where HBi denotes hydrogen bond energy between residue i and the ligand; djk denotes the distance between ligand atom j and atom k of the protein residue i, where j and k are the heavy atoms of hydrogen bond donor or acceptor. 2.4. PLEIC-SVM Model Construction. Support Vector Machine69 (SVM) is an attractive machine learning method in drug design.70−73 SVM is a classifier based on the structural risk minimization principle, which is less affected by duplicated data 1796


Article


Figure 3. Comparison of performances for different interaction type combinations. Performance of Gscore is also pictured as a comparison. (a) AUC; (b) ROC1.0%. The red and blue histograms denote the average value obtained from the 36 targets, and the black error bars represent the standard deviation separately.

Figure 4. Influence of hyperparameters on performance of PLEIC-SVM model. Grid search and 5-fold cross validation is used to tune the parameters for AKT1 and VGFR2. The black squares represents the default parameter pairs in the software libsvm.

parameters for SVM training on the training set; (vi) using the best parameter pair to generate the final classifier; (vii) validating the model using the external test set and calculating performance metrics; (viii) predicting novel compounds using the trained model. Detailed information on selected targets studied is shown in Table 1. Systems we studied cover most of the popular drug target types, including protease, GPCR, nuclear receptor, kinase, and ion channel. Residues around the binding pocket are all included in constructing the interaction matrix. However, not all selected pocket residues are included in the final PLEIC. Some interactions, which are of rare occurrences or small values, are excluded. Dimensions of PLEIC and numbers of residues in the binding pocket of each target used for SVM training are listed in Table 2. Multiple factors can affect the performance of a prediction model. Here, we examine the influence of different combinations of interaction types and of parameter pairs used in SVM model training. 3.2. Combination of Different Interaction Types. As is shown in Figure 3, six different combinations of interaction types are tested. The performances of different combinations of interaction types are compared. Two common enrichment metrics based on receiver−operator characteristic curves are

curve) criterion to evaluate models. The SVM algorithm implemented in the libsvm74 package is employed for the model construction. After the best parameter pair (C, γ) is found, the whole training set is trained again to generate the final classifier. The C and γ values are designed exponentially growing against 2, namely, 2n (C = 2−3, 2−2.5, 2−2, ..., 25; γ = 2−3, 2−2.5, 2−2, ..., 23). The grid space is set to 0.5 for both C and γ. Different (C, γ) parameters may all lead to the best AUC values. In such a situation, we picked the pair which has the smallest C value. Due to the unbalance of the actives and decoys (∼1:50), a higher weight 5 is set for the actives to balance the classification in the SVM training process.

3. RESULTS AND DISCUSSION 3.1. Overview of the PLEIC-SVM Model. The workflow of the PLEIC-SVM method is shown in Figure 2. The generic workflow for each target includes the following steps: (i) collecting actives and decoys from database (actives and decoys from DUD-E here); (ii) Glide-based docking and scoring of all the actives and decoys; selecting the best poses ranked by Gscore; (iii) calculating the decomposed empirical interactions and using the filter criteria to exclude the rarely occurring components; (iv) constructing the protein−ligand empirical interaction components; (v) using grid searching to tune the 1797


Article

Journal of Chemical Information and Modeling Table 3. Performance Statistics for the Test Sets Using ROC Analysisa protein aa2ar aces ada17 adrb1 adrb2 akt1 andr bace1 cah2 cdk2 dhi1 dpp4 drd3 dyr egfr esr1 esr2 fa10 fnta gcr hivpr hivrt lck mk14 mmp13 parp1 pde5a pgh2 ppara ppard pparg

AUCb

95% CIc

ROC0.5%d

ROC1.0%d

ROC2.0%d

ROC5.0%d

0.86 0.95 0.77 0.93 0.81 0.97 0.84 0.95 0.88 0.92 0.81 0.94 0.84 0.90 0.84 0.91 0.56 0.88 0.91 0.91 0.69 0.88 0.84 0.93 0.72 0.92 0.85 0.97 0.79 0.93 0.90 0.97 0.93 0.98 0.95 0.98 0.75 0.87 0.79 0.91 0.79 0.90 0.85 0.89 0.91 0.93 0.72 0.91 0.75 0.96 0.96 0.99 0.87 0.89 0.83 0.90 0.82 0.92 0.61 0.86 0.84

0.84−0.88 0.93−0.97 0.73−0.81 0.91−0.95 0.78−0.84 0.96−0.98 0.80−0.88 0.92−0.97 0.84−0.92 0.89−0.96 0.77−0.85 0.91−0.96 0.80−0.89 0.86−0.95 0.80−0.87 0.88−0.94 0.52−0.59 0.85−0.91 0.89−0.93 0.88−0.93 0.65−0.73 0.85−0.91 0.82−0.86 0.92−0.95 0.69−0.75 0.90−0.94 0.81−0.89 0.95−0.99 0.77−0.82 0.90−0.95 0.87−0.94 0.95−0.99 0.91−0.96 0.96−0.99 0.94−0.97 0.97−0.99 0.72−0.78 0.85−0.90 0.74−0.84 0.88−0.94 0.76−0.82 0.88−0.93 0.82−0.88 0.85−0.92 0.89−0.93 0.91−0.95 0.69−0.76 0.89−0.93 0.73−0.78 0.94−0.97 0.95−0.97 0.98−0.99 0.84−0.90 0.86−0.92 0.79−0.86 0.87−0.93 0.80−0.85 0.90−0.95 0.55−0.67 0.80−0.91 0.82−0.87

0.10 0.64 0.23 0.55 0.33 0.78 0.17 0.63 0.30 0.60 0.12 0.56 0.17 0.46 0.10 0.58 0.01 0.39 0.34 0.57 0.08 0.30 0.08 0.55 0.01 0.50 0.18 0.74 0.23 0.66 0.60 0.79 0.60 0.80 0.57 0.81 0.06 0.31 0.20 0.43 0.18 0.53 0.16 0.38 0.26 0.59 0.21 0.60 0.07 0.58 0.41 0.89 0.17 0.32 0.28 0.59 0.05 0.51 0.01 0.38 0.08

0.17 0.68 0.27 0.65 0.36 0.82 0.22 0.67 0.41 0.66 0.16 0.64 0.22 0.51 0.18 0.63 0.03 0.48 0.43 0.63 0.12 0.39 0.14 0.61 0.02 0.55 0.26 0.78 0.27 0.71 0.66 0.85 0.67 0.84 0.64 0.86 0.10 0.36 0.23 0.47 0.22 0.57 0.24 0.43 0.36 0.63 0.27 0.64 0.11 0.67 0.56 0.90 0.20 0.41 0.34 0.62 0.08 0.59 0.02 0.48 0.19

0.28 0.73 0.32 0.69 0.39 0.85 0.28 0.72 0.50 0.76 0.20 0.71 0.33 0.57 0.28 0.66 0.05 0.55 0.55 0.68 0.15 0.47 0.23 0.69 0.05 0.64 0.34 0.81 0.32 0.74 0.72 0.86 0.73 0.87 0.74 0.90 0.14 0.45 0.31 0.54 0.28 0.63 0.35 0.56 0.50 0.66 0.33 0.71 0.16 0.76 0.69 0.92 0.30 0.49 0.43 0.66 0.14 0.66 0.04 0.54 0.31

0.48 0.82 0.40 0.76 0.47 0.91 0.45 0.81 0.64 0.81 0.37 0.78 0.48 0.73 0.42 0.70 0.08 0.64 0.66 0.73 0.20 0.59 0.39 0.78 0.11 0.72 0.49 0.88 0.41 0.79 0.80 0.89 0.82 0.93 0.84 0.93 0.23 0.56 0.39 0.62 0.39 0.73 0.51 0.66 0.64 0.75 0.40 0.76 0.26 0.84 0.81 0.96 0.45 0.60 0.51 0.74 0.30 0.77 0.12 0.63 0.46

1798


Article

Journal of Chemical Information and Modeling Table 3. continued protein prgr src thrb try1 vgfr2 ave SDf

AUCb

95% CIc

ROC0.5%d

ROC1.0%d

ROC2.0%d

ROC5.0%d

0.92 0.92 0.96 0.84 0.93 0.91 0.95 0.90 0.95 0.77 0.95 0.82 0.93 0.09 0.03

0.90−0.95 0.89−0.94 0.94−0.98 0.81−0.86 0.91−0.95 0.89−0.93 0.94−0.97 0.87−0.92 0.93−0.97 0.73−0.81 0.93−0.96 0.79−0.85 0.90−0.95 0.09−0.08 0.04−0.03

0.60 0.23 0.65 0.13 0.58 0.46 0.71 0.41 0.80 0.21 0.55 0.22 0.58 0.16 0.15

0.67 0.37 0.68 0.19 0.62 0.55 0.73 0.49 0.82 0.26 0.62 0.28 0.64 0.18 0.14

0.72 0.52 0.74 0.28 0.73 0.61 0.77 0.60 0.84 0.33 0.70 0.35 0.69 0.19 0.12

0.77 0.61 0.82 0.42 0.79 0.69 0.84 0.68 0.87 0.42 0.78 0.47 0.77 0.19 0.10

a

For every protein target, the line above is the performance for Gscore, and the line below is the performance for PLEIC-SVM. bArea under the ROC curve. c95% confidence interval for AUC. dTrue positive rates at three “early” false positive rate values (0.5%, 1%, 2%, and 5%). eThe mean values of different metrics. fStandard deviation of the mean.

Figure 5. ROC curves for 9 DUD-E targets. Red and blue real lines represent the ROC curves for PLEIC-SVM and Gscore separately.

1799


Article


calculate the confidence interval. As is shown in Table 3, the AUC values using our model are all higher than or equal to the corresponding results derived from Glide SP scoring. In general, the PLEIC-SVM model has demonstrated a higher average AUC value at 0.93 versus Gscore at 0.82. An AUC value of about 0.8 indicates that there are a substantial number of decoys outscoring the known DUD actives. Robust performance is indicated by an AUC score that is in excess of 0.9. We plotted the ROC curves for nine targets in Figure 5. The average values of ROC0.5%, ROC1.0%, ROC2.0%, and ROC5.0% increase from 0.22, 0.28, 0.35, and 0.47 to 0.58, 0.64, 0.69, and 0.77 by using PLEIC-SVM model. A scatter plot of ROC1.0% values gained from PLEIC-SVM versus ROC1.0% values gained from Glide scoring is shown in Figure S1. These results indicate that there is more opportunity to find potential hits with much lower false positive rate by using PLEIC-SVM. From Table 3, we can find that the average value of AUC gained from Glide SP scoring is 0.82. In some targets, the AUC values obtained from Gscore are also very high. Among 36 targets, the AUC values of 7 targets (cdk2, esr2, fa10, lck, parp1, prgr, thrb) exceed 0.90. When using PLEIC-SVM, the overall performance gets a systematic improvement. For example, the AUC value of target HIVPR increased from 0.79 to 0.90, target GCR increased from 0.79 to 0.91, target PGH2 increased from 0.83 to 0.91, and target PPARG increased from 0.84 to 0.92. However, in some cases, improvements are not so obvious, for example, CDK2, ANDR, ADRB2, LCK, HIVRT, and PDE5A. In these cases, the 95% confidence interval for AUC values obtained from Gscore and PLEIC-SVM overlapped. Nevertheless, from Table 3, we can find that the performance of other metrics for these systems, including ROC0.5%, ROC1.0%, ROC2.0%, and ROC5.0%, gets an improvement. For example, for target LCK, the AUC values gained from Glide scoring and PLEICSVM are 0.91 and 0.93. The 95% confidence interval of two AUCs gained from Glide scoring and PLEIC-SVM overlapped, which are 0.89−0.93 and 0.91−0.95. The ROC0.5%, ROC1.0%, ROC2.0%, and ROC5.0% for target LCK increased from 0.26, 0.36, 0.50, and 0.64 to 0.59, 0.63, 0.66, and 0.75 by using PLEIC-SVM. These results suggest that PLEIC-SVM is a useful tool for ranking of actives and decoys. We also calculate ROC early enrichments for the test set of every target at several early false positive rates (ROC enrichment @ 0.5%, 1.0%, 2.0%, 5%, and 10% FPR) as suggested by Jain and Nicholls.76 Enrichment at 1% (ROC [email protected]%FPR) is the fraction of actives seen along with the top 1% of known decoys (multiplied by 100). This metric removes the dependence on the ratio of actives and inactives and directly quantifies early enrichment. In Table S3, EF@X%FPR for the test set of 36 targets and the average values for different target families are listed. The average values of [email protected]%FPR, [email protected]%FPR, [email protected]%FPR, [email protected]%FPR, and [email protected]% FPR for all 36 targets increased from 43.49, 27.81, 17.73, 9.35, and 5.79 to 116.21, 63.44, 34.71, 15.38, and 8.27 by using PLEIC-SVM. From Table S3, we can see that the enrichments obtained from PLEIC-SVM at all early false positive rates for each target are better than those from Gscore. Several previous studies also used ROC enrichment metric in their work to assess the performance of different docking softwares and scoring functions. We list direct comparison of ROC enrichments with two previous studies in Tables S4 and S5. In 2009, Cross et al. evaluated the performance of several molecular docking programs (DOCK, FlexX, Glide, ICM,

employed to evaluate the performance, including the area under the ROC curve (AUC) and the value on the y-axis at the 1.0% points of the x-axis of the ROC plot (ROC1.0%). ROC1.0% is the true positive rate when the false positive rate is at 1.0%. The average values and the standard deviations of AUCs and ROC1.0% obtained from the test sets of the 36 targets are calculated and plotted in Figure 3. As is shown in Figure 3, all kinds of combinations, except the combination of single interaction type HB, perform better than Glide scoring. The combination of all three energy components performs best. Both the ROC1.0% and AUC performance can be ranked as (HB) < Gscore < (VDW) < (HYD) < (VDW + HB) < (HYD + HB) < (VDW + HYD + HB). The model using only van der Waals or hydrophobic interaction components can also get a better performance than Gscore, partly due to the fact that most of the binding affinity comes from hydrophobic matching between the ligand and protein. 3.3. Hyperparameters Used in SVM. As has been discussed above, two hyperparameters (C and γ) may significantly affect the prediction performance of the model when using the RBF kernel function. Just using the default parameters (such as in libsvm C = 1.0 and γ = 1/n features) may miss the best model. By using the grid searching approach, C and γ values are optimized by finding the grid with the highest AUC value. As an illustration of the influence of different parameters on performance, we selected two targets (ATK1 and VGFR2) using all the data and tested a wide range of C and γ. For target AKT1, 68 interaction descriptors are used in SVM model training, and the default log2 γ is −6.1 and log2 C is 0. Similarly, for target VGFR2, the default log2 γ is −6.5 and log2 C is 0. As is shown in Figure 4, default parameters do not perform best in either of the two systems, which are highlighted by black squares. It is essential to tune the parameters for the model construction. Considering the time cost and avoiding overfitting, we selected a narrower range of C and γ in the current study. The optimized hyperparameters for each system are listed in Table 2. 3.4. Performance of the PLEIC-SVM Model. In order to quantitatively assess the performance of this approach, we use four kinds of enrichment metrics based on receiver operating characteristic, which are widely used computational evaluations of docking methods.75,76 Receiver operating characteristic (ROC), or ROC curve, is a graphical plot that illustrates the performance of a binary classifier system as its discrimination threshold is varied. We quantitatively compared the values of ROC0.5%, ROC1.0%, ROC2.0%, and ROC5.0%. ROC1.0% is the value on the y-axis at the 1.0% points of the x-axis of the ROC plot. The value of ROC1.0% = 0.25 indicates that 25% of the known actives are recovered by the time 1.0% of the known decoy compounds are screened. Besides, the area under the curve (AUC) is also used as an evaluator which is equal to the probability that a classifier will rank a randomly chosen positive instance higher than a randomly chosen negative one (assuming that “positive” ranks higher than “negative”). To give a comparison, the AUC values, the corresponding 95% confidence interval, ROC0.5%, ROC1.0%, ROC2.0%, and ROC5.0% of the test sets obtained from Gscore and the PLEIC-SVM model using all three interaction types are calculated for all 36 targets. There are several methods to estimate the confidence intervals (95%) for AUC values.77,78 Here we use the “pROC” package79 implemented in R to 1800


Article

Journal of Chemical Information and Modeling PhDOCK, and Surflex) with 40 DUD protein targets.80 We compare performance of overlapped target classes with results from Cross et al. in Table S4. From Table S4, we can find that PLEIC-SVM perform better than traditional methods in every target class. In 2010, Brooijmans et al. evaluated early ROC enrichments for 8 kinases at several early false positive rates.81 We compare performance of overlapped kinases with results from Brooijmans et al. in Table S5. For all of the four kinases, PLEIC-SVM performs better than the results of Brooijmans in each early enrichment. In 2015, Gaudreault and Najmanovch evaluated performance of FlexAID in virtual screening.82 The averaged ROC enrichment at 1% FPR from their study is 2.9. The averaged ROC EF@1% of PLEIC-SVM is 63.4. All these results indicate that PLEIC-SVM performs better than traditional screening methods. PLEIC-SVM is a postdocking method using active and inactive information from a specific target for model building. In traditional docking methods, the scoring functions were built for all targets. Since target specific information is included in our postdocking process, PLEICSVM is robust in capturing important features for better virtual screening. The Boltzmann-enhanced discrimination of ROC (BEDROC) is also used to evaluate the PLEIC-SVM model. BEDROC score is a good metric that is adapted to the “early recognition” problem. The BEDROC score is calculated as n

BEDROC =

∑i = 1 e−αri / N

(

Rα +

1 − e−α eα / N − 1

)

×

AUC values and other enrichment metrics indicate that our method is effective in enriching active molecules in virtual screening. 3.5. Important Interaction Patterns for Model Construction. Distinguishable features may exist between positive and negative samples to construct an effective classification model. In most cases, it is difficult to explain the classification model with one single feature or a simple linear combination of several features. A more complicated explanation is needed for discrimination. Here we present six examples (EGFR, VGFR2, AA2AR, TRY1, ESR1, and LCK) to get a deeper understanding of why the new method performs better by analyzing the importance of different features. 3.5.1. EGFR. The epidermal growth factor receptor (EGFR) is the cell-surface receptor for members of the epidermal growth factor family of extracellular protein ligands. The identification of EGFR as an oncogene has led to the development of anticancer therapeutics directed against EGFR. For this target, ROC1.0% gained from Gscore and PLEIC-SVM is at 0.27 and 0.71, respectively. The area under the ROC curve increases from 0.79 to 0.93 by using PLEICSVM. We use the feature selection tool implemented in libsvm to get the importance of each feature used in PLEIC-SVM model construction. The top 10 important interaction pairs extracted from the PLEIC-SVM model of EGFR are listed in Table S1. From Table S1, we can find that hydrogen bonds between residue MET973 and small molecules rank first, which implies that this interaction pattern is very important for discriminating actives and decoys in the PLEIC-SVM model. Then we analyze all crystal structures of EGFR in the PDB database which are structured with small molecule in the same pocket as the PDB structure 2RGP. There are 76 structures in total, and all of the 76 structures form hydrogen bonds between the small molecule and residue MET793. We superimposed all of the 76 crystal structures to PDB 2RGP, as is shown in Figure 6b. Statistical

R α sinh(α /2) cosh(α /2) − cosh(α /2 − αR α)

1 1 − eα(1 − R α)

(4)

where n is the number of actives, N is the total number of actives and decoys, Rα is the ratio of active in all active and decoy, and ri is the rank of the ith active molecule. BEDROC can be interpreted as the probability that a ranked active randomly selected will be positioned before a randomly selected compound distributed following an exponential of parameter α. The parameter α embodies the degree of “early recognition”. Here we calculate BEDROC with α at 20, 80.5, and 160.9. The details of the calculated BEDROC scores for the test sets using Gscore and PLEIC-SVM are listed in Table S6. As shown in Table S6, BEDROC scores also indicate that the early enrichment of actives improves by using the target specific classification model PLEIC-SVM. Several other molecular docking programs have been tested on DUD-E or DUD. In 2012, Repasky et al. evaluated the Glide SP performance on 40 targets in the DUD data set. The average AUC value for these targets is 0.80.83 In 2013, Sastry et al. evaluated the influence of parameters and protocols on virtual screening enrichments with GLIDE program using 36 targets in the DUD database.84 They found that more complete preparation does produce better virtual screening enrichments. In 2016, Chaput et al. compared the performance of four docking programs, GOLD, GLIDE, SURFLEX, and FlexX, using the 102 targets in the DUD-E database.85 By measuring the BEDROC scores with α = 80.5, they concluded that, on the overall database, GLIDE outperformed for 30 targets, GOLD for 27, FlexX for 14, and Surflex for 11. These studies indicate that Glide SP performs relatively well in docking studies. In our current study, we use Glide as a first step for the postdocking process to improve the enrichment of actives. Detailed comparisons of PIELC-SVM’s performance with previous virtual screening methods are listed in Tables S7 and S8.

Figure 6. The binding site of EGFR. (a) Superposition of a high Gscore decoy molecule with the reference ligand. The active reference liagnd is colored green. The decoy molecule is colored yellow. Hbonds between MET793 and the reference ligand are pictured by the green dotted line. (b) Superposition of all of the EGFR crystal structure ligands on structure 2RGP.

information derived from crystal structures of EGFR complexes also indicates that this hydrogen bond is of great importance for binding. These results suggest that important interaction patterns detected by the PLEIC-SVM procedure are reasonable and are helpful for discriminating false positives in structurebased virtual screening. One decoy molecule of EGFR named ZINC01894479 is picked for detailed analysis. Molecule ZINC01894479 gets a Gscore of −11.03 kcal/mol, which is ranked very high in Glide docking. In the PLEIC-SVM method, ZINC01894479 is recognized as an inactive molecule. We superimposed 1801


Article


insomnia, pain, depression, drug addiction, and Parkinson’s disease. For this target, ROC1.0% gained from Gscore and PLEIC-SVM is 0.17 and 0.68, respectively. The area under the ROC curve increases from 0.86 to 0.95 by using PLEIC-SVM. The binding site of target AA2AR is shown in Figure 8(a). The small molecule shown in green is the native ligand in

ZINC01894479 and the reference active molecule in the binding pocket and pictured them in Figure 6a. As is shown in Figure 6a, there are two H-bonds between MET793 and the reference ligand. However, there is no hydrogen bond between ZINC01894479 and MET793. Although ZINC01894479 gets a very high predicted binding affinity from Gscore, it lacks critical binding elements which may make it an inactive ligand. 3.5.2. VGFR2. Vascular endothelial growth factor receptor-2 (VGFR2) is a high-affinity receptor for vascular endothelial growth factor-A (VEGF-A) and mediates most of the endothelial growth and survival signals from VEGF-A. It is also an important anticancer target. For this target, ROC1.0% gained from Gscore and PLEIC-SVM is at 0.26 and 0.62, respectively. The area under the ROC curve increases from 0.77 to 0.95 by using PLEIC-SVM. The same analysis procedure is performed for VGFR2 to get the most important features. The top 10 important interaction pairs extracted from the PLEIC-SVM model for VGFR2 are listed in Table S2. From Table S2, we can find that hydrogen bonding between the small molecules and residue CYS919 is recognized as the most important feature. A total of 29 crystal structures in the PDB database have the same pocket as the PDB structure 2P2I which is used for model construction. Among them, 28 structures have hydrogen bonds between small molecules and the residue CYS919. We superimposed all of the 29 crystal structures on PDB 2P2I, which is shown in Figure 7b, and the important hydrogen bond is shown by a

Figure 8. The binding site of AA2AR. (a) Superposition of a high Gscore decoy molecule with the reference ligand. The active reference ligand is colored green. The decoy molecule is colored yellow. H-bond between ASN253 and the reference ligand is pictured by the green dotted line. (b) Superposition of all of the AA2AR crystal structure ligands on structure 3EML.

crystal structure 3EML which is used for the classification model construction for target AA2AR. The native ligand has two aromatic rings which are tightly packed in the hydrophobic pocket formed by residues PHE168, MET270, ILE274, and ASN253, and there is a hydrogen bond between the planar nitrogen atoms attached to the aromatic rings of the native ligand and residue ASN253, which is favorable to the binding free energy. A total of 26 crystal structures can be found in the PDB database with a ligand occupying the same pocket as 3EML. All structures form hydrogen bonds with residue ASN253. This suggests that this hydrogen bond is essential for binding. We superimposed all of the 26 crystal structures on PDB 3EML, which is shown in Figure 8b, and the important hydrogen bond is shown by a green dotted line. To check whether our model catches this important hydrogen bond, we analyzed the frequency of hydrogen bonds along with the top N molecules ranked by the active probabilities predicted by the PLEIC-SVM model. As shown in Figure S2, the higher the ranking for one molecule, the higher the chance of the occurrence of the hydrogen bond. For example, 95% of the top 100 molecules form hydrogen bonds with residue ASN253. This indicates that this feature is important for the classification model construction and helpful to the better ranking of actives and decoys. A molecule named ZINC39024671 is picked and is shown in sticks and colored yellow in Figure 8a. This molecule is a decoy molecule, scored at −10.39 kcal/mol by Gscore and predicted to be inactive by the PLEIC-SVM model. A part of this molecule occupies the hydrophobic pocket formed by residues PHE168, MET270, ILE274, and ASN253. The absence of the hydrogen bond between this molecule and residue ASN253 may decrease the probability of this molecule to be active. 3.5.4. TRY1. Trypsin is a serine protease from the PA clan superfamily. It can hydrolyze proteins, which is used for numerous biotechnological processes. In order to prevent the action of active trypsin in the pancreas, which can be highly damaging, many inhibitors are developed. For this target, ROC1.0% gained from Gscore and PLEIC-SVM is 0.49 and 0.82, respectively. The area under the ROC curve increases from 0.90 to 0.95 by using PLEIC-SVM. The binding site of target TRY1 is shown in Figure 9a. The small molecule shown in green is the native ligand in crystal

Figure 7. The binding site of VGFR2. (a) Superposition of a high Gscore decoy molecule with the reference ligand. The active reference liagnd is colored green. The decoy molecule is colored yellow. H-bond between CYS919 and the reference ligand is pictured by the green dotted line. (b) Superposition of all of the VGFR2 crystal structure ligands on structure 2P2I.

green dotted line. Both PLEIC-SVM and crystal structure analysis suggest that this hydrogen is critical for the binding of small molecules. One decoy molecule of VGFR2 named ZINC63710099 is picked for detailed analysis. Molecule ZINC63710099 gets a Gscore of −12.45 kcal/mol, which is also ranked very high in Glide docking results. In the PLEIC-SVM method, ZINC63710099 is recognized as an inactive molecule. We superimposed ZINC63710099 and the reference active molecule in the binding pocket and pictured them in Figure 7a. As is shown in Figure 7a, there is a hydrogen bond between CYS919 and the reference ligand, which is absent for ZINC63710099. On the opposite, there is a carbon atom in ZINC63710099 at the place where an acceptor atom is needed. The lacking of the important hydrogen bond may be the reason why ZINC63710099 is recognized as inactive. 3.5.3. AA2AR. The adenosine A2A receptor is a member of the G protein-coupled receptor (GPCR) family. The adenosine A2A receptor plays a critical role in regulating myocardial oxygen consumption and coronary blood flow. It is also a potential therapeutic target for treatment of conditions such as 1802


Article


through experimental binding data and 3D structure data. Such methods can only do add-on and cannot do subtract. From medicinal chemists’ experience, we know that reduction of a single functional group from the ligand may lead to total loss of activity in certain circumstances, which cannot be captured by current add-only empirical scoring functions. The add-only character may be responsible for high false positive rate of widely used empirical scoring functions. For the PLEIC-SVM method, it captures important recognition patterns between active ligands and the target through training of active and inactive data sets. Significant recognition patterns may include directed hydrogen bonds or hydrophobic packing. Lack of certain binding elements may lead to inactive. Structural analysis from SVM training results may also provide some useful information for lead optimization. It is helpful for pointing out which interactions should be kept and which one can be modified in drug design. There are several points that may contribute to the improved performance of the PLEIC-SVM model. First, the model we build is target specific, while general docking scoring functions are fitted from many kinds of targets. It is difficult for a general scoring function to cover the balance of all kinds of special interactions for all protein targets. In some binding pockets, coupling of hydrogen bond and hydrophobic interaction can make a special contribution to the binding. Some scoring functions have considered some specific contributions by experts’ knowledge. Target specific scoring functions or classification models are considered to be more attractive for dealing with special contributions. Second, our model is based on the known binding activities for one target. The automatic machine learning process helps us make full use of knowledge of available activity information from experiments and uncovering critical interaction patterns for specific protein− ligand recognition. Third, we decomposed three kinds of protein−ligand interactions and mapped them on pocket residues. Detailed interactions mapped on each residue form a set of fingerprints, which is suitable for machine learning. Finally, we use the RBF kernel function in machine learning which transfers sample data into a higher dimensional space and it is suitable for the case when the relation between the class labels and attributes is nonlinear. Compared to target specific scoring functions, the classification model for specific target can make use of more experimentally derived information. Both active and inactive information derived from HTS can be used for model construction. For target specific scoring functions, only active compounds are used.

Figure 9. The binding site of TRY1. (a) Superposition of a high score decoy molecule with the reference ligand. The active reference ligand is colored green. The decoy molecule is colored yellow. Salt bridge between ASP189 and the reference ligand is pictured by the green dotted line. (b) Superposition of all of the TRY1 crystal structure ligands on structure 2AYW.

structure 2AYW which is used for the classification model construction. As shown in Figure 9a, the binding site is very shallow. In the relatively deeper pocket, the native ligand forms a salt bridge like an anchor with residue ASP189. A total of 298 crystal structures in the PDB database bind with the same pocket as the structure 2AYW. 293 ligands form a salt bridge with residue ASN253. This indicates that this salt bridge is important for binding. We superimposed all of the 298 crystal structures on PDB 2AYW, which is shown in Figure 9b. The important salt bridge is shown by green dotted lines. Like the analysis of AA2AR, we analyze the frequency of the hydrogen bonds (salt bridge) along with the top N molecules ranked by the probabilities predicted by the PLEIC-SVM model. As is shown in Figure S3, the higher the ranking for one molecule, the higher the chance of the occurrence of this salt bridge. For example, the top 162 molecules all form hydrogen bonds with residue ASP189, and 83% of the top 300 molecules form hydrogen bonds with residue ASP189. This suggests that this feature contributes a lot to the classification model and the better ranking of actives and decoys. A molecule, named ZINC16697102, is picked and is shown in sticks and colored yellow in Figure 9a. This molecule is a decoy molecule, scored at −9.17 kcal/mol by Gscore and predicted to be inactive by PLEIC-SVM. As shown in Figure 9a, a part of this molecule aligns well with the native ligand, but there is no donor group or positively charged group near the carboxyl group of ASP189. The absence of the salt bridge between this molecule and residue ASP189 may decrease the probability of this molecule to be active. 3.5.5. ESR1 and LCK. Statistical analysis of each interaction pattern is informative for uncovering the distinguishable features. We analyze the average van der Waals hydrophobic and hydrogen bond interaction between each pocket residue and ligand for ESR1 and LCK. Figures S4 and S5 show the relative difference between active and decoy for each interaction feature in ESR1 and LCK. From Figures S4 and S5, we can find that the there is a large gap between active and decoy in several interaction patterns. In ESR1, three of the most distinguishable features are the hydrogen bonds formed with residue GLU353, ARG394, and HIS524. In LCK, the hydrogen bond formed with residue MET319 shows significant differences between active and decoy molecules. Those distinguishable features are very helpful in discriminating actives from decoys. Traditional empirical scoring functions estimate binding affinity through summing up all potential binding elements without judging which one is critical and which one is trivial. The contribution of each binding element is usually trained

4. CONCLUSIONS In this paper, a generally applicable structure-based target specific discrimination model for postdocking processing which uses knowledge of the binding interaction of the active and decoy compounds is developed. PLEIC-SVM is a residue-based interaction decomposition method combined with machine learning algorithms (SVM). Overall, this method demonstrates significant improvements compared to the GlideScore. This study shows that PLEIC-SVM is a useful tool for ranking and filtering virtual screening docking results. The PLEIC-SVM method can capture critical interaction patterns for specific target through training of information from active and inactive compounds and use it to judge whether a new molecule is active, while traditional empirical scoring functions that use add-only rules are weak in discriminating false positives in virtual screening. To some extent, the accuracy of our model 1803


Article


(6) Venhorst, J.; Núñez, S.; Terpstra, J. W.; Kruse, C. G. Assessment of Scaffold Hopping Efficiency by Use of Molecular Interaction Fingerprints. J. Med. Chem. 2008, 51, 3222−3229. (7) Yan, C.; Liu, D.; Li, L.; Wempe, M. F.; Guin, S.; Khanna, M.; Meier, J.; Hoffman, B.; Owens, C.; Wysoczynski, C. L.; et al. Discovery and characterization of small molecules that target the Ral GTPase. Nature 2014, 515, 443−447. (8) Cheeseright, T. J.; Mackey, M. D.; Melville, J. L.; Vinter, J. G. FieldScreen: Virtual Screening Using Molecular Fields. Application to the DUD Data Set. J. Chem. Inf. Model. 2008, 48, 2108−2117. (9) Liu, X.; Jiang, H.; Li, H. SHAFTS: A Hybrid Approach for 3D Molecular Similarity Calculation. 1. Method and Assessment of Virtual Screening. J. Chem. Inf. Model. 2011, 51, 2372−2385. (10) Sastry, G. M.; Inakollu, V. S. S.; Sherman, W. Boosting Virtual Screening Enrichments with Data Fusion: Coalescing Hits from TwoDimensional Fingerprints, Shape, and Docking. J. Chem. Inf. Model. 2013, 53, 1531−1542. (11) Hou, X.; Du, J.; Fang, H.; Li, M. 3D-QSAR Study on a Series of Bcl-2 Protein Inhibitors Using Comparative Molecular Field Analysis. Protein Pept. Lett. 2011, 18, 440−449. (12) Xue, X.; Wei, J.-L.; Xu, L.-L.; Xi, M.-Y.; Xu, X.-L.; Liu, F.; Guo, X.-K.; Wang, L.; Zhang, X.-J.; Zhang, M.-Y.; et al. Effective Screening Strategy Using Ensembled Pharmacophore Models Combined with Cascade Docking: Application to p53-MDM2 Interaction Inhibitors. J. Chem. Inf. Model. 2013, 53, 2715−2729. (13) Desjarlais, R. L.; Seibel, G. L.; Kuntz, I. D.; Furth, P. S.; Alvarez, J. C.; Ortiz de Montellano, P. R.; DeCamp, D. L.; Babe, L. M.; Craik, C. S. STRUCTURE-BASED DESIGN OF NONPEPTIDE INHIBITORS SPECIFIC FOR THE HUMAN IMMUNODEFICIENCY VIRUS-1 PROTEASE. Proc. Natl. Acad. Sci. U. S. A. 1990, 87, 6644− 6648. (14) Grüneberg, S.; Stubbs, M. T.; Klebe, G. Successful Virtual Screening for Novel Inhibitors of Human Carbonic Anhydrase: Strategy and Experimental Confirmation. J. Med. Chem. 2002, 45, 3588−3602. (15) Kitchen, D. B.; Decornez, H.; Furr, J. R.; Bajorath, J. Docking and scoring in virtual screening for drug discovery: methods and applications. Nat. Rev. Drug Discovery 2004, 3, 935−949. (16) Li, S.; Gao, J. M.; Satoh, T.; Friedman, T. M.; Edling, A. E.; Koch, U.; Choksi, S.; Han, X. B.; Korngold, R.; Huang, Z. W. A computer screening approach to immunoglobulin superfamily structures and interactions: Discovery of small non-peptidic CD4 inhibitors as novel immunotherapeutics. Proc. Natl. Acad. Sci. U. S. A. 1997, 94, 73−78. (17) Powers, R. A.; Morandi, F.; Shoichet, B. K. Structure-based discovery of a novel, noncovalent inhibitor of AmpC beta-lactamase. Structure 2002, 10, 1013−1023. (18) Shoichet, B. K.; Stroud, R. M.; Santi, D. V.; Kuntz, I. D.; Perry, K. M. STRUCTURE-BASED DISCOVERY OF INHIBITORS OF THYMIDYLATE SYNTHASE. Science 1993, 259, 1445−1450. (19) Song, H.; Wang, R. X.; Wang, S. M.; Lin, J. A low-molecularweight compound discovered through virtual database screening inhibits Stat3 function in breast cancer cells. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 4700−4705. (20) Kuntz, I. D.; Blaney, J. M.; Oatley, S. J.; Langridge, R.; Ferrin, T. E. A geometric approach to macromolecule-ligand interactions. J. Mol. Biol. 1982, 161, 269−288. (21) Krovat, E. M.; Langer, T. Impact of Scoring Functions on Enrichment in Docking-Based Virtual Screening: An Application Study on Renin Inhibitors. J. Chem. Inf. Model. 2004, 44, 1123−1129. (22) Leach, A. R.; Shoichet, B. K.; Peishoff, C. E. Prediction of Protein−Ligand Interactions. Docking and Scoring: Successes and Gaps. J. Med. Chem. 2006, 49, 5851−5855. (23) Moitessier, N.; Englebienne, P.; Lee, D.; Lawandi, J.; Corbeil, C. R. Towards the development of universal, fast and highly accurate docking/scoring methods: a long way to go. Br. J. Pharmacol. 2008, 153, S7−S26.

depends on the docking position generated by Glide. In some cases, the docking position may not be good enough for interaction analysis. One limitation of the current study is that we did not use our model to predict new compounds for specific targets and test them by experiments. Combination with experimental studies can help us improve the model.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.7b00017. Important interactions, BEDROC and ROC1.0% studies, H-bond frequencies, and average values of interactions for actives and decoys in target ESR1 and LCK (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Zhaoxi Sun: 0000-0001-8311-7797 John Z. H. Zhang: 0000-0003-4612-1863 Changge Ji: 0000-0001-5477-4894 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank the National Natural Science Foundation of China (Grants 21003048 and 21433004), National Key Research and Development Plan (2016YFA0501700), Shanghai Natural Science Foundation (Grant 14ZR1411900), and Shanghai Putuo District (Grant 2014-A-02) for financial support. C.J. is also supported by the Fundamental Research Funds for the Central Universities and Open Research Fund of the State Key Laboratory of Precision Spectroscopy, East China Normal University. C.J. acknowledges the support of the NYU-ECNU Center for Computational Chemistry at NYU Shanghai. We also thank the Supercomputer Center of ECNU for providing us with computational time.

■

ABBREVIATIONS IF, interaction fingerprint; DUD-E, Database of Useful Decoys: Enhanced; SVM, support vector machine; ROC, receiver operating characteristics

■

REFERENCES

(1) Keseru, G. M.; Makara, G. M. The influence of lead discovery strategies on the properties of drug candidates. Nat. Rev. Drug Discovery 2009, 8, 203−212. (2) Liu, J.; Wang, K.; Xu, F.; Tang, Z.; Zheng, W.; Zhang, J.; Li, C.; Yu, T.; You, X. Synthesis and photovoltaic performances of donor−π− acceptor dyes utilizing 1,3,5-triazine as π spacers. Tetrahedron Lett. 2011, 52, 6492−6496. (3) Liu, X.; Xie, H.; Luo, C.; Tong, L.; Wang, Y.; Peng, T.; Ding, J.; Jiang, H.; Li, H. Discovery and SAR of Thiazolidine-2,4-dione Analogues as Insulin-like Growth Factor-1 Receptor (IGF-1R) Inhibitors via Hierarchical Virtual Screening. J. Med. Chem. 2010, 53, 2661−2665. (4) Nitsche, C.; Klein, C. D. Aqueous microwave-assisted one-pot synthesis of N-substituted rhodanines. Tetrahedron Lett. 2012, 53, 5197−5201. (5) Ripphausen, P.; Nisius, B.; Peltason, L.; Bajorath, J. Quo Vadis, Virtual Screening? A Comprehensive Survey of Prospective Applications. J. Med. Chem. 2010, 53, 8461−8467. 1804


Article

Journal of Chemical Information and Modeling (24) Sousa, S. F.; Fernandes, P. A.; Ramos, M. J. Protein−ligand docking: Current status and future challenges. Proteins: Struct., Funct., Genet. 2006, 65, 15−26. (25) Wang, R.; Lu, Y.; Wang, S. Comparative Evaluation of 11 Scoring Functions for Molecular Docking. J. Med. Chem. 2003, 46, 2287−2303. (26) Warren, G. L.; Andrews, C. W.; Capelli, A.-M.; Clarke, B.; LaLonde, J.; Lambert, M. H.; Lindvall, M.; Nevins, N.; Semus, S. F.; Senger, S.; et al. A Critical Assessment of Docking Programs and Scoring Functions. J. Med. Chem. 2006, 49, 5912−5931. (27) Mooij, W. T. M.; Verdonk, M. L. General and targeted statistical potentials for protein−ligand interactions. Proteins: Struct., Funct., Genet. 2005, 61, 272−287. (28) Kennard, O.; Watson, D. G.; Town, W. G. Cambridge Crystallographic Data Centre. I. Bibliographic File. J. Chem. Doc. 1972, 12, 14. (29) Springer, C.; Adalsteinsson, H.; Young, M. M.; Kegelmeyer, P. W.; Roe, D. C. PostDOCK: A Structural, Empirical Approach to Scoring Protein Ligand Complexes. J. Med. Chem. 2005, 48, 6821− 6831. (30) Amari, S.; Aizawa, M.; Zhang, J.; Fukuzawa, K.; Mochizuki, Y.; Iwasawa, Y.; Nakata, K.; Chuman, H.; Nakano, T. VISCANA: Visualized Cluster Analysis of Protein−Ligand Interaction Based on the ab Initio Fragment Molecular Orbital Method for Virtual Ligand Screening. J. Chem. Inf. Model. 2006, 46, 221−230. (31) Crisman, T. J.; Sisay, M. T.; Bajorath, J. Ligand-Target Interaction-Based Weighting of Substructures for Virtual Screening. J. Chem. Inf. Model. 2008, 48, 1955−1964. (32) Stiefl, N.; Zaliani, A. A Knowledge-Based Weighting Approach to Ligand-Based Virtual Screening. J. Chem. Inf. Model. 2006, 46, 587− 596. (33) Gohlke, H.; Klebe, G. DrugScore Meets CoMFA: Adaptation of Fields for Molecular Comparison (AFMoC) or How to Tailor Knowledge-Based Pair-Potentials to a Particular Protein. J. Med. Chem. 2002, 45, 4153−4170. (34) Ortiz, A. R.; Pisabarro, M. T.; Gago, F.; Wade, R. C. Prediction of Drug Binding Affinities by Comparative Binding Energy Analysis. J. Med. Chem. 1995, 38, 2681−2691. (35) Matter, H.; Will, D. W.; Nazaré, M.; Schreuder, H.; Laux, V.; Wehner, V. Structural requirements for factor Xa inhibition by 3oxybenzamides with neutral P1 substituents: combining X-ray crystallography, 3D-QSAR, and tailored scoring functions. J. Med. Chem. 2005, 48, 3290−3312. (36) Da, C.; Kireev, D. Structural protein-ligand interaction fingerprints (SPLIF) for structure-based virtual screening: method and benchmark study. J. Chem. Inf. Model. 2014, 54, 2555−61. (37) Deng, Z.; Chuaqui, C.; Singh, J. Structural Interaction Fingerprint (SIFt): A Novel Method for Analyzing Three-Dimensional Protein−Ligand Binding Interactions. J. Med. Chem. 2004, 47, 337− 344. (38) Deng, Z.; Chuaqui, C.; Singh, J. Knowledge-Based Design of Target-Focused Libraries Using Protein−Ligand Interaction Constraints. J. Med. Chem. 2006, 49, 490−500. (39) Marcou, G.; Rognan, D. Optimizing fragment and scaffold docking by use of molecular interaction fingerprints. J. Chem. Inf. Model. 2007, 47, 195−207. (40) Mpamhanga, C. P.; Chen, B.; McLay, I. M.; Willett, P. Knowledge-Based Interaction Fingerprint Scoring: A Simple Method for Improving the Effectiveness of Fast Scoring Functions. J. Chem. Inf. Model. 2006, 46, 686−698. (41) Sato, T.; Honma, T.; Yokoyama, S. Combining Machine Learning and Pharmacophore-Based Interaction Fingerprint for in Silico Screening. J. Chem. Inf. Model. 2010, 50, 170−185. (42) Pérez-Nueno, V. I.; Rabal, O.; Borrell, J. I.; Teixidó, J. APIF: A New Interaction Fingerprint Based on Atom Pairs and Its Application to Virtual Screening. J. Chem. Inf. Model. 2009, 49, 1245−1260. (43) Hou, T.; Xu, Z.; Zhang, W.; McLaughlin, W. A.; Case, D. A.; Xu, Y.; Wang, W. Characterization of Domain-Peptide Interaction

Interface: A Generic Structure-based Model to Decipher the Binding Specificity of SH3 Domains. Mol. Cell. Proteomics 2009, 8, 639−649. (44) Hou, T.; Li, N.; Li, Y.; Wang, W. Characterization of Domain− Peptide Interaction Interface: Prediction of SH3 Domain-Mediated Protein−Protein Interaction Network in Yeast by Generic StructureBased Models. J. Proteome Res. 2012, 11, 2982−2995. (45) Sun, H.; Pan, P.; Tian, S.; Xu, L.; Kong, X.; Li, Y.; Li, D.; Hou, T. Constructing and Validating High-Performance MIEC-SVM Models in Virtual Screening for Kinases: A Better Way for Actives Discovery. Sci. Rep. 2016, 6, 24817. (46) Beveridge, D. L.; Dicapua, F. M. FREE-ENERGY VIA MOLECULAR SIMULATION - APPLICATIONS TO CHEMICAL AND BIOMOLECULAR SYSTEMS. Annu. Rev. Biophys. Biophys. Chem. 1989, 18, 431−492. (47) Müller, K.-R.; Rätsch, G.; Sonnenburg, S.; Mika, S.; Grimm, M.; Heinrich, N. Classifying ‘Drug-likeness’ with Kernel-Based Learning Methods. J. Chem. Inf. Model. 2005, 45, 249−253. (48) Tian, S.; Wang, J.; Li, Y.; Xu, X.; Hou, T. Drug-likeness Analysis of Traditional Chinese Medicines: Prediction of Drug-likeness Using Machine Learning Approaches. Mol. Pharmaceutics 2012, 9, 2875− 2886. (49) Agrafiotis, D. K.; Cedeño, W.; Lobanov, V. S. On the Use of Neural Network Ensembles in QSAR and QSPR. J. Chem. Inf. Model. 2002, 42, 903−911. (50) Byvatov, E.; Fechner, U.; Sadowski, J.; Schneider, G. Comparison of Support Vector Machine and Artificial Neural Network Systems for Drug/Nondrug Classification. J. Chem. Inf. Model. 2003, 43, 1882−1889. (51) Chen, B.; Harrison, R. F.; Papadatos, G.; Willett, P.; Wood, D. J.; Lewell, X. Q.; Greenidge, P.; Stiefl, N. Evaluation of machinelearning methods for ligand-based virtual screening. J. Comput.-Aided Mol. Des. 2007, 21, 53−62. (52) Ehrman, T. M.; Barlow, D. J.; Hylands, P. J. Virtual Screening of Chinese Herbs with Random Forest. J. Chem. Inf. Model. 2007, 47, 264−278. (53) Guha, R.; Jurs, P. C. Interpreting Computational Neural Network QSAR Models: A Measure of Descriptor Importance. J. Chem. Inf. Model. 2005, 45, 800−806. (54) Kauffman, G. W.; Jurs, P. C. QSAR and k-Nearest Neighbor Classification Analysis of Selective Cyclooxygenase-2 Inhibitors Using Topologically-Based Numerical Descriptors. J. Chem. Inf. Model. 2001, 41, 1553−1560. (55) Plewczynski, D.; Spieser, S. A. H.; Koch, U. Assessing Different Classification Methods for Virtual Screening. J. Chem. Inf. Model. 2006, 46, 1098−1106. (56) Sato, T.; Matsuo, Y.; Honma, T.; Yokoyama, S. In Silico Functional Profiling of Small Molecules and Its Applications. J. Med. Chem. 2008, 51, 7705−7716. (57) Schneider, G.; Wrede, P. Artificial neural networks for computer-based molecular design. Prog. Biophys. Mol. Biol. 1998, 70, 175−222. (58) Svetnik, V.; Liaw, A.; Tong, C.; Culberson, J. C.; Sheridan, R. P.; Feuston, B. P. Random Forest: A Classification and Regression Tool for Compound Classification and QSAR Modeling. J. Chem. Inf. Model. 2003, 43, 1947−1958. (59) Winkler, D. A. Neural networks as robust tools in drug lead discovery and development. Mol. Biotechnol. 2004, 27, 139−167. (60) Lei, T.; Li, Y.; Song, Y.; Li, D.; Sun, H.; Hou, T. ADMET evaluation in drug discovery: 15. Accurate prediction of rat oral acute toxicity using relevance vector machine and consensus modeling. J. Cheminf. 2016, 8, 6. (61) Mitchell, J. B. O. Machine learning methods in chemoinformatics. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2014, 4, 468−481. (62) Huang, N.; Shoichet, B. K.; Irwin, J. J. Benchmarking Sets for Molecular Docking. J. Med. Chem. 2006, 49, 6789−6801. (63) 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 1805


Article

Journal of Chemical Information and Modeling and Scoring. 1. Method and Assessment of Docking Accuracy. J. Med. Chem. 2004, 47, 1739−1749. (64) Halgren, T. A.; Murphy, R. B.; Friesner, R. A.; Beard, H. S.; Frye, L. L.; Pollard, W. T.; Banks, J. L. Glide: A New Approach for Rapid, Accurate Docking and Scoring. 2. Enrichment Factors in Database Screening. J. Med. Chem. 2004, 47, 1750−1759. (65) 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, 6177−6196. (66) Wang, R.; Lai, L.; Wang, S. Further development and validation of empirical scoring functions for structure-based binding affinity prediction. J. Comput.-Aided Mol. Des. 2002, 16, 11−26. (67) Eldridge, M. D.; Murray, C. W.; Auton, T. R.; Paolini, G. V.; Mee, R. P. Empirical scoring functions: I. The development of a fast empirical scoring function to estimate the binding affinity of ligands in receptor complexes. J. Comput.-Aided Mol. Des. 1997, 11, 425−445. (68) Murray, C. W.; Auton, T. R.; Eldridge, M. D. Empirical scoring functions. II. The testing of an empirical scoring function for the prediction of ligand-receptor binding affinities and the use of Bayesian regression to improve the quality of the model. J. Comput.-Aided Mol. Des. 1998, 12, 503−519. (69) Cortes, C.; Vapnik, V. Support-vector networks. Mach. Learn. 1995, 20, 273−297. (70) Burbidge, R.; Trotter, M.; Buxton, B.; Holden, S. Drug design by machine learning: support vector machines for pharmaceutical data analysis. Comput. Chem. 2001, 26, 5−14. (71) Doucet, J. P.; Barbault, F.; Xia, H.; Panaye, A.; Fan, B. Nonlinear SVM Approaches to QSPR/QSAR Studies and Drug Design. Curr. Comput.-Aided Drug Des. 2007, 3, 263−289. (72) Warmuth, M. K.; Liao, J.; Rätsch, G.; Mathieson, M.; Putta, S.; Lemmen, C. Active Learning with Support Vector Machines in the Drug Discovery Process. J. Chem. Inf. Model. 2003, 43, 667−673. (73) Zernov, V. V.; Balakin, K. V.; Ivaschenko, A. A.; Savchuk, N. P.; Pletnev, I. V. Drug Discovery Using Support Vector Machines. The Case Studies of Drug-likeness, Agrochemical-likeness, and Enzyme Inhibition Predictions. J. Chem. Inf. Model. 2003, 43, 2048−2056. (74) Chang, C. C.; Lin, C. J. LIBSVM: A library for support vector machines. ACM Trans. Intell. Syst. Technol. 2011, 2, 389−396. (75) Jain, A. N. Bias, reporting, and sharing: computational evaluations of docking methods. J. Comput.-Aided Mol. Des. 2008, 22, 201−212. (76) Jain, A. N.; Nicholls, A. Recommendations for evaluation of computational methods. J. Comput.-Aided Mol. Des. 2008, 22, 133− 139. (77) Brown, C. D.; Davis, H. T. Receiver operating characteristics curves and related decision measures: A tutorial. Chemom. Intell. Lab. Syst. 2006, 80, 24−38. (78) Nicholls, A. Confidence limits, error bars and method comparison in molecular modeling. Part 1: The calculation of confidence intervals. J. Comput.-Aided Mol. Des. 2014, 28, 887−918. (79) Robin, X.; Turck, N.; Hainard, A.; Tiberti, N.; Lisacek, F.; Sanchez, J.-C.; Müller, M. pROC: an open-source package for R and S + to analyze and compare ROC curves. BMC Bioinf. 2011, 12, 77. (80) Cross, J. B.; Thompson, D. C.; Rai, B. K.; Baber, J. C.; Fan, K. Y.; Hu, Y.; Humblet, C. Comparison of Several Molecular Docking Programs: Pose Prediction and Virtual Screening Accuracy. J. Chem. Inf. Model. 2009, 49, 1455−74. (81) Brooijmans, N.; Cross, J. B.; Humblet, C. Biased retrieval of chemical series in receptor-based virtual screening. J. Comput.-Aided Mol. Des. 2010, 24, 1053−1062. (82) Gaudreault, F.; Najmanovich, R. J. FlexAID: Revisiting Docking on Non-Native-Complex Structures. J. Chem. Inf. Model. 2015, 55, 1323. (83) Repasky, M. P.; Murphy, R. B.; Banks, J. L.; Greenwood, J. R.; Tubert-Brohman, I.; Bhat, S.; Friesner, R. A. Docking performance of the glide program as evaluated on the Astex and DUD datasets: a complete set of glide SP results and selected results for a new scoring

function integrating WaterMap and glide. J. Comput.-Aided Mol. Des. 2012, 26, 787−799. (84) Sastry, G. M.; 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, 221−234. (85) Chaput, L.; Martinez-Sanz, J.; Saettel, N.; Mouawad, L. Benchmark of four popular virtual screening programs: construction of the active/decoy dataset remains a major determinant of measured performance. J. Cheminf. 2016, 8, 56.

1806


Protein-Ligand Empirical Interaction Components for Virtual Screening

Recommend Documents